📖 单细胞RNA测序(scRNA-seq)分析看似复杂，其实有一套标准流程。本教程按照真实项目顺序，从数据读取、质控、标准化到聚类、注释、可视化，一步步带你完成首个单细胞分析！
📋 分析流程概览 1.数据读取 → 2.质控 → 3.标准化 → 4.降维 → 5.聚类 → 6.注释 Step 0: 环境准备 # 安装Seurat包 (最新版v5) install.packages ( "Seurat" ) # 加载必要的包 library (Seurat) library (dplyr) library (patchwork) # 用于拼图 # 检查版本 packageVersion ( "Seurat" ) # 建议使用 v4 或 v5 💡 内存要求：单细胞数据通常较大，建议至少16GB内存。
如果是大型数据集(>100k细胞)，建议32GB以上。
Step 1: 数据读取 # 10x Genomics 数据 (最常用) # 数据目录应包含: barcodes.tsv, genes.tsv, matrix.mtx pbmc.data <- Read10X (data.dir = "path/to/filtered_feature_bc_matrix/" ) # 查看数据维度 dim (pbmc.data) # [1] 13714 2700 # 13714个基因, 2700个细胞 # 创建Seurat对象 pbmc <- CreateSeuratObject ( counts = pbmc.data, project = "PBMC" , # 项目名 min.cells = 3 , # 基因至少在3个细胞表达 min.features = 200 # 细胞至少表达200个基因 ) # 查看对象摘要 pbmc 📌 数据过滤参数：• min.cells：去除低表达基因（可能是噪音） • min.features：去除低质量细胞（可能是空液滴） • 这两个参数根据具体数据调整 Step 2: 质量控制 (QC) # 计算线粒体基因比例 (关键QC指标) pbmc[[ "percent.mt" ]] <- PercentageFeatureSet (pbmc, pattern = "^MT-" ) # 计算核糖体基因比例 (可选) pbmc[[ "percent.rb" ]] <- PercentageFeatureSet (pbmc, pattern = "^RP[SL]" ) # 可视化QC指标 VlnPlot (pbmc, features = c ( "nFeature_RNA" , "nCount_RNA" , "percent.mt" ), ncol = 3 ) # 查看QC指标分布 head ([email protected]) ⚠️ QC指标含义：• nFeature_RNA：检测到的基因数，过低可能是空液滴，过高可能是双细胞 • nCount_RNA：UMI总数，反映测序深度 • percent.mt：线粒体基因比例，高比例表示细胞可能已死亡或受损 # 过滤低质量细胞 (根据实际数据调整阈值) pbmc <- subset (pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5 " ) # 查看过滤后剩余细胞数 ncol (pbmc) Step 3: 数据标准化 # 方法1: 标准LogNormalize (经典方法) pbmc <- NormalizeData ( pbmc, normalization.method = "LogNormalize" , scale.factor = 10000 # CPM scaling ) # 方法2: SCTransform (推荐，v5默认) # pbmc <- SCTransform(pbmc, vars.to.regress = "percent.mt") # 识别高变基因 (HVGs) pbmc <- FindVariableFeatures ( pbmc, selection.method = "vst" , # 推荐方法 nfeatures = 2000 ) # 可视化高变基因 VariableFeaturePlot (pbmc) + theme (legend.position = "none" ) 💡 高变基因的重要性：高变基因是细胞间差异表达的主要来源，它们定义了细胞的异质性。
后续降维和聚类只使用这些基因，可以降低噪音并提高计算效率。
Step 4: 降维分析 # 数据缩放 (Scale) all.genes <- rownames (pbmc) pbmc <- ScaleData (pbmc, features = all.genes) # 线性降维: PCA pbmc <- RunPCA ( pbmc, features = VariableFeatures (object = pbmc), npcs = 50 # 计算前50个主成分 ) # 可视化PCA (Elbow Plot) ElbowPlot (pbmc, ndims = 50 ) 📌 如何选择PC数量？
：• 查看 Elbow Plot，在"拐点"之前选择 • 通常选择前10-30个PC • 可以使用 ElbowPlot 辅助判断 • 不同数据集可能需要不同数量的PC # 非线性降维: UMAP 和 t-SNE pbmc <- RunUMAP (pbmc, dims = 1 : 20 ) pbmc <- RunTSNE (pbmc, dims = 1 : 20 ) # 可视化降维结果 DimPlot (pbmc, reduction = "umap" ) Step 5: 细胞聚类 # 构建KNN图 pbmc <- FindNeighbors ( pbmc, dims = 1 : 20 # 使用前20个PC ) # 聚类 (Louvain算法) pbmc <- FindClusters ( pbmc, resolution = 0.5 # 聚类分辨率 ) # 查看聚类结果 table (pbmc$seurat_clusters) # UMAP上显示聚类 DimPlot (pbmc, reduction = "umap" , label = TRUE ) 💡 resolution 参数：控制聚类粒度：• 低值(0.1-0.5)：少量大cluster • 高值(0.8-2.0)：更多小cluster • 可以尝试多个值，选择生物学意义最合理的 Step 6: 寻找Marker基因 # 寻找每个cluster的marker基因 pbmc.markers <- FindAllMarkers ( pbmc, only.pos = TRUE , # 只找上调基因 min.pct = 0.25 , # 最少25%细胞表达 logfc.threshold = 0.25 # logFC阈值 ) # 查看top markers top.markers <- pbmc.markers %>% group_by (cluster) %>% top_n (n = 5 , wt = avg_log2FC) # 可视化top markers DotPlot (pbmc, features = c ( "CD3D" , "CD8A" , "CD4" , "MS4A1" , "LYZ" )) + RotatedAxis () 📌 Marker基因的作用：Marker基因是定义细胞类型的"身份证"。
通过查找文献或数据库(如CellMarker、PanglaoDB)，可以确定每个cluster对应的细胞类型。
Step 7: 细胞类型注释 # 根据marker基因手动注释 # 需要结合生物学知识和文献 # 定义新的细胞类型标签 new.cluster.ids <- c ( "Naive CD4 T" , "Memory CD4 T" , "CD14+ Mono" , "B" , "CD8 T" , "FCGR3A+ Mono" , "NK" , "DC" , "Platelet" ) # 重命名cluster names (new.cluster.ids) <- levels (pbmc) pbmc <- RenameIdents (pbmc, new.cluster.ids) DimPlot (pbmc, reduction = "umap" , label = TRUE , pt.size = 0.5 ) + NoLegend () 💡 自动注释工具：如果手动注释困难，可以使用：• SingleR：基于参考数据集的自动注释 • scCATCH：自动匹配marker到细胞类型 • CellTypist：机器学习注释工具 Step 8: 高级可视化 # 1. FeaturePlot: 基因表达图 FeaturePlot (pbmc, features = c ( "MS4A1" , "CD3D" )) # 2. VlnPlot: 小提琴图 VlnPlot (pbmc, features = c ( "MS4A1" , "CD3D" ), pt.size = 0.1 ) # 3. RidgePlot: 山脊图 RidgePlot (pbmc, features = c ( "MS4A1" , "CD3D" )) # 4. Heatmap: marker热图 top10 <- pbmc.markers %>% group_by (cluster) %>% top_n (n = 10 , wt = avg_log2FC) DoHeatmap (pbmc, features = top10$gene) + NoLegend () 📌 发表级图表建议：• Figure 1：QC violin plot + UMAP聚类图 • Figure 2：Marker基因热图 + Dot plot • Figure 3：FeaturePlot展示关键基因 • 使用 ggsave() 保存高分辨率图片 📝 完整工作流代码总结 # 一键运行完整流程 library (Seurat) # 1. 读取数据 data <- Read10X ( "data/dir" ) pbmc <- CreateSeuratObject (data) # 2. 质控 pbmc[[ "percent.mt" ]] <- PercentageFeatureSet (pbmc, "^MT-" ) pbmc <- subset (pbmc, subset = percent.mt < 5 ) # 3. 标准化 + HVG pbmc <- SCTransform (pbmc) # 4. 降维 pbmc <- RunPCA (pbmc) pbmc <- RunUMAP (pbmc, dims = 1 : 20 ) # 5. 聚类 pbmc <- FindNeighbors (pbmc, dims = 1 : 20 ) pbmc <- FindClusters (pbmc) # 6. 可视化 DimPlot (pbmc, reduction = "umap" , label = TRUE )