单细胞RNA测序(scRNA-seq)分析看似复杂，其实有一套标准流程。 本教程按照真实项目顺序， 从数据读取、质控、标准化到聚类、注释、可视化，一步步带你完成首个单细胞分析！
📋 分析流程概览
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(pbmc@meta.data)
⚠️ 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
- • 可以使用
PCTechPlot辅助判断 - • 不同数据集可能需要不同数量的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)