这一部分课程 正式开始学习 单细胞转录组测序的数据分析，大家无需担心自己的 R语言 方面的编程基础。本次课程将从分析工具包的 安装 、 加载 开始，对脚本的参数进行 全方位 的讲解。具体分析内容方面，此次课程将详细的给大家介绍 Seurat对象 的创建、过滤、数据结构特征及理解、存储、线粒体含量计算、分群、降维、Marker计算、各类可视化方法。并同时提醒大家不同物种处理时需要进行一些改动。
此次课程为整个《 单细胞转录组测序数据分析 》课程的 核心内容，其中代码的运行方式、数据结构的理解与处理均十分重要，后续的进阶分析都需以此课程作为基础。关于此讲内容，Seurat V4及V5之间并无明显差别。代码框中内容可以直接复制至 Rstudio 中运行 1. 视频及测试文件 1.1 视频教程 1.2 测试文件 代码中包含测试文件下载方式，如果网络不佳可在推送中获取。
下载本手册的同时本地也包含配套测试文件~ 请在推送中获取:手把手教你做单细胞测序数据分析
（三）:单样本分析 2. 代码 2.1 R包安装与加载 if ( ! require (multtest))BiocManager :: install ( "multtest" ) if ( ! require (Seurat)) install.packages ( "Seurat" ) if ( ! require (dplyr)) install.packages ( "dplyr" ) if ( ! require (patchwork)) install.packages ( "patchwork" ) if ( ! require (R.utils)) install.packages ( "R.utils" ) # gc 函数用于触发垃圾回收（Garbage Collection, GC），回收不再使用的内存，并以矩阵形式返回当前的内存使用情况，其中包含两种类型的内存信息：# Ncells（表示 R 语言中的对象，如列表、函数等） # Vcells（表示向量数据，如数值向量、字符向量等） gc # 设置并行计算，提升下面流程的计算速度 makecore <- function (workcore,memory){ if ( ! require (Seurat)) install.packages ( 'Seurat' ) if ( ! require (future)) install.packages ( 'future' ) plan ( "multicore" , workers = workcore) options ( future.globals.maxSize= memory * 1024 * 1024 ** 2 ) } 核心,允许future任务占用 最多2GB内存为例 makecore ( 4 , 2 ) 2.2 数据读入并创建 Seurat对象 2.2.1 下载并解压数据 # 下载速度受到网络稳定性影响，如果下载失败可以多尝试几次 download.file ( 'https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz' , 'my.pbmc.tar.gz' ) # 解压".tar.gz"结尾的文件完成之后在当前目录即data目录下面会出现一个filtered_gene_bc_matrices文件夹 untar ( "my.pbmc.tar.gz" ) 2.2.2 读入数据并创建Seurat分析对象 pbmc.data <- Read10X ( data.dir = "./filtered_gene_bc_matrices/hg19/" ) 如果你无法正常读入，请尝试理解数据的格式，以下几个教程可供参考：为什么Read10X也会报错? 如何一次性读取一个目录下的cellranger输出文件？
没有barcode文件的单细胞数据要怎么读取 四个单细胞样本只给了一套文件怎么读 pbmc <- CreateSeuratObject ( counts = pbmc.data, project = "pbmc3k" , min.cells = 3 , min.features = 200 ) # 初步查看Seurat对象：pbmc # An object of class Seurat # 13714 features across 2700 samples within 1 assay # Active assay: RNA (13714 features, 0 variable features) # 1 layer present: counts # 查看对象中包含多少个细胞：ncol (pbmc) # [1] 2700 ncol (pbmc.data) # [1] 2700 这里需要注意，Seurat V5中引入了 layer 的概念：所以下面这样无法取出 counts 的稀疏矩阵：# 以下代码会报错：try ({ lalala <- as.data.frame (pbmc[[ "RNA" ]] @ counts) }) # Error in h(simpleError(msg, call)) : # error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': no slot of name "counts" for this object of class "Assay5" 现在的做法应该是：# 这是一种方法：lalala <- pbmc @ assays[[ "RNA" ]] @ layers[[ "counts" ]] # 或者试用GetAssayData函数也是可以的：lalala <- GetAssayData (pbmc, assay= "RNA" , slot= 'counts' ) # 查看一下矩阵的部分数据：lalala[ 1 : 4 , 1 : 4 ] # 4 x 4 sparse Matrix of class "dgCMatrix" # AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 # AL627309.1 . . . # AP006222.2 . . . # RP11-206L10.2 . . . # RP11-206L10.9 . . . # AAACCGTGCTTCCG-1 # AL627309.1 . # AP006222.2 . # RP11-206L10.2 . # RP11-206L10.9 . write.table (lalala, './mycount.txt' , sep = ' \t ' ) # 删除对象，减少内存占用 rm (lalala) 2.3 质控数据及可视化 # if(!require(patchwork))install.packages("patchwork") # 线粒体基因表达比例统计 # 人源的数据为MT，鼠源的需要换成mt pbmc[[ "percent.mt" ]] <- PercentageFeatureSet (pbmc, pattern = "^MT-" ) # 新增percent.mt列 head (pbmc @ meta.data, 5 ) # orig.ident nCount_RNA nFeature_RNA percent.mt # AAACATACAACCAC-1 pbmc3k 2419 779 3.0177759 # AAACATTGAGCTAC-1 pbmc3k 4903 1352 3.7935958 # AAACATTGATCAGC-1 pbmc3k 3147 1129 0.8897363 # AAACCGTGCTTCCG-1 pbmc3k 2639 960 1.7430845 # AAACCGTGTATGCG-1 pbmc3k 980 521 1.2244898 VlnPlot (pbmc, features = c ( "nFeature_RNA" , "nCount_RNA" , "percent.mt" ), ncol = 3 ) # FeatureScatter 是 Seurat 包中的一个函数，它用于绘制散点图（Scatter Plot），可视化两个指标之间的关系 plot1 <- FeatureScatter (pbmc, feature1 = "nCount_RNA" , feature2 = "percent.mt" ) plot2 <- FeatureScatter (pbmc, feature1 = "nCount_RNA" , feature2 = "nFeature_RNA" ) # CombinePlots这步需要你的绘图窗口足够大 CombinePlots ( plots = list (plot1, plot2)) # percent.mt 随 nCount_RNA 增加而增加，可能有一些高 percent.mt 的异常点（低质量细胞） # nCount_RNA 和 nFeature_RNA 成正比，如果 nCount_RNA 很高但 nFeature_RNA 低，可能是低质量细胞或技术噪声 # 注意：percent.mt不是越低越好，部分组织本身代谢高，例如心脏percent.mt相对其他组织数值更高，属于正常现象 # 过滤低质量细胞 pbmc # An object of class Seurat # 13714 features across 2700 samples within 1 assay # Active assay: RNA (13714 features, 0 variable features) # 1 layer present: counts pbmc <- subset (pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5 ) # 看看现在的细胞数：pbmc # An object of class Seurat # 13714 features across 2638 samples within 1 assay # Active assay: RNA (13714 features, 0 variable features) # 1 layer present: counts 2.4 细胞分群 pbmc <- NormalizeData (pbmc, normalization.method = "LogNormalize" , scale.factor = 10000 ) # 寻找高变基因，默认nfeatures=2000 pbmc <- FindVariableFeatures (pbmc, selection.method = "vst" , nfeatures = 2000 ) # 查看十个高变基因 top10 <- head ( VariableFeatures (pbmc), 10 ) top10 # [1] "PPBP" "LYZ" "S100A9" "IGLL5" "GNLY" "FTL" "PF4" "FTH1" # [9] "GNG11" "S100A8" # 可视化单细胞 RNA-seq 数据中的高变基因，x轴基因表达均值；
y轴基因标准变异度。大部分基因点为灰色（背景基因）。高变基因（HVGs）被标记红色。这些高变基因在后续降维（PCA、UMAP）和聚类分析中至关重要。
plot1 <- VariableFeaturePlot (pbmc) # 在 VariableFeaturePlot 上 标注前 10 个高变基因，repel = TRUE避免标签重复，plot2 <- LabelPoints ( plot = plot1, points = top10, repel = TRUE ) plot1 + plot2 pbmc <- ScaleData (pbmc, features = rownames (pbmc)) dim (pbmc @ assays $ RNA @ layers $ scale.data) # [1] 13714 2638 <- ScaleData(pbmc) 默认用高变基因来scale a <- ScaleData (pbmc) dim (a @ assays $ RNA @ layers $ scale.data) # [1] 2000 2638 rm (a) # 提取高变基因 hvg <- VariableFeatures ( object = pbmc) length (hvg) # [1] 2000 # PCA(Principal Component Analysis)主成分分析,主成分分析是一种无监督的降维技术，它能将原始的高维数据转化为一组新的、互不相关的综合变量，即主成分。
这些主成分按照方差大小排序，方差越大，包含的原始数据信息就越多。通过保留前面几个方差较大的主成分，就能在尽量保留原始数据信息的前提下，显著降低数据的维度。单细胞转录组数据通常具有高维度的特点，每个细胞可能有数千个基因的表达值。高维度数据不仅会增加计算的复杂度，还可能导致过拟合等问题。PCA 可以将高维的基因表达数据投影到低维的主成分空间，从而减少数据的维度，降低计算成本。
pbmc <- RunPCA (pbmc, features = hvg) # 输出每个前五个主成分top5基因 print (pbmc[[ "pca" ]], dims = 1 : 5 , nfeatures = 5 ) # PC_ 1 # Positive: CST3, TYROBP, LST1, AIF1, FTL # Negative: MALAT1, LTB, IL32, IL7R, CD2 # PC_ 2 # Positive: CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DQB1 # Negative: NKG7, PRF1, CST7, GZMB, GZMA # PC_ 3 # Positive: HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPB1 # Negative: PPBP, PF4, SDPR, SPARC, GNG11 # PC_ 4 # Positive: HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1 # Negative: VIM, IL7R, S100A6, IL32, S100A8 # PC_ 5 # Positive: GZMB, NKG7, S100A8, FGFBP2, GNLY # Negative: LTB, IL7R, CKB, VIM, MS4A7 # 展示 PCA 降维过程中，每个基因对指定主成分（PC1和PC2）的贡献程度，即基因的载荷值。
基因载荷反映了基因表达值与主成分之间的相关性，通过可视化载荷值，可以帮助我们理解每个主成分所代表的生物学意义。x轴：表示基因在对应主成分上的载荷值。载荷值的绝对值越大，说明该基因对该主成分的贡献越大；
正载荷值表示基因表达与主成分呈正相关，负载荷值表示基因表达与主成分呈负相关 VizDimLoadings (pbmc, dims = 1 : 2 , reduction = "pca" ) # 基于PCA计算结果展示细胞降维分布 DimPlot (pbmc, reduction = "pca" ) # 热图展示指定主成分（PC）与细胞之间的关系。
热图可以帮助我们直观地了解每个主成分在不同细胞中的表达模式，进而识别出对特定主成分贡献较大的细胞 # dims指定要展示的主成分范围，这里表示展示第 1 到第 15 个主成分 # 指定要在热图中展示的细胞数量，这里随机选取 500 个细胞进行展示 # 表示会平衡选取不同簇的细胞，以确保每个簇的细胞在热图中都有适当的代表 DimHeatmap (pbmc, dims = 1 : 15 , cells = 500 , balanced = TRUE ) # JackStraw分析评估主成分显著性，通过对数据进行多次重采样"这里是100次"，模拟出随机数据的主成分分布，然后将真实数据的主成分与随机数据的主成分进行比较，从而判断每个主成分是否具有统计学意义 # 这一步计算时间较长,约4min pbmc <- JackStraw (pbmc, num.replicate = 100 ) # ScoreJackStraw 函数用于对JackStraw 分析的结果进行打分。
它会计算每个主成分的p值，用于评估该主成分的显著性 # dims指定要评估的主成分范围，这里表示评估第 1 到第 20 个主成分 pbmc <- ScoreJackStraw (pbmc, dims = 1 : 20 ) # 该函数用于绘制JackStraw分析的结果图。图形会展示每个主成分的p值分布，帮助我们直观地判断哪些主成分是显著的，哪些是随机噪声。x轴经验值，y轴理论值黑色斜线表示零假设下的均匀p值分布。
其作用是作为判断主成分（PC）显著性的基准线，具体含义如下：零假设的参考JackStraw分析通过随机置换数据生成“零分布”（即假设主成分中基因的表达模式是随机的），黑色虚线即代表这种均匀分布的p值预期。如果主成分的实际p值分布（实线）显著偏离虚线（尤其是位于虚线之上），说明该主成分包含的基因表达模式并非随机，具有统计学意义。实线高于虚线：表示该主成分中低p值的基因富集，说明其包含真实的生物学信号，应保留用于后续分析（如聚类、可视化）。
实线接近或低于虚线：说明主成分的 p 值分布与随机情况无显著差异，可能主要反映噪声，可忽略。JackStrawPlot (pbmc, dims = 1 : 20 ) # ElbowPlot 函数用于绘制肘部图。肘部图展示了每个主成分的方差解释率随主成分数量的变化情况。
通常，我们会在肘部图中找到一个“肘部”位置，即方差解释率下降趋势变缓的点，该点对应的主成分数量可以作为后续分析中保留的主成分数量 # 例如这里可以选择前10个主成分 ElbowPlot (pbmc) # FindNeighbors 函数用于在低维空间（这里是 PCA 空间）中计算细胞之间的邻接关系。
它会根据指定的主成分（这里是第 1 到第 10 个主成分）计算细胞之间的距离，并构建一个 K 近邻图（K-nearest neighbor graph），用于后续的聚类分析。pbmc <- FindNeighbors (pbmc, dims = 1 : 10 ) # FindClusters 函数用于对细胞进行聚类分析。
它基于 FindNeighbors 函数构建的 K 近邻图，使用图聚类算法（如 Louvain 算法）将相似的细胞聚为不同的簇。# resolution = 0.5：指定聚类的分辨率。分辨率越高，聚类的粒度越细，会得到更多的小簇；分辨率越低，聚类的粒度越粗，会得到较少的大簇。
pbmc <- FindClusters (pbmc, resolution = 0.5 ) # Modularity Optimizer version 1.3.0 by Ludo Waltmand Nees Jan van Eck # # Number of nodes: 2638 # Number of edges: 95965 # # Running Louvain algorithm... # Maximum modularity in 10 random starts: 0.8723 # Number of communities: 9 # Elapsed time: 0 seconds 如果你在上面的FindNeighbors中发生了 object 'CsparseMatrix_validate' not found 报错，可以参考下面的方法解决：答读者问| 22.object ‘CsparseMatrix_validate’ not found。
如果你用的是 Seurat V5 那么你大概率不会遇到上面的问题。2.5 分群后的可视化 2.5.1 tsne+umap # 原理：# UMAP：基于流形学习的思想，假设数据分布在一个低维流形上。它通过构建数据点之间的局部连接图来近似流形结构，然后优化低维空间中的点的位置，使得低维空间中的点之间的连接关系尽可能与高维空间中的一致。其核心是使用模糊拓扑来保留数据的全局和局部结构。# t-SNE：通过概率分布来表示数据点之间的相似性。
在高维空间中，使用高斯分布来计算数据点之间的相似度；在低维空间中，使用t分布来计算相似度。t-SNE的目标是最小化高维空间和低维空间中数据点对之间的概率分布差异，通常使用KL散度来衡量这种差异。# 计算效率：# UMAP：通常比t-SNE 更快，尤其是在处理大规模数据集时。UMAP的优化过程相对简单，并且可以使用近似算法来加速计算，能够在较短的时间内完成大规模数据的降维任务。
# t-SNE：计算复杂度较高，尤其是对于大规模数据集，计算时间会显著增加。它的优化过程是迭代的，并且需要计算所有数据点对之间的相似度，因此在处理大数据集时效率较低。
pbmc <- RunUMAP (pbmc, dims = 1 : 10 ) DimPlot (pbmc, reduction = "umap" ) pbmc <- RunTSNE (pbmc, dims = 1 : 10 ) DimPlot (pbmc, reduction = "tsne" ) 2.5.2 寻找marker基因并对cluster进行重命名 # 计算cluster5与cluster(0+3)群之间的差异基因 cluster5.markers <- FindMarkers (pbmc, ident.1 = 5 , ident.2 = c ( 0 , 3 ), min.pct = 0.25 ) head (cluster5.markers) # p_val avg_log2FC pct.1 pct.2 p_val_adj # FCGR3A 2.150929e-209 6.832372 0.975 0.039 2.949784e-205 # IFITM3 6.103366e-199 6.181000 0.975 0.048 8.370156e-195 # CFD 8.891428e-198 6.052575 0.938 0.037 1.219370e-193 # CD68 2.374425e-194 5.493138 0.926 0.035 3.256286e-190 # RP11-290F20.3 9.308287e-191 6.335402 0.840 0.016 1.276538e-186 # IFI30 1.575100e-181 6.382751 0.796 0.014 2.160093e-177 # 依次计算每一个亚群与剩余所有亚群的差异基因 pbmc.markers <- FindAllMarkers (pbmc, only.pos = TRUE , min.pct = 0.25 , logfc.threshold = 0.25 ) if ( ! require (dplyr)) install.packages ( "dplyr" ) # 每个亚群按照avg_log2FC排序得到的top2差异基因 # wt 参数代表权重变量（weight variable），它指定了用于排序的变量 pbmc.markers %>% group_by (cluster) %>% top_n ( n = 2 , wt = avg_log2FC) # # A tibble: 18 × 7 # # Groups: cluster [9] # p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene # <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <chr> # 1 1.17e- 83 2.37 0.435 0.108 1.60e- 79 0 CCR7 # 2 3.28e- 49 2.10 0.333 0.103 4.50e- 45 0 LEF1 # 3 3.10e-139 7.28 0.3 0.004 4.25e-135 1 FOLR3 # 4 1.63e-121 6.74 0.277 0.006 2.23e-117 1 S100A12 # 5 2.61e- 59 2.11 0.424 0.111 3.58e- 55 2 AQP3 # 6 1.94e- 35 1.90 0.267 0.069 2.66e- 31 2 CD40LG # 7 2.40e-272 7.38 0.564 0.009 3.29e-268 3 LINC00926 # 8 2.75e-237 7.14 0.488 0.007 3.76e-233 3 VPREB3 # 9 4.93e-169 4.37 0.595 0.056 6.76e-165 4 GZMK # 10 3.51e- 93 3.59 0.437 0.06 4.82e- 89 4 GZMH # 11 1.69e-212 5.43 0.506 0.01 2.32e-208 5 CDKN1C # 12 8.23e-168 5.88 0.37 0.005 1.13e-163 5 CKB # 13 1.05e-265 6.01 0.986 0.071 1.44e-261 6 GZMB # 14 7.57e-183 6.21 0.493 0.014 1.04e-178 6 AKR1C3 # 15 1.48e-220 7.63 0.812 0.011 2.03e-216 7 FCER1A # 16 1.46e-207 8.03 0.5 0.002 2.00e-203 7 SERPINF1 # 17 0 14.3 0.571 0 0 8 LY6G6F # 18 4.36e-206 13.8 0.357 0 5.98e-202 8 RP11-879F14.2 top10 <- pbmc.markers %>% group_by (cluster) %>% top_n ( n = 10 , wt = avg_log2FC) DoHeatmap (pbmc, features = top10 $ gene) + NoLegend # 细胞类型marker基因FeaturePlot p1 <- FeaturePlot (pbmc, features = c ( "CCR7" , "CD14" , "LYZ" , Mono "MS4A1" , "CD8A" , "FCGR3A" , Mono "NKG7" , "FCER1A" , "PF4" )) p2 <- DimPlot (pbmc) library (gridExtra) p1 | p2 可以在这个 CellMarker 中查看细胞类型对应的marker # 定义细胞名称 new.cluster.ids <- c ( "Naive CD4 T" , "CD14+ Mono" , "Memory CD4 T" , "B" , "CD8 T" , "FCGR3A+ Mono" , "NK" , "DC" , "Platelet" ) new.cluster.ids # [1] "Naive CD4 T" "CD14+ Mono" "Memory CD4 T" "B" "CD8 T" # [6] "FCGR3A+ Mono" "NK" "DC" "Platelet" levels (pbmc) # [1] "0" "1" "2" "3" "4" "5" "6" "7" "8" # 将细胞名和分群信息一一对应 names (new.cluster.ids) <- levels (pbmc) new.cluster.ids # 0 1 2 3 4 # "Naive CD4 T" "CD14+ Mono" "Memory CD4 T" "B" "CD8 T" # 5 6 7 8 # "FCGR3A+ Mono" "NK" "DC" "Platelet" # 将new.cluster.ids作为新的细胞名称 pbmc <- RenameIdents (pbmc, new.cluster.ids) pbmc $ celltype <- Idents (pbmc) DimPlot (pbmc, reduction = "umap" , label = TRUE , pt.size = 0.5 ) + NoLegend # 保存为rds文件，便于后续使用 # saveRDS(pbmc,'../result_seuratV5/pbmc.rds') pbmc <- readRDS ( '../result_seuratV5/pbmc.rds' ) DimPlot (pbmc, reduction = "umap" , label = TRUE , pt.size = 0.5 ) + NoLegend