单细胞测序的系列推文写到第五篇前面的数据质控、标准化、高变基因筛选这些“地基”工程大家都已经打好了今天终于到了整个流程里最出彩、也最让人兴奋的部分——把几万个细胞在二维平面上画出来然后找出每个群的“身份证基因”也就是我们常说的marker基因。这一步做完你手里的数据才算真正开始“说话”细胞类型注释、功能分析、拟时序分析等后续工作才能铺开。t-SNE聚类和marker基因寻找可以说是单细胞数据分析从“数据预处理”跨入“生物学解读”的转折点。很多刚接触单细胞的小伙伴做到这一步经常一头雾水RunTSNE、FindClusters、FindAllMarkers这些函数到底按什么顺序跑参数一堆到底怎么调跑出来的图细胞都糊成一团怎么办marker基因一大堆该怎么挑这篇文章我会把整个流程从原理到实操完整走一遍把每一步为什么这么做、参数怎么选、踩过的坑有哪些全都摊开讲清楚。这篇内容适合两类读者一类是已经跑通了上游流程、正准备做聚类注释的单细胞入门选手另一类是做了几次但总觉得结果不太对、想回头搞清楚原理的分析人员。我会以Seurat包为例主流程跑一遍t-SNE聚类和marker基因鉴定同时穿插一些R语言代码细节让大家可以直接在本地复现。1. 整体设计与思路拆解1.1 单细胞数据为什么要降维和聚类单细胞转录组测序scRNA-seq一次实验通常能拿到数千到数万个细胞每个细胞我们检测到的基因数在小几千到上万之间。如果把每个细胞看成一个样本每个基因看成一个维度那我们手里的数据就是一个“细胞数 x 基因数”的高维矩阵。拿10x平台的数据来说常见的矩阵规模是5000个细胞 x 20000个基因这个维度处理起来问题很大。第一个问题是计算量。直接在两万维空间里算细胞两两之间的距离矩阵大小是5000x5000还能忍但如果是十万个细胞十万乘十万的距离矩阵就是100亿个元素内存直接爆掉。第二个问题更重要——生物学信号被淹没在噪声里。两万个基因里真正能区分不同细胞类型的可能只有几百个剩下的大部分基因表达水平在各类细胞间差异不大它们的存在反而会稀释真正的差异信号。所以我们需要两步操作第一步是降维把高维数据压缩到几十个“主成分”里把最关键的信息提取出来第二步才是聚类在压缩后的低维空间里找到细胞之间的“天然分群”。t-SNE就是最常用的降维可视化方法之一配合聚类算法在Seurat里默认是Louvain算法使用能够在二维平面上直观展示细胞群体结构。1.2 为什么用t-SNE而不是直接PCA或UMAP很多新手会问PCA不是已经能降维了吗为什么还要t-SNE这个问题的关键在于PCA是线性降维t-SNE是非线性降维两者的目标完全不同。PCA追求的是“保留数据整体的方差”它擅长捕捉主要的线性结构比如不同处理组之间的整体差异但对局部结构不敏感。打个比方PCA像站在山顶看城市全貌能看清城市的大致轮廓和主干道但看不到街道里具体哪户人家。t-SNE则相反它更关注“相邻细胞的局部关系”目标是让高维空间中相近的点在二维平面上也靠得近原本离得远的点在图上就远。这就像走进街区虽然看不到城市全貌但每条巷子、每栋房子的关系清清楚楚。实际分析中两个都是用但分工不同PCA的结果作为t-SNE和聚类的基础t-SNE或者UMAP则负责最后可视化呈现。至于为什么不直接用UMAP我的看法是t-SNE在单细胞数据里被验证的时间最长社区积累的经验最丰富你的数据团分得开分不开、marker基因找得对不对很多老经验都是基于t-SNE图形总结出来的。UMAP速度快、全局结构保真度更高也可以作为补充展示但分析流程的主干用t-SNE更稳妥尤其是对新手来说遇到问题更容易找到参考资料。2. 核心细节解析与实操要点2.1 t-SNE的四个关键参数t-SNE算法本身并不复杂但参数对结果影响极大。很多人跑t-SNE时图不好看八成不是数据问题而是参数没调对。perplexity困惑度这是最核心的参数它控制算法认为每个细胞周围有多少个“邻居”。Seurat默认值是30但对单细胞数据来说30往往偏小。如果perplexity太小结果容易出现小碎群看起来像“撒了一把豆子”如果偏大则各群可能挤在一起分不开。我的习惯是先用30跑一遍如果图形过于破碎增加到50如果太糊就降到20一般三轮以内能找到合适值。需要注意perplexity值不应大于细胞总数的三分之一否则算法会报错。iterations迭代次数默认值是1000。t-SNE是通过不断迭代优化代价函数来获得最终布局的迭代太早停止图形还没稳定细胞排列会很乱看起来像一团噪声。如果跑了1000轮仍然不稳定可以从输出信息看到cost function是否还在下降可以增加到2000甚至3000。learning rate学习率Seurat里对应的是eta参数默认是200。学习率决定了每次迭代的步长值太大细胞会挤成一条线或塌缩成一团值太小收敛极慢。如果图看起来一团乱麻可以考虑把eta调小到100或者调大到500试试。theta近似度这是Barnes-Hut加速算法的参数。theta为0时是精确计算但速度很慢默认0.5时计算速度快但会有一定近似。十万级以上的细胞用0.5没问题小样本几千细胞建议设成0出图更精确。参数调起来确实有“手感”成分我的建议是不要追求一次到位先跑一版默认参数看一眼图再定向微调。这比从零摸索快得多。2.2 FindClusters的resolution到底怎么选聚类这一步的核心函数是FindClusters里面最关键的参数是resolution分辨率。这个值直接决定了你最终会得到多少个细胞群。Seurat的聚类流程是基于共享最近邻SNN图再用Louvain算法对图进行划分。resolution控制的是划分的“粒度”值越大分出来的群越多。默认值是0.8得到的群数通常在10~20个。对一般组织样本来说这个范围是合理的因为一个组织里常见的细胞类型也就十几种。但具体到自己的数据你需要根据marker基因的结果反复试验。分辨率太小不同类型的细胞可能被合并成一个群比如T细胞和NK细胞表面marker基因表达模式接近cluster分辨率不足时很容易被划到同一群里分辨率太大同一个细胞类型会被人为拆成好几个群增加了下游注释负担。我的经验是做一个“分辨率扫描”分别用0.4、0.6、0.8、1.0、1.2跑一遍把结果导出来对比每个群的marker找出那种“既能把已知细胞类型分开、又不会把同一类细胞拆碎”的分辨率。虽然多花点时间但这一步做好了后面注释会顺利很多。2.3 marker基因的统计检验逻辑找marker基因的标准做法是FindAllMarkers但很多人只关注返回了哪些基因却不清楚这些检验结果的生物学含义。这里有个重要概念marker基因的判定是靠统计检验完成的。Seurat默认使用Wilcoxon秩和检验非参数检验比较目标群细胞的基因表达量与其他所有细胞的差异。输出结果里有几个关键列avg_log2FC表示两组间表达差异的log2倍数值p_val是检验的显著性p值p_val_adj是多重假设检验校正后的p值pct.1和pct.2分别是该基因在目标群和其余群中检测到的细胞比例。这里特别要注意的是p_val_adj。因为每个基因都做一次检验两万个基因就要做两万次在这么多重检验下p值本身几乎没有参考价值必须看校正后的p_val_adj。我一般只保留p_val_adj 0.05且avg_log2FC 0.58也就是表达倍数变化大于1.5倍的基因作为候选marker同时要求pct.1和pct.2的差距足够大。为什么要求这个差距简单说如果一个基因在目标群80%细胞里表达但在其他群67%细胞里也表达即便差异检验显著用它来定义细胞类型也不可靠因为特异性不够。3. 实操过程与核心环节实现3.1 从Seurat对象到PCA降维为了让大家有完整的参照我先假设你已经拿到了质量控制后的Seurat对象名字叫scobj。这个对象里有原始计数矩阵、经过标准化后的数据以及高变基因信息。如果你是从上游直接跑过来的代码会类似这样library(Seurat) scobj - NormalizeData(scobj) scobj - FindVariableFeatures(scobj, nfeatures 2000) scobj - ScaleData(scobj)单细胞分析的后续步骤如PCA、聚类、t-SNE等通常只针对高变基因进行计算。为什么要这样因为高变基因是细胞间差异最大的那批基因它们承载了区分细胞类型的主要信息而大量在几乎所有细胞里表达水平一致的“管家基因”在分群中起不到区分作用反而会拖慢计算、稀释信号。接下来跑PCAscobj - RunPCA(scobj, npcs 50, features VariableFeatures(scobj))RunPCA之后怎么决定用多少个主成分这个问题的答案直接决定了后续聚类和t-SNE的效果。如果主成分数太少会丢失稀有细胞类型太多又会引入噪声让聚类结果变乱。Seurat提供了两种辅助判断方式一是ElbowPlot(scobj)会画出每个主成分解释方差的百分比找“拐点”——曲线急剧下降后趋于平缓的位置一般取拐点之前的主成分数量二是JackStraw方法通过置换检验找出显著的主成分对小数据集更可靠但计算量大一些。实际操作中我用PC数一般在10~30之间。十来个主成分的时候t-SNE图的每个群还相对松散增加到二十个以上分群更清晰稀有细胞群也更容易被识别出来。如果数据集的细胞类型比较丰富比如免疫细胞、基质细胞、上皮细胞共存的肿瘤样本我建议往高里选用20或更多。确定好PC数后就进入主题流程。3.2 t-SNE聚类分析批量跑法在Seurat里面聚类和t-SNE的代码并不长但先后顺序有讲究必须先聚类再跑t-SNE。因为FindClusters是基于PCA降维后的数据在“原来的高维空间”计算细胞间的SNN图而t-SNE只是将聚类结果可视化到二维平面它不会影响聚类结果本身。# 先聚类resolution可以先给0.8 scobj - FindNeighbors(scobj, dims 1:20) scobj - FindClusters(scobj, resolution 0.8) # 再跑t-SNE scobj - RunTSNE(scobj, dims 1:20, perplexity 30, seed.use 42) # 可视化 DimPlot(scobj, reduction tsne, label TRUE)这里有个细节FindNeighbors的dims参数需要和RunTSNE的dims一致否则你聚类和可视化用到的数据不是同一套逻辑上就不对。另外t-SNE算法的随机性很强每次跑结果可能略有不同所以一定要设置随机种子。我用的是seed.use 42这个数字本身没有特殊意义关键是固定下来保证结果可复现。跑完之后用DimPlot画图正常情况下你会看到若干个细胞团。这里我要提醒一句第一眼看到t-SNE图觉得“很不错”或“很糟”都没关系因为t-SNE是可视化工具不是结果判定的唯一标准。真正判断分群好坏要靠marker基因这个后面详细说。值得注意的是Seurat还有一个函数叫RunTSNE默认返回的是二维结果如果你需要三维t-SNE图可以设置dim.embed 3然后用DimPlot绘制。三维图在某些场景下确实信息量更大比如有多个相似的细胞亚群相互嵌套时二维投影可能会把它们叠在一起而三维可以帮助区分。不过目前主流发表文章还是以二维展示居多三维图更多用于内部辅助判断。3.3 聚类效果评估与可视化组合聚类跑完之后第一步不是急着找marker而是先做一轮“人工检查”看看t-SNE群和细胞数量的分布。# 看看每个群有多少细胞 table(Idents(scobj)) # 看看有没有某两个群离得特别近 DimPlot(scobj, reduction tsne, label TRUE, label.size 5)这一步能帮我发现几个潜在问题。比如某个群细胞数极少比如只有几十个这可能是一个稀有细胞类型也可能是数据质量差导致的假群需要后续重点检查再比如两个群在图上几乎是合并状态边界完全模糊那可能是分辨率太低也可能是这两个群本身就是同一种细胞的连续状态。此时我会去检查它们之间是否有差异表达的marker基因再决定是否需要提高分辨率重新聚类。同时建议用FeaturePlot把已知的细胞类型marker画在t-SNE图上做一次“快速对表”。比如你要分析的是外周血样本那就可以检查CD3DT细胞、CD14单核细胞、MS4A1B细胞、NKG7NK细胞等经典的marker在t-SNE上的表达分布。如果某个已知marker只在一个群或者少数几个群中高表达并且图案清晰边界分明说明聚类结果可信如果marker到处都是或者多个群都强表达那就要警惕聚类过度或数据混入异常细胞。FeaturePlot(scobj, features c(CD3D, CD14, MS4A1, NKG7), reduction tsne, ncol 2)FeaturePlot画出来的点图上每个点代表一个细胞颜色深浅代表该基因表达量的高低。在t-SNE图里基因高表达的细胞如果呈现出“两头翘”的分布通常是正常的但如果高表达细胞分散在各个群、没有明显聚集这个基因可能不是一个好的群marker需要谨慎使用。3.4 FindAllMarkers与FindMarkers的实操代码当你对聚类结果满意之后进入找marker基因环节。最常用的函数是FindAllMarkers它会遍历每个cluster把该群与其他所有群进行比较找出每个群的特异性高表达基因。# 找到每个cluster的marker基因 all_markers - FindAllMarkers(scobj, only.pos TRUE, min.pct 0.25, logfc.threshold 0.58)only.pos TRUE表示只保留在目标群中上调的基因往下游注释通常只关注阳性的特异性marker所以这个参数建议打开。min.pct 0.25表示基因至少在25%的细胞中检测到才会被纳入分析太低会把很多零散表达的基因也算进来加大噪声。logfc.threshold 0.58是差异表达的阈值对应1.5倍表达变化低于这个阈值不算显著富集。跑完会发现返回的数据框很大每个群有几十到几百个候选marker。接下来要做的是按p_val_adj和avg_log2FC排序挑出每个群的前几个marker用于注释library(dplyr) top_markers - all_markers %% group_by(cluster) %% arrange(p_val_adj, desc(avg_log2FC)) %% slice_head(n 10)FindMarkers和FindAllMarkers的区别是FindMarkers用于比较指定的两个群或一个群与另一个群的差异FindAllMarkers则是一次性做所有群的循环比较。实际分析中当你想比较两个相近亚群比如CD4 T和CD8 T的差异或者需要将一个群与特定参照组做比较时用FindMarkers更灵活。# 比较cluster1和cluster2 markers_1_vs_2 - FindMarkers(scobj, ident.1 1, ident.2 2, only.pos TRUE, min.pct 0.25)3.5 从marker基因列表热图到细胞类型注释marker基因列表出来后几种常用的可视化手段可以把结果表达得更加直观清晰。热图是最直观的用DoHeatmap把每个群top marker的表达量画出来可以一次看清每个群“专属高表达”的基因组合。DoHeatmap(scobj, features top_markers$gene, size 3)热图的每一行是一个基因每一列是一个细胞默认按cluster排列。理想状态下你应该看到沿对角线方向的明显“方块状”色块也就是每个群的marker在自身群体的细胞中集中高表达。如果热图上各群之间没有明显的颜色区分marker的表达不够特异就需要回去调整聚类参数。热图看完就到了真正的生物学判断环节——根据marker基因的已知功能与文献信息判断每个cluster的细胞类型。这一步是整个流程中最需要生物学知识沉淀的地方也是很多新手最头疼的部分。比如cluster 0高表达CD3D、CD3E基本可以判定是T细胞cluster 1高表达MS4A1、CD79A大概率是B细胞cluster 2高表达LYZ、CD68可能是巨噬细胞或者单核细胞。但有些群没这么友好比如高表达CCL5、NKG7、GNLY的细胞可能是NK细胞但也有可能是活化后的CD8 T细胞。这时候我会结合更多marker比如CD8A是否阳性如果CD8A也是阳性的那就进一步倾向于细胞毒性T细胞而不是经典NK必要时还需要借助权威的文献列表来确认。注释还有一个技巧千万不要只看top1的marker就下结论至少要看3~5个marker联合判断。比如一个群高表达CD3D和CD4但同时也高表达FOXP3和CTLA4免疫抑制相关基因那这个群很可能是调节性T细胞而非普通CD4 T细胞两者在功能上差异巨大不能混为一谈。4. 常见问题与排查技巧实录4.1 t-SNE图分不开或挤成一团的排查这是被问得最多的问题。t-SNE图分不开可能的原因有很多我把排查过程整理成一个速查表。表现可能原因处理办法所有细胞挤成一团几乎没有结构PCA主成分数太少增大dims从10调到20-30再试细胞群破碎严重出现大量小群perplexity太小增大perplexity到40-60群与群边界模糊互相渗透分辨率偏低增大resolution或用更高分辨率重新聚类图上出现明显的“直线”或“弧形”排列学习率eta偏大将eta从200下调到100同一样本不同次运行结果差异大未固定随机种子设seed并保持一致需要特别说明的是“挤成一团”还有一种常见原因是数据里存在大量的重复细胞比如同一个细胞被捕获了两次或细胞周期效应没有消除。这种情况不是单纯调参能解决的可能需要回到上游做去重处理或用ScaleData时加入细胞周期评分作为协变量回归掉。4.2 FindAllMarkers结果为空或top基因不特异偶尔会得到完全空结果的FindAllMarkers输出或者找到的标记基因其实在其他群体里也广泛表达。这里有个隐藏的“坑”logfc.threshold和min.pct这两个参数是互相关联的如果设置太严格比如logfc.threshold设成2对应4倍表达差异不少细胞群的marker基因可能全被过滤掉。此时可以降低到0.25~0.58试试。另一个常见情境是你分析的数据包含多个细胞类型差异并不大的亚群比如CD4 T细胞不同亚型它们之间的表达差异本身就是渐进式的。这时候不必强求找到“干净”的marker可以更关注基因组合或特异度较高的调控因子。我还遇到过一种情况某些cluster的top marker居然是线粒体基因或者核糖体基因比如MT-、RPS/RPL开头。这通常意味着数据质量有问题高比例线粒体基因往往代表细胞凋亡或损伤核糖体基因高表达可能是细胞应激的产物。遇到这种情况最稳妥的做法是回到质控阶段重新过滤而不是强行用这些基因作为细胞类型marker。4.3 从marker基因到细胞类型注释的实操心得最后分享一下我在实际项目里积累的几个注释心得。第一注释要有参考依据。最可靠的参考是已发表的单细胞文献和人/小鼠细胞图谱。比如做人的外周血样本可以参考Human Cell Atlas人类细胞图谱计划的相关数据做肿瘤样本可以搜索该癌种的单细胞图谱论文。还有一个更直接的办法如果你的样本与某个公开数据集接近直接下载公开数据的Seurat对象用FindTransferAnchors做label transfer标签迁移可以快速获得初步的注释结果再人工校正。第二善用R包的自动注释工具做初筛但不要盲信。目前比较成熟的有SingleR、Garnett、scCATCH等它们基于不同参考数据库给出细胞类型预测。对于初学者先用SingleR跑一版注释结果与手动注释对照能节省不少时间。但自动注释的准确率和参考库的质量高度相关不能完全替代人工判断。第三marker基因做到“闭眼能说出来”才算合格。我自己的习惯是每个注释好的细胞群至少总结3~5个特异性强、功能清晰的marker基因把它们背下来之后画FeaturePlot、检查群间相似性时脑海里立刻能涌现出对应的表达模式。这一步熟练之后速度会越来越快也不容易出低级错误。第四注意物种差异。人、小鼠的基因命名规则不同人源基因全大写如CD3D小鼠源基因首字母大写其余小写如Cd3d。如果你用人的marker列表去套小鼠数据代码里features参数写不对FeaturePlot会画不出任何东西而且低版本Seurat还可能直接报错。遇到这种问题先查Gene Symbol大小写这也是我踩过最多次的坑之一。4.4 给分析流程收尾的几个建议到这里t-SNE聚类和marker基因寻找的主流程就算完整跑完了。按照惯例最后补充几个新手上路最容易忽略的点。算力资源方面如果你的样本量很大超过5万细胞t-SNE会比较耗时需要预留足够的内存和运行时间。这种情况下可以考虑先降到10万细胞再跑t-SNE用于初步探索最终分析再用全量数据。GPU加速方面目前Seurat的t-SNE和UMAP实现并没有利用GPU大规模数据建议改用Python的scanpy配合leiden聚类速度优势明显。但这是另一套体系初学者先把手头流程跑透更重要。临时保存中间结果。找到合适的resolution和perplexity参数后建议用saveRDS把Seurat对象保存下来后续继续做拟时序、细胞通讯或其他分析直接在保存的对象上展开不需要重复跑聚类和marker寻找。我自己现在做数据每跑完一步t-SNE都会顺手把当次参数的图形导出保存哪怕之后的参数更理想也保留旧图备查。这个习惯帮我复盘过不少次也让最终选定的结果更有说服力。所有分析与判断都记录得清清楚楚将来写论文或者合作者复盘时每一张图、每个参数都有据可查这对科研数据管理的意义不言而喻。t-SNE聚类这步做完单细胞分析最难的“从数据到生物学”的跨越就算完成了一大半。剩下的事情无论是细胞类型鉴定、差异分析、拟时序分析还是细胞间通讯都是在这个基础上往上搭建的。只要这一步的基础打牢了下游哪怕再复杂也都有了可靠的前提。