写一篇关于单细胞数据分析第四个核心环节的实战教程是时候把“降维、聚类、注释”这三个最深也最常让人卡壳的步骤讲透了。特别是刚跑完前面三步预处理质控、标准化、高变基因识别正准备往深水区走的朋友这篇文字就是你的过渡桥。老规矩这不是复读文档而是贴近真实分析现场的经验记录。从PCA到底选多少个主成分到UMAP和tSNE怎么选到Leiden聚类分辨率调到多少才算“对”再到注释细胞类型时怎么把marker基因表用活——我会按实操顺序一点点捋。1. 分析起点预处理后的adata到底长什么样如果你是从系列第一篇跟过来的此刻手里的adata对象大概已经完成了这几步原始counts矩阵经过了pp.filter_cells和pp.filter_genes的粗过滤用pp.normalize_total做了文库大小校正默认target_sum1e4又通过pp.log1p做了对数化然后再用pp.highly_variable_genes筛出了一批高变基因。这个阶段的adata结构大致是这样adata.shape # (细胞数, 基因数) 例如 (8291, 18924) 这类规模 adata.var[highly_variable].value_counts() # True ~2000-5000 # False ~15000注意高变基因数量默认是2000-5000。为什么所有教程都强调这一步因为后边的PCA、聚类这些计算全部只在高变基因子集上进行。这个操作有两层含义第一计算上更干净。高变基因是那些在细胞间表达差异显著的基因它们才携带了区分细胞类型的主要信号。那些在所有细胞里都恒定表达或几乎不表达的基因管家基因、零表达基因放进去只会稀释信号增加噪声。第二速度上更快。单细胞数据动辄几千上万个细胞基因数近两万全量跑PCA和邻居图计算代价很大。缩减到两三千个高变基因计算量降了一个量级而且结果往往更稳。不过有件事得提醒是否所有数据都必须用高变基因不一定。如果你的研究就盯着某个特定基因集比如免疫相关的几百个基因完全可以绕过全局高变基因筛选直接基于自定义基因集做后续分析。但这种做法的前提是你对研究假设非常笃定。对通用探索性分析来说老老实实跑一遍高变基因筛选是最不容易出错的路线。另外一个小细节很多人会忽略在筛选高变基因后非常建议把adata.raw存一份原始对数化矩阵。因为之后你做的任何降维、聚类、标记基因识别一旦想回头查看某个基因在所有细胞中的原始表达情况从adata.raw取最稳妥。这行代码几乎是我每个项目必写的adata.raw adata.copy()如果你是从scanpy现成的教程数据集比如pbmc3k入门的这些步骤都有标准答案。但换成自己的真实数据各种边角问题就冒出来了。下面正式进入本篇的主题——降维、聚类、注释。2. PCA降维为什么单细胞分析绕不开这一步2.1 PCA不是为了出图而是为了去噪和压缩很多新手看到PCA第一反应是这玩意儿画出来又看不懂为什么不直接拿高变基因去做UMAP和聚类这个问题问得非常到位。PCA在这里的角色可以类比成给一堆嘈杂的声音做“降噪抽取主旋律”。高变基因有两千多个但大多数基因的表达模式是有相关性的。比如某个细胞亚群会同时上调20个基因这20个基因并不算20个独立信号本质上是一组共变模块。PCA做的事就是把这些共变模块拆成一个个互相正交的“主成分”每个主成分捕获的方差逐次递减。scanpy里跑PCA很简单import scanpy as sc sc.tl.pca(adata, n_comps50, svd_solverarpack)这里n_comps50的意思是先算50个主成分后面用几个另说。svd_solver有两种选项arpack适合矩阵规模中等、内存有限的情况数据量特别大的话可以换randomized速度更快但精度稍微有点损失——实际使用中我很少感受到两者结果有明显差异。2.2 该保留多少个主成分三个实用判断方法这是本环节最核心的问题。PCA算出来了50个主成分但后边聚类如果全部用上等于又把噪声带回来了。一般推荐的做法是先跑一个elbow plot看看拐点sc.pl.pca_variance_ratio(adata, n_pcs50) # 看前30个主成分累计解释了多少方差 print(cumsum_ratio[29])经验法则。基于我实际跑过的几十个样本大部分外周血、PBMC数据保留15-25个主成分效果最好组织样本因为细胞类型更混杂信号维度更多经常需要到20-35个。少于10个通常偏少说明数据里的异质性没被充分捕获多于50个往往过度拟合聚类容易把一些技术性伪影当成分群信号。2.3 PCA结果怎么快速检查跑完PCA别急着往下走。先看一眼主成分的热图判断前几个主成分是不是押在已知的细胞类型或批次信号上sc.pl.pca_loadings(adata, components1,2,3)这个图会展示每个主成分上贡献最大的那些基因。如果你做的是免疫细胞样本PC1上出现CD3D、CD3E这类T细胞标志基因或者LYZ、S100A8这类髓系标志基因那就说明主要变异信号确实是细胞类型差异。反过来如果前几个主成分主要贡献来自线粒体基因或核糖体基因一堆MT-开头和RPS/RPL开头那就要警惕了。这说明数据里残留的技术性噪音如破碎细胞比例过高还没有被完全清理干净可能得回头调整质控阈值。另一个快速检查维度是看PCA图上有没有非常明显的某个来源标签分群比如测序批次、样本编号、建库日期。如果分群太整齐往往提示数据中有强批次效应——这个问题在后边介绍批次整合的方法时再展开。3. 邻居图与UMAP/tSNE可视化背后的计算逻辑3.1 为什么scanpy要先构建邻居图再画UMAP很多人把UMAP当成一个普通的“降维画图函数”总觉得sc.tl.umap这一步是纯可视化的。这个想法不能说错但它会让你错过理解整个分析逻辑的机会。scanpy的做法是两步走sc.pp.neighbors(adata, n_neighbors15, n_pcs20) sc.tl.umap(adata, min_dist0.5)第一步pp.neighbors构建的是细胞间的“连接关系图”。它基于PCA空间默认取你指定的前20个主成分计算细胞间的距离然后用KNN算法找出每个细胞最近的15个邻居最终生成一个带权重的图结构——每个细胞是图上的一个点细胞间的连线权重代表它们的相似度。第二步tl.umap才是把这张高维图嵌入到二维平面。注意UMAP在这个阶段已经不再直接使用PCA坐标而是在邻居图的基础上做“图布局”优化。这也是为什么如果你改了n_neighbors或n_pcs哪怕PCA结果一模一样UMAP图也完全不同。3.2 n_neighbors和min_dist到底怎么调这两个参数是画图时最常被反复调整的。先说结论它们直接影响的是“分辨率”和“局部结构保留”的trade-off。n_neighbors控制UMAP在计算时考虑多大的局部邻域。值设得越小比如5-10算法越关注细胞之间的精细局部结构画出来的图往往更“碎”、更多小簇值设得越大比如20-30算法越关注全局结构画出来的图更“拢”、大群之间边界更清楚。min_dist控制最终图上点与点之间的最小间距。值越小点挤得越紧簇和簇之间的分离越明显值越大点排列越松散但可能牺牲一些局部结构的可读性。给一个参考起点目前的默认min_dist0.5其实对很多免疫细胞场景略微偏大。我在PBMC数据上常用的配置是sc.pp.neighbors(adata, n_neighbors15, n_pcs20) sc.tl.umap(adata, min_dist0.3)这个组合在区分T细胞亚群CD4 vs CD8naive vs memory时表现不错。如果你主要目标是看大类群比如五大免疫细胞谱系可以试试n_neighbors20, min_dist0.5画面会干净利落不少。3.3 tSNE什么情况下还能派上用场这个话题几乎每逢单细胞分析都会被拿出来对比。现在的实际状态是UMAP基本是主流tSNE被提起得越来越少。但不能因此就认定“tSNE没用”。tSNE有一个特性是所有数据可视化工具都比不上的——它对局部结构的保持极其敏感。同一个数据UMAP可能已经能把预想的细胞亚群分开tSNE则能把这些亚群内部的层次关系也显示出来。代价是tSNE计算极慢而且相空间里远距离细胞的几何关系几乎不具有任何解释性。我个人的建议是最终发表图用UMAP但探索数据阶段可以同时画一张tSNE做交叉验证。尤其是当你发现UMAP上某个亚群边界模糊、疑似中间态细胞很多的时候tSNE往往能告诉你这些中间态细胞到底更贴近哪一头。scanpy里tSNE的实现也很简单sc.tl.tsne(adata, n_pcs20, perplexity30) sc.pl.tsne(adata, colorcell_type)perplexity可以理解为tSNE的邻居数一般设在15-50之间。计算慢是老问题所以我现在基本只在关键样本上用。3.4 邻居图的中间产物怎么帮助判断还有一个容易被忽视但很有用的图sc.pl.neighbors(adata)会把构建出来的邻居图以非布局方式画出来主要用于观察图的连通性。如果图中存在大量完全孤立、不跟任何簇相连的“孤岛”往往是低质量细胞或doublet的体现。这类细胞在后边聚类时会单独扯出一个簇干扰注释。看到这种情况我通常会把这些孤岛标出来回到原始QC指标总counts、基因数、线粒体比例上检查。如果确认是低质量细胞可以考虑在聚类前再过滤一轮比聚类后手动删稳妥。4. Leiden聚类分辨率参数和结果评估4.1 从Louvain到Leiden聚类算法到底换了个什么scanpy早期版本教程里默认是Louvain算法现在新版本基本默认Leiden。这两个算法都是社区发现算法——把图里的节点细胞划分成内部连接紧密、外部连接稀疏的社团细胞群。Louvain在2019年被发现一个病理性的缺陷它可能产生“断开的社区”。意思是算法把一些明明没有联系的节点强行划分为同一社区。这个bug在单细胞场景下危害不小因为你可能看到一个“细胞群”里其实混合了两种完全无关的细胞。Leiden算法正是针对这个问题提出的改进它做了一步额外的细化refinement确保输出的社区是内部连通的。因此在scanpy里我从新版开始基本直接用sc.tl.leiden(adata, resolution0.8, key_addedleiden_0_8)注意key_added的作用把每次不同resolution的聚类结果存到不同key里。这样后面你想对比分辨率1.0和0.5的聚类差异只需要来回切换颜色标签不必重跑聚类。4.2 resolution参数别再用单一数值走天下了分辨率是聚类的核心超参数也是新手最困惑的地方。它的直观意义是分辨率越高聚类出来的群越细、越多分辨率越低群越粗、越少。但是这个“多少算合适”没有统一答案取决于你的生物学问题。举两个典型例子如果你研究的目标是“患者外周血里的免疫细胞组成”那分辨率0.5-1.0就够——它能把T、B、NK、Monocyte这些大群分出来。但如果你要深挖T细胞里的亚型Naive、CM、EM、EMRA、Treg分辨率经常得调到1.2-2.0不然这些亚型会被压在一个超大簇里看不出来。我通常的做法是跑一系列分辨率for res in [0.1, 0.3, 0.5, 0.8, 1.0, 1.5, 2.0]: sc.tl.leiden(adata, resolutionres, key_addedfleiden_{res})然后把leiden_0.8设定为主要分析标签再根据后续marker基因在各res下的表现决定要不要微调。这里的关键判断不是“哪个分辨率更漂亮”而是“哪个分辨率下每个簇都能找到相对清晰的marker基因身份”。4.3 聚类结果怎么看质量三个切入点聚类跑完不能直接去注释先做三件事。第一看各簇的细胞数和占比adata.obs[leiden_0.8].value_counts()理想情况下不应该有特别大比如占到总数30%以上的单簇也不该有特别小比如只有几十个细胞的碎簇。如果有别急着注释先回到数据检查。第二看PCA空间里聚类的分布。聚类是基于邻居图做的而邻居图又基于PCA空间。回头在PCA图上上色sc.pl.pca(adata, colorleiden_0.8)如果出现两个不同的leiden簇在PCA空间中几乎完全重叠说明聚类分辨率偏低或这两个簇之间的边界本来就很模糊。如果同一簇在PCA空间中明显断开成两块说明聚类偏碎。第三检查是否存在“批次标签”与聚类结果强相关。如果你有多个样本合并成一个adata用批次信息上色看看sc.pl.umap(adata, color[leiden_0.8, batch])如果每个leiden簇内都基本只来自同一个批次那说明聚类的驱动力很可能不是生物学信号而是批次效应。这种情况比较麻烦后续要么做批次整合如Harmony、scVI要么至少得把批次标签作为混杂因素在注释和差异分析时加以说明。4.4 聚类树状图和marker热图确定最终分辨率这是一个很多人不知道的小技巧。scanpy里可以用sc.tl.dendrogram对聚类结果做层次聚类看哪些簇在表达谱上更相似sc.tl.dendrogram(adata, groupbyleiden_0.8) sc.pl.dendrogram(adata, groupbyleiden_0.8)如果两个相邻的leiden簇在树状图上非常接近并且合并后marker基因模式也没有显著差异说明当前分辨率把它们分得过头了。反过来如果树状图上距离很远的两个簇被归在一起说明存在不合理的分裂。同时用marker基因热图辅助判断marker_genes { T: [CD3D, CD3E], B: [MS4A1, CD79A], NK: [NKG7, GNLY], Monocyte: [LYZ, S100A8], } sc.pl.dotplot(adata, marker_genes, groupbyleiden_0.8)哪个分辨率下每个候选簇在对应marker上出现干净的唯一表达模式就选哪个。这个流程比单纯盯着“多少个簇最好看”要有依据得多。5. 细胞注释从marker基因到自动注释一个都不能少5.1 先区分你是做“谱系注释”还是“状态注释”细胞注释看起来只是一步其实包含两个层次。谱系注释指的是确定细胞属于哪个大类T细胞、B细胞、NK细胞、髓系细胞、上皮细胞、成纤维细胞等。这是所有单细胞分析的第一步相对容易因为大类之间有非常清晰的marker基因。状态注释指的是在谱系确定后再细分比如T细胞里的CD4 vs CD8naive vs memory vs exhausted巨噬细胞里的M1 vs M2样。这一步难度大得多因为同一谱系内不同状态的表达差异往往是连续的没有绝对的黑白边界。两种注释策略完全不同。谱系注释可以依赖经典marker表 表达矩阵快速搞定状态注释则需要更细致的打分、单基因成像、甚至结合轨迹分析才能下结论。5.2 经典marker基因速查表人外周血/免疫组织实战中最常用的marker基因整理一份我自己的速查表细胞类型核心marker人补充markerT细胞泛CD3D, CD3ECD2, TRACCD4 TCD3D CD4IL7RCD8 TCD3D CD8ACD8BTregFOXP3, IL2RACD25CTLA4, IKZF2NK细胞NKG7, GNLY, KLRD1NCAM1CD56B细胞MS4A1CD20, CD79ACD19, BANK1浆细胞MZB1, SDC1CD138XBP1, IGHG1单核细胞经典CD14, LYZ, S100A8FCGR3A低单核细胞非经典FCGR3ACD16, MS4A7CD14低或无树突状细胞CLEC9A, BATF3cDC1CD1C, FCER1AcDC2LILRA4pDC巨噬细胞CD68, CSF1R, MRC1CD206C1QA, C1QB中性粒细胞FCGR3B, CSF3RS100A8/S100A9上皮细胞EPCAM, KRT8, KRT18KRT19成纤维细胞COL1A1, COL3A1, DCNPDGFRA内皮细胞VWF, PECAM1CD31CLDN5, FLT1这份表投射到scanpy里注释的做法是这样先用sc.tl.rank_genes_groups给每个leiden簇找差异基因然后手动检查每个簇的top差异基因是否落到上表的对应组合里。sc.tl.rank_genes_groups(adata, groupbyleiden_0.8, methodwilcoxon) sc.pl.rank_genes_groups(adata, n_genes20, shareyFalse)rank_genes_groups的输出会按每个簇列出差异最显著的基因默认按adjusted p-value排名。这里要小心一个常见误区差异表达排名第一的基因并不等于是这个簇最有特异性的marker。因为wilcoxon检验看的是差异显著性而显著性同时受效应量和样本量影响。一个高变基因在一个大簇里可能p值极小但实际表达量差异很小。所以更稳妥的做法是结合sc.tl.score_genes做一个基因集打分专门评估每个簇在核心marker组合上的得分sc.tl.score_genes(adata, gene_list[CD3D, CD3E, CD2], score_nameT_score) sc.tl.score_genes(adata, gene_list[MS4A1, CD79A, CD19], score_nameB_score) sc.pl.umap(adata, color[T_score, B_score])哪个簇T_score高而其他基因评分低它就大概率是T细胞。这个思路比单看一个marker要稳健得多因为单基因可能出现dropout某些细胞该基因没有捕获到表达但多个基因联合打分能大幅降低这种假阴性风险。5.3 自动注释工具怎么用SingleR和CellTypist的取舍手动注释虽然直观但到大样本几万甚至几十万个细胞时非常费人力而且不同人的注释标准未必一致。这时候自动注释工具是很好的加速器。SingleR是较早也较流行的参考集注释方法基于参考转录组数据集计算每个测试细胞与参考集中各细胞类型的相关性挑最匹配的标签。用法大致是import SingleR # 参考数据需要是已注释好的表达矩阵 ref_data ... ref_labels ... pred SingleR.annotate(testadata.to_df().T, refref_data, labelsref_labels)但SingleR有个问题就是它非常吃参考集的代表性。参考集里没有的细胞类型它也没法识别只能硬给一个最接近的标签。而且它基于「相关性」做判断对状态连续的高可塑性细胞如巨噬细胞极化中间态容易给错。CellTypist是近两年比较受欢迎的新选择基于深度学习模型在免疫细胞注释上准确率表现不错。它有个网站版也有Python API。使用逻辑比较简单from celltypist import models, classify model models.Model.load(modelImmune_All_High.pkl) predictions classify(adata.to_df().T, modelmodel)我觉得最务实的方案是**“自动工具初筛 手动复核”双轨制**先用CellTypist或SingleR给出候选注释然后在UMAP上用marker基因打分和rank_genes_groups验证每个标签是否站得住。或者反过来——先手动注释再用自动工具检查有没有漏掉的稀有种群。这两种路径我都用过推荐后者的次数更多因为手动注释能保留你对数据的直觉自动工具只作为安全网。5.4 注释结果怎么展示和落地注释最终会回写成一个列adata.obs[cell_type] adata.obs[leiden_0.8].map(cluster_to_celltype)映射表一般是一个字典比如{0: CD4 T, 1: CD8 T, 2: B, ...}。如果某个簇的注释你迟疑不决建议明确标成“未知”或“暂定”。硬要给出一个标签后续下游分析可能会被这个小错误污染。注释完成后的展示组合有很多。我常用的几组sc.pl.umap(adata, color[cell_type, leiden_0.8, CD3D, CD79A]) sc.pl.dotplot(adata, marker_genes, groupbycell_type) sc.pl.heatmap(adata, marker_genes, groupbycell_type)到这里一个单样本的“降维—聚类—注释”主流程就闭环了。至于差异表达、富集分析、细胞通讯这些下游内容全部建立在这个注释结果之上——注释错了后面全得推倒重来。所以哪怕前面几步跑得再快注释这一步也值得花最多时间反复确认。6. 这几个坑我在真实数据上踩过提前替你试试本章是我的私货时间全是从自己跑过的PBMC、肿瘤和自身免疫病样本里总结出来的经验教训不一定写在任何官方文档里但基本都会遇到。6.1 doublet混入导致“神秘新细胞类型”场景重现聚类注释时发现一个簇marker表达看起来很“杂”——既表达NK细胞的NKG7、GNLY又高表达T细胞的CD3D。在UMAP图上它往往位于T和NK两个大簇的中间地带。这就是典型的doublet两个细胞被一个液滴捕获特征。处理方式有两个方向分析流前方加scvi的SOLO或scrublet做doublet预测把高分doublet细胞直接过滤掉。或者在注释阶段看到这种“不伦不类”的簇别强行命名先把它标为“疑似doublet”后续单独剔除。我比较推荐两手都做前期过滤掉明显的后期聚类看到可疑的再用DoubletFinder或scrublet复核。6.2 高变基因筛选和降维的基因集不一致这个问题通常出现在你做了两轮高变基因筛选之后。比如第一轮筛选跑出了3000个高变基因后续你用adata.var[highly_variable]做了PCA然后你换了个参数又跑了一遍高变基因筛选。此时adata.var里的布尔标签被覆盖了但PCA结果还是基于旧的高变基因集合——两者对不上。如果后续别人或者三个月后的你拿到这个adata重新跑聚类用新的高变基因子集做PCA聚类结果会完全不同。scanpy本身不会帮你检查这点。我的习惯是高变基因筛选只做一次确定参数后记录在案不到万不得已不变更。如果确实要重跑严格重头开始不保留之前的PCA结果。6.3 注释时只信“新手套餐”marker很多入门教程都只给几个最泛化的marker比如CD3、CD19、LYZ这在PBMC数据里够用。但换成组织样本就远不够了。举个例子组织里的巨噬细胞和树突状细胞都高表达CSF1R和C1Q单看这两基因很难分开内皮细胞和某些成纤维细胞也都可能低表达PECAM1。经验法则是每个候选细胞类型至少准备3-5个marker其中至少包含一个正marker和一个负marker。比如鉴定CD4 T细胞不只要看CD4的表达还要排除CD8A、NKG7鉴定cDC2不能只看CD1C还要确认LILRA4阴性排除pDC。6.4 批次效应在聚类上表现“完美”但令人警惕有些样本做批次整合后UMAP看起来非常干净——每个批次完全混合在一起。但这时候要多个心眼过强的批次整合可能把真实的生物学差异也抹掉了。比如同一个个体治疗前后的细胞状态变化如果批次和生物学变量完全混淆治疗前批次1治疗后批次2那么批次整合模型会把治疗效应也当作批次效应给纠正掉。在没有明确实验设计能分离“批次”和“生物学差异”的情况下我通常先不着急整合而是分别跑单样本聚类注释然后用跨样本聚类比较的方式分析。只有确认批次效应成为分析瓶颈后才考虑Harmony或scVI。6.5 UMAP图的“距离”——别当它是真实的生物学距离新手最容易产生的误解是认为UMAP上两个簇离得远就代表它们生物学差异大。严格来说这个说法只在很粗糙的层面成立。UMAP做了非线性变换图上的欧几里得距离不能直接解读为转录组距离。所以当你在UMAP上看到某个过渡态的细胞群时别急着下结论说“这是分化中间态”建议用sc.tl.draw_graph基于ForceAtlas2布局或直接回到PCA距离做量化验证。6.6 聚类分辨率审稿人问起怎么答这算个加分项。如果投稿时审稿人问“为什么选resolution0.8”你不能只回答“试了几个觉得顺眼”。比较好的做法是给出一系列分辨率的聚类结果然后展示你选择特定分辨率时的marker基因表达一致性检查以及关键簇在新分辨率下是否保持稳定。这个稳定性验证的过程也是你自己对注释质量的一次压力测试。从PCA到聚类再到注释这条流程跑通不难跑好很难。难点不在函数调用而在每一步你都能为自己的参数选择给出生物学上的解释。降维时知道为什么选20个PC聚类时知道为什么选0.8的分辨率注释时知道为什么给这个簇冠上“CD8 resident memory T”的名称——这些“为什么”才是分析质量的分水岭。后面如果要继续扩展可以在这个基础上做细胞类型特异的差异表达、拟时序分析、细胞通讯推断。但前提都是本篇文章里这个基础打得足够扎实。祝顺利拿到数据别急着出图先把每个细胞的身份搞清楚。