单细胞聚类详解:Seurat的FindNeighbors与FindClusters原理与调参指南
2026/9/17 17:59:47 网站建设 项目流程

1. 为什么说聚类是单细胞分析绕不开的关口

只要做过单细胞转录组数据分析,就一定绕不开FindClusters这一步。拿到表达矩阵之后,我们面对的是几万个细胞、两万多个基因的庞大表格,单纯靠肉眼或者看某个基因的表达量根本没有办法判断数据里到底有哪些细胞类型。聚类分析本质上就是把“表达谱相似”的细胞归到一起,让每个cluster对应一种细胞状态或细胞类型——这一步的质量,直接决定了后面所有差异分析、拟时序分析、细胞通讯分析靠不靠谱。

Seurat的聚类流程被封装成了两个连续的函数:FindNeighborsFindClusters。很多初学者直接把这两行代码当成“标准动作”,跑完就急着去看UMAP图,结果换一个数据集、换一个resolution参数,得到的细胞群数量就完全不一样了,于是开始怀疑是不是自己的数据有问题。我自己带过不少学生和合作者,遇到这种困惑的次数特别多。其实问题往往不在数据本身,而是你没有搞清楚这两个函数在数学层面到底做了什么。

这篇内容我就围绕FindNeighborsFindClusters把原理拆开讲。我会先讲清楚Seurat为什么要先降维再聚类,然后分别解释这两个函数内部的算法逻辑,最后结合实际操作讲参数怎么调、常见坑怎么避开。内容偏原理但不会堆公式,尽量用直观的方式把“为什么这么做”讲明白,适合正在使用Seurat做分析、想深入理解聚类机制而不是只满足于跑通代码的读者。

2. 聚类之前的降维:为什么非得先跑PCA

在真正进入FindNeighbors之前,有一个前置步骤经常被忽略,那就是PCA降维,以及dims这个参数的选择。很多人不理解:FindNeighbors的输入本质上是一个细胞与细胞之间的距离矩阵,为什么不能直接用原始表达矩阵算距离,非要先降维?

原因有二。第一,原始表达矩阵的维度是细胞数乘基因数,一般都有两万个以上的基因维度,在高维空间里计算距离会出现“维度灾难”。简单说,维度一旦高了,所有细胞对之间的距离都会趋向于接近,距离的区分度变得很差,聚类算法很难找到有意义的近邻关系。第二,单细胞表达矩阵非常稀疏,包含了大量噪声,比如测序深度差异、批次效应、dropout事件造成的零值膨胀。这些噪声如果直接参与距离计算,会把真正的生物学信号淹没掉。

PCA的作用就是在这种情况下做一次“信号浓缩”。它把两万多个基因的表达模式压缩成几十个主成分,每个主成分都代表一组协同变化的基因模块,也就是一种潜在的表达程序。比如说,某些基因一起高表达可能是因为它们都属于T细胞激活程序,另一些基因一起变化可能是因为它们都受同一个转录因子调控。PCA把这些信息提取出来,保留前几十个主成分,就相当于只保留了数据里最稳定、最具生物学可解释性的信息,把细碎的噪声留在了后面被截断的分量里。

这也解释了为什么dims这个参数值得认真选而不是默认填一个1:10就完事。如果dims选得太小,比如只选前5个主成分,可能会丢失稀有细胞类型的信息;选得太大,比如选了30个,又会把噪声带回距离计算中。常规的做法是先跑ElbowPlot,看主成分方差贡献率的拐点落在哪里,再结合实际关注的目标细胞类型来定。有时候我也会参考JackStraw的显著性检验结果,但说实话,在绝大多数分析场景下,拐点图加领域判断已经足够用了。

选好dims之后,FindNeighbors拿到的就是每个细胞在PCA空间里的坐标矩阵。在PCA空间里计算距离,信息密度高、噪声低,这才是整个聚类流程能“聚得出来”的前提。

3. FindNeighbors在做什么:从PCA坐标到细胞关系网络

3.1 KNN图:先找到每个细胞最像的K个邻居

FindNeighbors的第一步是在PCA空间里为每个细胞找到距离最近的K个邻居,这个K由k.param控制,默认值是20。距离度量默认是欧氏距离,脚本里如果你不额外指定distance.matrix参数,用的就是欧氏距离。

所谓KNN图,可以理解为“给每个细胞发一张朋友名单”:每个细胞只跟自己最近的20个细胞建立连接关系,形成一个邻居圈。这一步的计算量并不小,几万个细胞两两计算距离的复杂度是O(n²),所以Seurat内部用了RANN包做近似最近邻搜索,速度上有很大的优化。实际跑几万细胞的时候,这一步通常只需要几十秒到几分钟,不会成为瓶颈。

这里有一个很容易被忽略但对后续聚类影响很大的细节:K值的设定决定了对数据局部结构的敏感程度。K比较小,比如10,每个细胞只连接少数几个邻居,网络会比较稀疏,对局部差异敏感,容易把小众细胞群单独聚出来,但也容易把同一个细胞类型因为微小的异质性拆成多个小cluster。K比较大,比如50,网络变密,聚类结果会更偏向大群结构,稀有细胞群很容易被“拉进”大群里消失不见。实际操作中,我一般会先用默认的20跑一遍主流程,如果发现稀有细胞群始终不出来,再尝试把K降到10或者15重新聚类,看结果是否稳定。

3.2 SNN图:为什么要给共享邻居加权重

如果只做KNN,聚类算法面对的是一个无权图,所有连接的重要性一视同仁。但生物学的直觉告诉我们:如果两个细胞不仅互相认识,而且共享了一大堆共同的朋友,那它们之间的关系应该比只有一条直接连接的细胞对更加紧密。SNN(Shared Nearest Neighbor)就是在这种直觉上建立起来的。

FindNeighbors默认在KNN基础上自动计算SNN矩阵。它的核心思想是:两个细胞之间的权重,取决于它们共享了多少个K近邻。共享的邻居越多,权重越高。Seurat官方的FindNeighbors文档里有一个nn.method参数,可以选择"rann""annoy",但无论选哪种近邻搜索方法,后面构建SNN的权重计算逻辑都是一致的。

SNN的引入有非常实际的生物学意义。单细胞数据里存在大量技术噪声,有些细胞虽然是真正的同类,但因为测序深度等因素,直接表达谱距离反而比跟异类细胞还要远。KNN对这种噪声是脆弱的,因为只要距离近就建边。而SNN通过“朋友的共识”来加权,两个细胞即使直接距离稍远,只要它们共享的邻居多,依然会被赋予高权重,这就相当于把局部的结构信息引入了图里,让聚类结果对噪声更稳健。

FindNeighbors的返回值包含两个矩阵:一个是RNA_snn,一个是RNA_nn。前者是加权后的SNN图,FindClusters实际上用的就是RNA_snn这个矩阵,后者是纯KNN邻接矩阵,主要用于后续可视化展示细胞间连接的时候使用。知道这一点之后,你再看FindClusters源码或者帮助文档,就不会再有“为什么聚类用的是snn而不是nn”的疑惑了。

3.3 从矩阵到图:找社区的前提是把数据变成网络

到这里,FindNeighbors做的事情可以总结成一句话:把细胞表达谱矩阵变成一个加权图。图中的节点是细胞,边是细胞间的近邻关系,权重代表关系的亲疏。图构建好之后,FindClusters的使命就非常清晰了——在这个图上做“社区发现”。

做一个类比来帮助理解:把每个细胞想象成社交网络里的一个用户,KNN是每个人主动加好友,SNN是系统根据共同好友数量给每条好友关系打分。最终的社交网络图谱里,兴趣爱好相同的人会形成密集的小团体,这些小团体就是细胞类型。FindClusters要做的事情,就是设计一套规则来“划分”出这些小团体。

4. FindClusters的核心算法:模块度优化与Louvain/Leiden算法

4.1 模块度是什么:比“谁和谁近”更高一层

FindClusters默认采用Louvain算法(新版本Seurat v5也支持Leiden算法),两者都属于模块度优化类算法。要说清楚这两个算法,先得把“模块度”这个概念讲明白。

模块度的英文是modularity,它衡量的是一个社区划分的质量:划分出来的社区内部边足够密集,而社区之间的边足够稀疏。公式长这样:Q = (1/2m) * Σ[A_ij - k_i*k_j/(2m)] * δ(c_i, c_j),其中A_ij是节点i和j之间的边的权重,k_i是节点i的总连接强度(所有边的权重之和),m是所有边的总权重,δ(c_i, c_j)表示i和j是否被分在同一个社区。

不用硬记公式,只需要抓住它的核心直觉:如果两个细胞之间的实际连接权重A_ij,明显高于随机情况下期望的连接权重k_i*k_j/(2m),那么这两个细胞放在同一个社区里,是对模块度有贡献的;反之,如果它们之间的连接比随机期望还弱,却硬被分在一起,模块度就会下降。所以Louvain算法的目标,就是找到一个划分,让总的模块度尽可能大。

这个“与随机期望比较”的设计非常妙。它本质上是在回答一个问题:这两个细胞的连接紧密程度,是显著超出了偶然水平,还是仅仅因为这两个细胞在整体网络里本身就很“活跃”?单细胞数据里,高表达基因多的细胞天然更容易跟其他细胞产生高权重连接,如果不做这种随机期望校正,这些“活跃”的细胞就会被错误地聚到一起。模块度的计算天然规避了这个陷阱。

4.2 Louvain算法的两步循环:局部贪心加全局聚合

Louvain算法的实现思路非常直观,分为两个阶段循环迭代。

第一阶段是局部移动。开始时,每个节点都被视为一个独立的社区。算法从左到右遍历所有节点,尝试把每个节点移动到它邻居所在的社区中,计算移动后模块度增量ΔQ是否大于0。如果大于0,就采纳这个移动,否则保持原状。如此反复遍历,直到任何节点的移动都无法再提升模块度为止。由于这个阶段的决策只看局部信息,算法跑得非常快,几万个节点通常几十秒就能收敛。

第二阶段是网络聚合。把第一阶段得到的社区视为新的“超级节点”,社区之间的连接权重等于所有跨社区边的权重之和,社区内部的连接也做相应的折叠,然后在这个压缩后的新网络上重新执行第一阶段。重复这两个阶段,直到模块度不再提升。

Louvain算法的优点是快、内存占用低,特别适合几十万甚至上百万细胞的数据集。但它有一个众所周知的缺点:分辨率限制,即无法识别出规模小于某个阈值的社区。这个阈值跟网络的规模有关,网络越大,能识别的最小社区规模也越大,这就是为什么大样本数据里一些稀有细胞类型特别容易被Louvain“吞掉”。

4.3 Leiden算法:解决Louvain的“社区连通性”问题

Seurat v5开始原生支持Leiden算法,通过algorithm参数指定为algorithm = 4即可调用。Leiden算法本质上是对Louvain的改进,在局部移动阶段增加了一个“细化子社区”的步骤:先把每个社区内部进一步划分为更小的子社区,再基于这些子社区决定是否移动节点。这个操作保证了最终划分出的社区是内部连通的,不会出现Louvain常见的一种问题——把一个实际上内部互不相连的松散群体硬包成一个社区。

Leiden方案的另一大优势是能保证聚类的连通性约束。有些社区的边界是模糊的,Louvain会基于贪心策略把它们归拢在一起,而Leiden会在细化步骤里把这些模糊连接拆开,从而保留更真实的群落结构。如果你的数据里存在发育轨迹这类连续性较强的细胞状态,比如从干细胞到分化终末阶段的连续谱系,Leiden往往能比Louvain给出更符合生物学直觉的划分。

对于小数据集或者对稀有群体特别关注的场景,我建议试试Leiden。它比Louvain的计算开销稍高,但现代机器的算力基本可以忽略这个差距。选择哪一种,关键还是看你的科学问题:如果在做免疫微环境的稀有亚群挖掘,Leiden通常更合适;如果只是常规的细胞注释、做个大群划分,Louvain已经足够用了。

4.4 resolution参数的本质:给模块度加一个“放大镜”

FindClusters里被问得最多的问题就是:resolution到底应该设多少?0.5还是1.0?要回答这个问题,得先理解resolution在算法里起什么作用。

Seurat的模块度函数里引入了resolution参数γ,实际优化的目标是 Q = (1/2m) * Σ[A_ij - γk_ik_j/(2m)] * δ(c_i, c_j)。γ乘在随机期望项上,它的作用相当于调节“多少倍于随机期望的连接强度才算是有效的社区内连接”。

当γ小于1时,随机期望项被缩小,相当于放松了社区划分的阈值,算法更容易把松散的节点合并成较大的社区,所以总cluster数会变少。当γ大于1时,随机期望项被放大,节点之间的连接需要比随机期望高出更多倍才会被归入同一个社区,于是小社区更倾向于保持独立,cluster数量变多。这就是为什么resolution从0.5调到1.2,UMAP图上cluster数量往往显著增加。

理解了这个机制,你就可以摆脱“套默认参数”的焦虑了。我通常的做法是跑一个resolution梯度,一般是从0.1、0.3、0.5、0.8、1.0、1.5这样递增,然后结合已知的marker基因来评估哪个分辨率下的cluster与已知生物学最吻合。有一个常用的辅助工具叫clustree,可以可视化不同分辨率下cluster的稳定性:能够跨多个分辨率稳定维持的cluster,大概率是真实存在的细胞类型;而总是来回分分合合的,往往是边界不清的过渡态细胞群。

提示:不同分辨率的聚类结果之间,没有“对”与“错”的绝对标准,只有“是否适合回答你的科学问题”这个相对标准。做差异分析时,偏大的resolution会给你更多细分亚群;做细胞注释时,偏小的resolution更容易给出清晰的、可命名的大类。

5. 实操参数组合与常见坑

5.1 一套可以参考的参数梯度策略

直接给一套我常用的参数策略,适合10x标准建库的PBMC或肿瘤组织样本:

分析目的k.paramalgorithmresolutiondims参考说明
常规大类注释201(Louvain)0.5依据ElbowPlot得到5~10个大群,便于marker注释
稀有亚群挖掘10~154(Leiden)0.8~1.2可适当多留几个PC对小群体更敏感,cluster更细碎
发育轨迹分析20~304(Leiden)0.3~0.5尽量少的PC保留连续过渡态,不要太碎裂

dims的选择是这些参数里最容易被忽视的一环。很多人习惯了直接用1:20或者1:30,并没有去检查到底有多少个PC承载着有意义的信号。如果PCA结果里第10个PC之后基本就是噪声,那1:30实际上就是把这些噪声当成了聚类依据,最终表现为cluster分不开、注释不清晰。反过来,如果把dims设得太小,一个稀有细胞亚群的信号恰好落在后面的PC里,就会被丢掉。所以每次拿到新数据,我都建议先跑一遍ElbowPlot,再跑一遍DimHeatmap看感兴趣PC的基因载荷是否具有生物学意义,再决定dims的取值,这个步骤花不了几分钟,但对后续聚类质量的提升非常明显。

5.2 聚类前必须检查的3个前置条件

很多人聚类跑完发现结果离谱,回过去查才发现是前面的数据质控就出了问题。这里有三个我踩过坑之后形成的检查习惯,列出来给大家参考。

一是PCA之前是否做了正确的数据标准化。Seurat的NormalizeData默认是LogNormalize方法,也就是log1p(counts/total_counts * 10000)。这一步必须在ScaleData之前做,顺序搞反会导致后续所有分析失真。ScaleData默认只对VariableFeatures里的基因做中心化和标准化,这一步的默认行为本身没问题,但要注意它默认回归掉的只是测序深度相关的变异,如果你知道数据里有强的批次效应,得在这里额外设置vars.to.regress,而不是逃避做一个更严谨的批次整合。

二是是否有明显的批次效应。如果你是把多个样本合并在一起分析的,建议聚类之前用Harmony或者Seurat自带的IntegrateData做批次校正。很多人忽略了一点:FindClusters聚类时的assay默认是RNA,如果你做了整合分析,生成的是integratedassay,必须在FindNeighbors里显式指定reduction = "pca",并且保证这个pca是基于integrated数据计算的。在Seurat v5里,整合后的默认reduction通常可直接使用,但旧版本里经常有人因为没指定reduction,导致聚类还是基于未校正的RNA数据来做的,批次效应直接带进了聚类结果。

三是双细胞的比例是否过高。FindClusters对双细胞很敏感,尤其是一些表达谱广泛、转录本量高的细胞类型(比如巨噬细胞、肿瘤细胞),经常会被错误地聚成一个高转录本的“垃圾群”。我一般会在聚类之前用DoubletFinder或者Scrublet跑一遍双细胞预测,把明显的高分双细胞先过滤掉,而不是等聚类之后再去猜测哪个群是双细胞混合群。

5.3 聚类之后如何验证结果是真的“对”

聚类做完之后不要急着往下游走,先做两个验证。

第一个是marker基因的验证。用FindAllMarkers跑出每个cluster的差异表达基因,然后跟已知的细胞类型marker对照。一个合格的结果是:每个cluster至少有一个明确的marker基因组合能支持它的身份注释。如果出现一个cluster的marker基因列表里全是rRNA、线粒体基因或者热休克蛋白,那这个cluster大概率来自低质量细胞,应该考虑回溯到质控环节去检查。

第二个是umap图的目视检查。UMAP只是用于可视化,不是聚类依据,但它能直观反映聚类的结构是否合理。我关注两件事:一是有没有“被强行拼接”的哑铃状cluster——两头是两种不同的细胞类型,中间只有少数几个细胞勉强连接,这种往往说明分辨率太低或者K值太大,需要用更高分辨率重新聚类;二是同一个已知细胞类型是否被拆成了很多小碎块——如果是多个cluster都表达同样的marker,差别只在于一些增殖相关基因的表达高低,那很可能不是真正的亚型,而是细胞周期的影响,需要在ScaleData时回归掉细胞周期相关基因,或者用CellCycleScoring检查一下。

5.4 一个典型的数据实战复盘

我用一个公开的PBMC数据集跑过一次完整的聚类流程,整个过程能帮助理解上面这些理论如何落地。

数据集大概有8000多个细胞。跑完PCA之后看ElbowPlot,前15个PC之后曲线基本上就平了,于是设dims = 1:15resolution = 0.5聚类,得到7个cluster。然后跑FindAllMarkers,根据已知marker很快注释出了CD14+单核细胞、FCGR3A+单核细胞、CD4 T细胞、CD8 T细胞、NK细胞、B细胞、树突状细胞,结果非常标准。

随后我尝试把resolution调到1.2,CD4 T细胞被拆分成了三个cluster:一个是初始T细胞(CCR7+、SELL+),一个是中央记忆T细胞(TCF7+、IL7R+),一个是效应记忆T细胞(GZMK+、IFNG+)。这种拆分在生物学上是合理的,因为CD4 T细胞的不同分化状态确实有着不同的表达特征。但如果我的研究问题只需要大类的细胞构成比例,拆分出来的细节反而是干扰。这就是为什么我一直强调,参数必须跟着问题走,没有一组参数是“万能最优”的。

6. 顺着这个思路还能做哪些扩展

理解了FindNeighborsFindClusters的原理之后,你实际上就掌握了一套通用的“网络聚类”思想,这个思想可以迁移到很多其他场景里。

比如空间转录组数据分析中,Seurat的FindSpatialClusters本质上也是基于类似的方法,只不过图的构建方式变成了“空间邻近”加“表达相似”的结合,KNN的邻居关系从细胞在PCA空间中的距离变成了空间坐标距离和表达谱距离的融合。有了细胞聚类的底子,再去理解空间域识别,上手会快很多。

又比如多组学整合时候的细胞类型对应问题。当同一个细胞类型在不同批次或不同平台的数据中都被聚出来了,你想判断这些cluster跨数据集是否对应,可以构建一个“簇级别的共现网络”,用网络聚类的思路去分群。这是我处理多批次数据时经常用到的技巧。

更直接的一个扩展是对聚类结果做更精细的注释。细胞类型注释的下一步是亚型识别,现在的流行思路是用FindAllMarkers后的top基因跟已知数据库比对,但更严谨的方式是先做基因集打分,比如用AddModuleScore计算一组已知转录因子打分,再结合聚类结果判断亚型的真实存在性。这类方法依赖的仍然是聚类结果的质量,所以把FindNeighborsFindClusters的原理吃透,是后续所有分析的地基。

7. 踩坑之后的几点体会

最近几次实操过程中,我有一个越来越强烈的感受:整个聚类流程,真正决定结果的往往不是FindClusters本身,而是它前面那个不起眼的dims参数,以及你有没有认真检查过降维和批次校正是否做对了。代码层面,FindNeighborsFindClusters合起来不过三四行,但每一行背后都有明确的数学逻辑和生物假设。

我自己吃了好几次亏之后,养成了一个固定的工作习惯:每次拿到新数据,先固定resolution梯度跑一遍clustree,同时把不同resolution下的umap排列出来,结合marker基因做一次统一的“人工审查”,这个流程看上去花时间,但往往能省下后面注释阶段反复返工的时间。另外一个经验是,如果最终注释出来的细胞类型跟文献里的预期相差很大,不要急着怀疑算法,先回头检查数据质量——很多分散的、无法解释的cluster,本质上是低质量细胞、双细胞或者批次效应在聚类图上的投影。

聚类只是单细胞分析这条漫漫长路的一个节点,但它决定了你的细胞注释能不能做准,也决定了后续所有统计检验的地基稳不稳。把这篇内容里的原理理解透了,你再来跑Seurat,应该会有一种“看得见算法在干什么”的感觉,而不是单纯地做一个代码执行者。这也是我在实际使用中最希望分享给每一位读者的东西。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询