1. 从零开始理解单细胞转录组分析的核心脉络如果你刚接触单细胞转录组分析面对海量的细胞和基因数据可能会觉得无从下手。别担心我刚开始做的时候也是一头雾水感觉像是在大海里捞针。但经过这些年的实战我总结出一条清晰的路径数据预处理 - 降维 - 聚类 - 注释。这四步就像盖房子的地基、框架、隔墙和装修每一步都至关重要一步走错后面的结果可能就全乱了。简单来说单细胞转录组分析就是给成千上万个细胞“测个序”看看每个细胞里哪些基因在活跃地“说话”表达。我们的目标是从这堆嘈杂的数据里把功能、状态相似的细胞找出来归成一类并搞清楚它们到底是什么类型的细胞。比如从一堆免疫细胞里把T细胞、B细胞、巨噬细胞区分开并研究它们各自在干什么。这个过程为什么需要Seurat和Harmony呢你可以把Seurat想象成一个功能极其强大的“细胞数据分析工厂”。从原始数据读入、质量控制、标准化到降维、聚类、可视化它提供了一条龙服务。而Harmony则是这个工厂里一个非常聪明的“批次效应校正工程师”。当我们把来自不同实验、不同时间、不同测序平台的数据也就是不同批次放在一起分析时细胞之间的差异可能不是生物学差异而是这些技术因素造成的“噪音”。Harmony的作用就是巧妙地消除这些噪音让不同批次的细胞能够基于真实的生物学信号进行公平比较和整合。这篇文章我就带你手把手走一遍这个实战流程。我会用我踩过的坑、试过的参数、验证过的经验让你不仅能跑通代码更能理解每一步背后的“为什么”。准备好了吗我们这就开始。2. 实战第一步数据读取与Seurat对象构建拿到数据第一步不是急着跑分析而是先“认识”你的数据。通常我们从公共数据库比如GEO下载的数据是压缩的矩阵文件每个样本一个文件。我的习惯是先在一个干净的目录里解压然后系统性地读取。2.1 环境准备与数据读取首先确保你的R环境里安装了必要的包。我建议直接用install.packages和BiocManager::install来装别在版本问题上纠结太久。# 清空环境避免旧变量干扰 rm(listls()) # 加载核心包 library(Seurat) # 单细胞分析主力 library(harmony) # 批次整合神器 library(patchwork) # 图片排版 library(tidyverse) # 数据操作全家桶用着顺手 library(data.table) # 大文件读取快假设你的数据文件叫GSE167297_RAW.tar里面是一堆.txt.gz文件。解压后用list.files获取所有文件名然后用lapply循环读取。这里有个关键点创建Seurat对象时给每个样本起个好记的project名。我习惯把文件名里冗余的“GSM数字”和“_CountMatrix.txt.gz”后缀去掉只保留有生物学意义的标识。# 解压数据 untar(GSE167297_RAW.tar, exdir GSE167297_RAW) # 设置数据目录 data_dir - /path/to/your/GSE167297_RAW/ sample_files - list.files(data_dir) print(sample_files) # 先看看有哪些文件心里有数 # 循环读取每个文件并创建Seurat对象 sce_list - lapply(sample_files, function(file_name) { # 读入数据data.table的fread速度很快 count_matrix - fread(file.path(data_dir, file_name), data.table FALSE) # 设置行名基因名通常在第一列 rownames(count_matrix) - count_matrix[, 1] count_matrix - count_matrix[, -1] # 去掉基因名列剩下纯表达矩阵 # 创建Seurat对象 # min.cells和min.features是初始质控过滤掉在太少细胞中表达的基因和基因数太少的细胞 seurat_obj - CreateSeuratObject(counts count_matrix, project gsub(_CountMatrix.txt.gz, , gsub(^GSM[0-9]*_, , file_name)), min.cells 5, min.features 300) # 注意原代码里的 subset(sce, downsample700) 是用于演示时减少计算量自己分析时务必去掉 return(seurat_obj) })2.2 合并数据与初步检查单个样本对象创建好后我们需要把它们合并成一个大的Seurat对象。merge函数可以轻松完成这个任务add.cell.ids参数可以为每个来源的细胞添加前缀防止细胞名重复。# 合并所有样本 sce_combined - merge(x sce_list[[1]], y sce_list[-1], add.cell.ids gsub(_CountMatrix.txt.gz, , gsub(^GSM[0-9]*_, , sample_files))) # 查看合并后的数据规模 dim(sce_combined) # 输出类似[1] 18207 9301表示有18207个基因9301个细胞到这里一个包含所有数据的Seurat对象就构建好了。它就像一个大仓库里面整齐地存放着每个细胞的基因表达信息以及它们来自哪个样本orig.ident的元数据。接下来我们就要对这个仓库里的“货物”进行质量检查把破损的、不合格的剔除掉。3. 至关重要的质量控制识别并过滤低质量细胞质量控制QC是分析中最容易被忽视但又最致命的一环。用不好的数据再高级的算法也出不来好结果。QC主要看几个指标每个细胞检测到的基因数nFeature_RNA、总分子数nCount_RNA、以及线粒体基因percent.mt、核糖体基因percent.rp、血红蛋白基因percent.hb的占比。3.1 计算QC指标并可视化线粒体基因高通常意味着细胞状态不好比如正在凋亡或已经死亡核糖体基因高可能和细胞活跃的翻译状态有关但也需结合看血红蛋白基因如果在非红细胞中高表达可能是红细胞污染。我们先计算这些指标sce - sce_combined # 赋值给新变量方便后续操作 # 计算线粒体基因比例人类基因前缀为MT-小鼠为mt- sce[[percent.mt]] - PercentageFeatureSet(sce, pattern ^MT-) # 计算核糖体基因比例 sce[[percent.rp]] - PercentageFeatureSet(sce, pattern ^RP[SL]) # 计算血红蛋白基因比例 sce[[percent.hb]] - PercentageFeatureSet(sce, pattern ^HB[^(P)])算出来不能光看数字要画图。小提琴图VlnPlot可以直观地看每个样本这些指标的分布。# 绘制QC指标的小提琴图按样本分组 qc_violin - VlnPlot(sce, features c(nFeature_RNA, nCount_RNA, percent.mt, percent.rp, percent.hb), ncol 3, pt.size 0.1, # 点调小些图更清晰 group.by orig.ident) ggsave(QC_VlnPlot.png, plot qc_violin, width 10, height 8, dpi 300)光看分布还不够要看指标间的相关性。比如nCount_RNA和nFeature_RNA应该呈较强的正相关测到的分子越多检测到的基因也越多。而nCount_RNA和percent.mt最好是弱相关或负相关如果强正相关说明大量RNA来自线粒体细胞质量堪忧。# 绘制指标间的散点图 plot1 - FeatureScatter(sce, feature1 nCount_RNA, feature2 percent.mt, group.by orig.ident) plot2 - FeatureScatter(sce, feature1 nCount_RNA, feature2 nFeature_RNA, group.by orig.ident) plot3 - FeatureScatter(sce, feature1 nCount_RNA, feature2 percent.rp, group.by orig.ident) plot4 - FeatureScatter(sce, feature1 percent.rp, feature2 percent.mt, group.by orig.ident) combined_scatter - plot1 plot2 plot3 plot4 ggsave(QC_FeatureScatter.png, plot combined_scatter, width 16, height 10, dpi 300)3.2 设定阈值并过滤细胞看图说话设定过滤阈值。这里没有绝对的金标准需要根据你的数据和生物学背景灵活调整。比如如果你的细胞代谢活跃如心肌细胞线粒体基因比例天生就高阈值就要放宽。我通常的起始参考范围是nFeature_RNA过滤掉基因数过少如200和过多可能是双细胞或多细胞如6000的细胞。nCount_RNA过滤掉总分子数异常高或低的细胞。percent.mt一般设定上限比如10%-20%。我的示例数据用了18%。percent.hb如果不是研究红细胞可以设一个很低的上限如1%过滤掉污染。# 根据散点图观察结果应用过滤条件 sce_filtered - subset(sce, subset nFeature_RNA 6000 nCount_RNA 30000 percent.mt 18 percent.rp 15 # 这个阈值根据情况调整 percent.hb 1) # 过滤后记得检查细胞数 dim(sce_filtered)踩坑提醒过滤不要太狠一次过滤掉超过50%的细胞就要回头检查数据或阈值是否合理。稳妥的做法是先宽后严分步过滤每次过滤后都重新画图看看分布变化。4. 数据标准化、高变基因与Harmony整合经过QC我们得到了相对干净的数据。但不同细胞之间的测序深度总分子数差异很大直接比较会引入偏差。所以我们需要标准化。4.1 数据标准化与高变基因筛选Seurat的NormalizeData函数默认使用“LogNormalize”方法将每个细胞的表达量除以该细胞的总表达量乘以一个缩放因子默认为10000然后进行对数转换。这能消除测序深度的影响。sce - sce_filtered sce - NormalizeData(sce, normalization.method LogNormalize, scale.factor 10000)接下来是找高变基因HVGs。不是所有基因都包含有用的信息很多基因在所有细胞里表达量都差不多管家基因或者表达量极低噪音。高变基因是那些在不同细胞间表达差异大的基因它们往往驱动着细胞间的生物学差异。我们只对这些基因进行后续分析可以大幅降维并聚焦信号。# 鉴定2000个高变基因vst方法是目前最常用的 sce - FindVariableFeatures(sce, selection.method vst, nfeatures 2000) # 看看高变基因都是谁 top10_hvgs - head(VariableFeatures(sce), 10) print(top10_hvgs) # 可视化高变基因 var_plot - VariableFeaturePlot(sce) labeled_plot - LabelPoints(plot var_plot, points top10_hvgs, repel TRUE) ggsave(Variable_Features.png, plot labeled_plot, width 10, height 6, dpi 300)找到高变基因后需要进行缩放Scaling。这一步的目的是让每个基因在所有细胞中的表达量均值为0方差为1。这样在后续计算细胞距离时高表达基因就不会因为绝对值大而占据主导地位。特别注意默认只缩放高变基因以节省计算资源但如果你计划后续做差异表达分析涉及所有基因最好缩放所有基因。# 获取所有基因名 all_genes - rownames(sce) # 缩放数据这里对所有基因进行缩放为后续分析做准备 sce - ScaleData(sce, features all_genes)4.2 PCA降维与Harmony批次校正缩放后的数据维度依然很高几千个高变基因。我们需要用主成分分析PCA进一步降维提取最主要的变异方向。# 基于高变基因进行PCA降维 sce - RunPCA(sce, features VariableFeatures(object sce), verbose FALSE) # 可以可视化前几个主成分 # VizDimLoadings(sce, dims 1:2, reduction pca) # DimPlot(sce, reduction pca, group.by orig.ident)如果数据来自多个批次PCA图可能会显示样本间明显的分离这很可能就是批次效应。这时就该Harmony出场了。它的原理简单理解就是在保留生物学差异的前提下把不同批次“对齐”到同一个低维空间。# 运行Harmony指定批次变量通常是样本来源orig.ident sce - RunHarmony(sce, group.by.vars orig.ident, reduction pca) # 可视化Harmony校正后的结果 harmony_plot - DimPlot(sce, reduction harmony, group.by orig.ident) ggsave(Harmony_Integration.png, plot harmony_plot, width 8, height 6, dpi 300)运行后你会看到sce对象里多了一个harmony降维结果。比较一下PCA和Harmony的图理想情况下Harmony图中不同颜色的点代表不同批次应该混合得更好而基于生物学类型的聚类结构应该更清晰。4.3 如何确定使用多少个主成分PCs降维后我们需要决定用前多少个PC来代表数据。用太少会丢失信息用太多会引入噪音。对于Harmony校正后的数据我们使用ElbowPlot肘部图来辅助判断。# 绘制肘部图查看每个PC解释的方差百分比 elbow_plot - ElbowPlot(sce, ndims 50, reduction harmony) # 通常看前50个PC就够了 ggsave(ElbowPlot_Harmony.png, plot elbow_plot, width 8, height 6, dpi 300)肘部图曲线拐点像手肘一样的位置通常是一个合理的PC数量选择。但这不是绝对的有时拐点不明显。我的经验是对于单细胞数据选择范围通常在10到50之间。你可以选择一个拐点附近的数值比如图中曲线在PC30之后趋于平缓那么dims 30可能是个不错的选择。也可以多试几个值看看对后续聚类稳定性的影响。5. 细胞聚类与可视化找到自然的细胞群体降维之后细胞在高维空间中的距离更近了。接下来我们要根据它们之间的相似性距离进行聚类把相似的细胞归为同一个群体。5.1 构建邻近图与聚类首先基于选定的Harmony维度计算细胞间的邻近关系图。# 假设我们根据肘部图选择前30个Harmony维度 selected_dims - 30 sce - FindNeighbors(sce, dims 1:selected_dims, reduction harmony, verbose FALSE)然后使用FindClusters函数进行聚类。这里有个关键参数resolution分辨率。分辨率值越大得到的聚类数越多划分越细值越小聚类数越少划分越粗。没有唯一正确的值需要根据你的生物学问题来探索。# 尝试一系列分辨率观察聚类数如何变化 sce - FindClusters(sce, resolution seq(0.1, 2, by 0.1), verbose FALSE) # 查看默认分辨率通常是最后一个下的聚类结果 head(Idents(sce), 5)5.2 使用Clustree选择最佳分辨率面对一串分辨率怎么选clustree包可以可视化不同分辨率下聚类结果的演变关系帮你做出更明智的选择。library(clustree) # 绘制聚类树 tree_plot - clustree(sce, prefix RNA_snn_res.) coord_flip() # 选择一个分辨率例如2查看其UMAP图 sce - FindClusters(sce, resolution 2, verbose FALSE) # 重新以选定分辨率运行 cluster_plot - DimPlot(sce, reduction harmony, label TRUE) combined_tree_plot - tree_plot cluster_plot plot_layout(widths c(3, 1)) ggsave(Clustree_Analysis.png, plot combined_tree_plot, width 20, height 12, dpi 300)在聚类树上寻找那些“树干”稳定、分支清晰的分辨率。避免选择细胞簇频繁分裂或合并的分辨率。选好后就用该分辨率重新运行FindClusters。5.3 UMAP/t-SNE可视化聚类结果需要在二维图上展示最常用的就是UMAP和t-SNE。UMAP通常能更好地保持全局结构我更喜欢用。# 基于相同的维度进行UMAP降维 sce - RunUMAP(sce, dims 1:selected_dims, reduction harmony, umap.method uwot, metric cosine) # 可视化聚类结果 umap_plot - DimPlot(sce, reduction umap, label TRUE, pt.size 0.5) ggsave(UMAP_Clusters.png, plot umap_plot, width 10, height 8, dpi 300)现在你得到了第一张单细胞图谱不同颜色的点代表不同的细胞簇。但这只是数字编号0,1,2...我们还需要知道每个簇具体是什么细胞类型。6. 细胞注释为细胞簇赋予生物学身份这是整个分析中最需要生物学知识的一步也是最有意思的一步。我们要把冷冰冰的聚类编号变成有生命的细胞类型。6.1 寻找每个簇的标记基因首先用FindAllMarkers函数找出每个簇相对于其他所有簇特异性高表达的基因标记基因。# 寻找所有簇的阳性标记基因 sce.markers - FindAllMarkers(sce, only.pos TRUE, # 只找高表达的 min.pct 0.25, # 基因至少在25%的该簇细胞中表达 logfc.threshold 0.25, # 对数倍变化阈值 verbose FALSE) # 查看结果 head(sce.markers) # 保存结果 write.csv(sce.markers, all_markers.csv, row.names TRUE) # 提取每个簇表达量最高的前5个基因 top5_markers - sce.markers %% group_by(cluster) %% top_n(n 5, wt avg_log2FC) write.csv(top5_markers, top5_markers_per_cluster.csv, row.names TRUE)6.2 人工注释基于已知标记基因这是最经典的方法。你需要根据已有的生物学知识查看每个簇高表达的标记基因来判断其细胞类型。比如看到CD3D、CD3E高表达很可能是T细胞CD19、MS4A1CD20高表达很可能是B细胞。我们可以用点图DotPlot来一次性可视化多个已知标记基因在所有簇中的表达模式。# 定义一组感兴趣的已知标记基因 known_markers - c(PTPRC, CD3D, CD3E, CD4, CD8A, # T细胞 CD19, CD79A, MS4A1, # B细胞 CD68, CD163, CD14, # 髓系细胞单核/巨噬 NKG7, GNLY, # NK细胞 EPCAM, KRT19, # 上皮细胞 PECAM1, VWF, # 内皮细胞 ACTA2, MYH11, # 成纤维细胞/平滑肌 MKI67, TOP2A) # 增殖细胞 dot_plot - DotPlot(sce, features known_markers, group.by seurat_clusters) theme(axis.text.x element_text(angle 45, hjust 1)) coord_flip() ggsave(Marker_DotPlot.png, plot dot_plot, width 12, height 10, dpi 300)根据点图结合top5_markers列表你就可以开始给每个簇贴标签了。例如一个簇高表达CD3D,CD8A低表达CD4那它很可能是CD8 T细胞。6.3 自动注释使用SingleR等工具对于新手或者面对不熟悉的组织自动注释工具能提供很好的参考。SingleR是一个常用的工具它通过将你的单细胞数据与已注释好的参考数据集如Blueprint、Human Primary Cell Atlas进行比较来预测细胞类型。library(SingleR) library(celldex) # 加载参考数据集以Blueprint为例 ref - BlueprintEncodeData() # 提取待注释数据的归一化表达矩阵 test_data - GetAssayData(sce, layer data) # Seurat v5使用 layerdata # 运行SingleR进行注释 pred - SingleR(test test_data, ref ref, labels ref$label.main, # 使用参考数据的主标签 clusters Idents(sce)) # 按我们已有的聚类结果进行注释 # 查看注释结果 print(pred$pruned.labels) # 将注释结果添加到Seurat对象 new_ids - pred$pruned.labels names(new_ids) - levels(sce) sce - RenameIdents(sce, new_ids) # 可视化自动注释结果 auto_anno_plot - DimPlot(sce, reduction umap, label TRUE, pt.size 0.5) ggsave(UMAP_AutoAnnotation.png, plot auto_anno_plot, width 10, height 8, dpi 300)重要提示自动注释的结果绝对不能盲从一定要结合人工检查标记基因的表达。参考数据库可能不包含你研究体系中的所有细胞类型或者注释粒度不符合你的需求。自动注释结果应该作为你人工注释的“第二意见”或起点。6.4 合并与精细化注释有时自动注释或初步人工注释后你会发现有些簇的标记基因很相似比如都是T细胞可能需要合并。或者一个簇内可能包含不止一种细胞类型需要提高分辨率重新分群。# 示例手动合并注释结果假设经过检查簇0和簇1都是CD4 T细胞 new_cluster_ids - c(CD4_T, CD4_T, CD8_T, B, Monocyte, NK, Fibroblast) names(new_cluster_ids) - levels(sce) sce - RenameIdents(sce, new_cluster_ids) # 或者将注释信息添加到元数据中方便后续分组分析 sce$celltype - Idents(sce) # 保存最终注释好的对象 saveRDS(sce, file annotated_seurat_object.rds)至此一个完整的单细胞转录组分析流程就走完了。从原始数据到一张标注好细胞类型的UMAP图谱你已经掌握了最核心的实战技能。记住参数不是固定的生物学背景才是根本。多看图多思考多结合文献你的注释才会越来越准。