Finding stable clusterings of single-cell RNA-seq data

该论文提出了一种基于分裂式层次谱聚类及归一化割度量的方法，通过构建嵌套聚类树并利用互补子样本的一致性检验，实现了对单细胞 RNA-seq 数据聚类结果及其稳定性的评估。

要点

这篇论文就像是在教我们如何给单细胞 RNA 测序数据（scRNA-seq）里的细胞“分家”，并且要确保这个分法经得起考验，不会今天分明天就变。

想象一下，你手里有一大袋混合了各种水果（细胞）的果篮，你想把它们按种类（比如苹果、香蕉、橙子）分开。但问题是，这些水果长得都很像，而且有些可能还坏了（数据噪音）。怎么分才靠谱？

1. 核心难题：分得对不对，怎么验证？

通常，科学家分完类后，很难知道分得准不准。这就好比你把果篮分好了，但如果你多买一袋同样的水果，重新分一次，结果会不会完全不一样？如果每次分的结果都大相径庭，那说明你的分类方法太“随性”了，不可靠。

这篇论文提出的“稳定性测试”就像是一个“分身术”实验：

• 传统想法（很难实现）： 如果你能拿到两倍多的水果，重新分一次，看看结果变没变？但这在现实中很难做到（因为实验成本太高，或者样本已经没了）。
• 作者的新思路（反其道而行之）： 既然不能拿“更多”的数据，那我们就拿“更少”的数据试试。
  • 把原本的一大袋水果，随机切成两半（比如 A 组和 B 组）。
  • 分别对 A 组和 B 组进行“分家”。
  • 然后，看看 A 组里分出来的“苹果”，是不是和原本整袋水果里分出来的“苹果”是同一拨人？
  • 结论： 如果不管怎么切，分出来的结果都差不多，那说明这个分类是稳定的，是靠谱的。如果切一半就乱套了，说明这个分类太脆弱，不可信。

2. 具体怎么做？（分家三部曲）

作者设计了一套流程，就像是一个智能分拣流水线：

第一步：清洗与“去噪”（过滤坏水果）

在分家之前，先要把那些明显坏了或者长得太奇怪的水果挑出去。

• 剔除“欧几里得异常值”： 想象在果篮里，大部分苹果都挤在一起，突然有一个苹果离群索居，离其他苹果十万八千里。这种“独行侠”很可能是个坏苹果（噪音数据），或者是个特殊的变异体，先把它们挑出来，以免干扰大局。
• 剔除“基因异常值”： 有些基因在某个细胞里表达量特别高，高得不正常，就像某个水果突然发了疯一样长。这种基因也会把分类搞乱，所以也要剔除。

第二步：画地图与“分叉路”（聚类算法）

把剩下的好细胞放到一个多维度的“地图”上（欧几里得空间）。

• 分家策略： 作者没有用那种“一刀切”的方法，而是用了一种**“分裂树”**的方法。
  • 先把所有细胞分成两堆。
  • 再看这两堆里，哪一堆还能继续分？继续分。
  • 就像一棵大树，树干不断分叉，长出树枝，树枝再分叉，直到分不出为止。
• 关键指标（Normalized Cut）： 每次分叉，都要看这两堆分得干不干净。如果两堆之间还有很多“藕断丝连”的细胞，说明分得不好；如果分得很干脆，那就保留这个分法。

第三步：稳定性大考（反复测试）

这是最精彩的部分。作者把数据切成了 40 份不同的“小样”（就像把果篮随机切了 40 次）。

• 对每一份小样都跑一遍上面的“分家”流程。
• 打分： 看看小样里的分类结果，和原本大样里的结果有多像。
  • 如果 90% 的情况下，分出来的“苹果”还是“苹果”，那这个分类就是稳定的。
  • 如果分出来的“苹果”有一半变成了“香蕉”，那这个分类就是不稳定的，要扔掉。

3. 实验结果：分得怎么样？

作者用 7 个真实的数据集（就像 7 个不同的果篮）做了实验：

• 完美的案例（Zhengmix4eq）： 这是一个已知答案的果篮（里面只有 4 种细胞）。作者的方法分出来的结果和标准答案几乎一模一样，而且非常稳定。
• 有挑战的案例（68k PBMC 血液细胞）： 这是一个巨大的果篮，有 6 万多个细胞。
  • 作者发现，有些分法（比如分成 12 类）是稳定的，但其中有两类特别“调皮”，怎么分都分不匀（不稳定）。
  • 作者提出：只要那些“调皮”的类别细胞数量很少（少于 500 个），整个分法还是可以接受的。 这就像果篮里混进了几个烂果子，只要烂果子不多，整个果篮的分类还是靠谱的。
• 失败的案例（乳腺癌数据）： 这个果篮太复杂，或者坏果子太多。无论怎么分，结果都不稳定。这说明在这个数据里，可能根本不存在那么清晰的分类，或者数据质量本身有问题。

4. 总结与启示

这篇论文的核心思想可以用一句话概括：“真金不怕火炼，好分类不怕切分。”

• 不要盲目相信结果： 以前科学家分完类就完了，现在我们要问：如果我只用一半的数据，结果还一样吗？
• 接受不完美： 一个分类里，大部分细胞分得很准，只有极少数细胞分得乱，这没关系。只要把那些“乱分”的细胞剔除或标记出来，剩下的结果依然有价值。
• 科学的态度： 通过这种“反复切分测试”，我们可以把那些不可靠的、随机的分类结果过滤掉，只留下那些真正反映生物规律的稳定分类。

打个比方：
以前我们分水果，可能凭感觉分一次，觉得像苹果就是苹果。
现在，作者说：“别急，把这袋水果随机分成两半，再分一半，再分一半。如果每次分出来的‘苹果’都是同一批人，那才是真苹果；如果每次分出来的‘苹果’都不一样，那说明你根本分不清，或者这袋水果里根本没有苹果。”

这就是这篇论文想告诉我们的：在单细胞测序的世界里，稳定性（Stability）比什么都重要。

技术摘要

这是一份关于单细胞 RNA 测序（scRNA-seq）数据稳定聚类方法的详细技术总结，基于 Victor Klebanoff 的论文《Finding stable clusterings of single-cell RNA-seq data》。

1. 研究背景与问题 (Problem)

• 核心问题：在单细胞 RNA 测序（scRNA-seq）数据分析中，如何确定聚类结果的稳定性（Stability）和可重复性（Replicability）？目前对于基于 UMI（Unique Molecular Identifier）计数的数据，尚无共识性的方法来寻找稳定的聚类。
• 现有挑战：
  • 聚类结果往往对数据采样、参数设置（如分辨率）和异常值敏感。
  • 传统的聚类方法（如 Leiden 算法）在调整参数以获得不同数量的聚类时，往往缺乏系统性，且难以保证不同样本间结果的一致性。
  • 一个聚类方案可能包含非常稳定的簇，同时也包含极不稳定的簇，导致下游分析（如差异表达分析）产生误导。
• 核心假设：如果将数据量加倍（即增加一倍细胞），聚类结果是否会改变？由于无法获得“加倍”的真实数据，作者提出了反向视角：如果将现有数据随机分为两半，分别进行聚类，其结果是否与全量数据在对应子集上的聚类结果一致？如果这种一致性在多次采样中成立，则认为聚类是稳定的。

2. 方法论 (Methodology)

作者提出了一套完整的分析流程（Pipeline），主要包含以下步骤：

2.1 数据预处理与特征选择

• 过滤：保留在至少 50 个细胞中有非零计数的基因；对于部分数据集（如 PBMC），排除线粒体基因占比过高的细胞。
• 基因变异度计算：使用泊松模型计算每个基因的皮尔逊残差（Pearson Residuals）的平方和（SSQ）。
• 分析基因选择：仅保留在全集和所有子样本中均属于前 2000 个高变基因的基因。
• 降维：将皮尔逊残差矩阵视为低秩加噪声矩阵，使用 optht 程序（基于 Gavish & Donoho 算法）估计矩阵秩，并通过奇异值分解（SVD）将细胞映射到低维欧几里得空间。

2.2 异常值剔除 (Outlier Exclusion)

• 欧几里得异常值：基于 k-最近邻（kNN）距离分布。如果某细胞到其第 k 个最近邻的距离超过“均值 + 3 倍标准差”的阈值，则视为异常值并剔除。
• 迭代剔除：
  • 细胞异常值：基于细胞对基因 SSQ 的贡献率。如果某细胞在某个样本中对基因变异度的贡献过大（无法在互补样本中找到对应影响），则剔除。
  • 基因异常值：检查基因在不同样本间的 SSQ 一致性。如果某基因在不同样本间的 SSQ 比值过大，则视为异常基因剔除。
  • 该过程进行了三轮迭代，逐步剔除异常数据。

2.3 聚类算法：分裂式层次谱聚类 (Divisive Hierarchical Spectral Clustering)

• 替代方案：放弃了 Leiden 算法，因为其在调整分辨率参数以获取特定数量聚类时存在单调性问题和计算效率低的问题。
• 核心算法：采用 Ng, Jordan, and Weiss 的谱聚类算法，但进行了修改：
  • 亲和度定义：不使用高斯函数，而是定义为两点间欧几里得距离的倒数（$1/distance$）。
  • 稀疏性：仅对 k-最近邻（k=64）定义非零亲和度。
  • 分裂策略：递归地将点集分裂为两个簇，直到满足停止条件（最大树深度或最小簇大小）。
• 树映射：将生成的聚类树映射为嵌套的聚类集合。分支长度定义为“归一化割（Normalized Cut）”，距离根节点的距离决定了聚类的层级。

2.4 稳定性评估指标

为了评估聚类及其子簇的稳定性，作者设计了以下指标：

• 样本划分：将细胞随机分为互补的两组（$C_1, C_2$），重复 20 次（共 40 个样本）。
• 误分类误差距离 (MED)：比较子样本的聚类结果与全量数据在子样本上的限制聚类结果。MED 经过随机打乱标签后的均值进行归一化。
  • 稳定聚类标准：归一化 MED 的第 90 百分位数 $\le 0.10$。
• 簇误分类误差率 (CMER)：针对每个簇，计算子样本中被错误分类的细胞比例。
  • 稳定簇标准：归一化 CMER 的第 90 百分位数 $\le 0.50$（即至少 90% 的样本中，该簇内少于一半细胞被误分）。
  • 非常稳定/极度稳定：CMER 第 90 百分位数 $\le 0.10$ 或 $\le 0.02$。
• 可接受性：如果一个聚类方案包含不稳定簇，但每个不稳定簇的细胞数少于 500 个，则仍视为可接受用于下游分析。

3. 关键贡献 (Key Contributions)

1. 稳定性评估框架：提出了一种基于“子样本一致性”的严格框架，不仅评估整体聚类的稳定性，还量化了每个独立簇的稳定性（通过 CMER）。
2. 新颖的聚类映射方法：提出了一种将分裂式层次谱聚类树映射为嵌套聚类集合的方法，利用归一化割作为分支长度，能够系统地生成不同数量的聚类。
3. 异常值处理机制：设计了一套基于统计贡献（SSQ 贡献率）和几何距离（kNN 距离）的迭代异常值剔除流程，显著提高了聚类的鲁棒性。
4. 对现有方法的批判性验证：通过对比已知“真实标签（Ground Truth）”和已发表文献的结果，验证了该方法在不同数据集上的表现，并揭示了某些看似合理的聚类（如与文献一致）实际上可能是不稳定的。

4. 实验结果 (Results)

作者在 7 个公开数据集上进行了测试（3 个小数据集，4 个大数据集）：

• Zhengmix4eq (4 种细胞类型)：
  • 4-聚类结果与真实标签高度一致。
  • 所有 4 个簇均为极度稳定（CMER 极低）。
• Zhengmix8eq (8 种细胞类型)：
  • 7-聚类结果与真实标签吻合度较好（调整 Rand 指数 0.74），所有簇非常稳定。
  • 8-聚类（分裂了其中一个簇）稳定性下降，部分簇变得不稳定。
• CD14 Monocytes (单细胞类型)：
  • 所有聚类（2-10 簇）均不稳定（MED 中位数 > 0.50），表明数据本质上是同质的，任何细分都是人为的。
• 68k PBMC (外周血单个核细胞)：
  • 12-聚类（第一轮迭代）：被判定为可接受。虽然包含两个完全不稳定的小簇（<500 细胞），但其余 10 个簇非常稳定。
  • 9-聚类（第三轮迭代，与文献 k-means 结果更接近）：虽然与文献结果更相似，但包含完全不稳定的簇（Cluster 0, 4, 5），CMER 极高。这证明了“与文献一致”并不等同于“统计稳定”。
• 25k Retinal (视网膜)：
  • 11-聚类结果可接受。部分已发表的细胞簇（如杆状细胞）被稳定地分裂为两个簇，提示可能需要通过差异表达分析进一步验证这种分裂的生物学意义。
• 65k Lung (肺)：
  • 16-聚类表现出极高的稳定性（MED 第 90 百分位数仅为 0.01），所有簇均极度或非常稳定。尽管文献报道了 57 种细胞类型，但该方法认为 16 个稳定簇是数据内在结构的最佳反映。
  • 19-聚类虽然整体稳定，但包含两个完全不稳定的簇。
• 100k Breast Cancer (乳腺癌)：
  • 未发现稳定的聚类。所有尝试的聚类方案（2-51 簇）的 MED 第 90 百分位数均超过 0.10 的阈值。
  • 即使与文献报道的 9 种主要细胞类型最接近的 9-聚类，也包含极度不稳定的簇（如 Cluster 4）。
  • 异常值剔除过程对浆母细胞（Plasmablasts）造成了不成比例的剔除，暗示数据分布的高度异质性。

5. 意义与结论 (Significance & Conclusion)

• 重新定义下游分析标准：论文强调，只有稳定的聚类才应被用于下游分析（如差异表达）。不稳定的聚类可能反映了噪声或数据中的偶然结构，而非真实的生物学亚群。
• 揭示“虚假”一致性：研究发现，某些与已发表文献（或金标准）高度一致的聚类方案，在稳定性测试中可能表现极差。这提示研究者不能盲目依赖文献中的聚类结果，而应进行稳定性验证。
• 数据异质性的警示：在乳腺癌数据中，未能找到稳定聚类可能暗示该数据集存在极端的批次效应、细胞状态连续变化或异常值干扰，简单的离散聚类可能无法捕捉其真实结构。
• 方法论的普适性：该框架不依赖于特定的聚类算法（虽然使用了谱聚类），任何聚类流程都可以套用此稳定性评估方案。
• 未来方向：作者指出需要进一步研究如何优化异常值剔除的迭代次数，以及如何更好地利用差异表达分析来验证不稳定簇的生物学意义。

总结：这篇文章提出了一套严谨的、基于重采样的统计框架，用于评估 scRNA-seq 聚类结果的可靠性。它通过量化每个簇的稳定性，帮助研究者区分真实的生物学亚群和由噪声或采样偏差引起的虚假聚类，为单细胞数据分析提供了重要的质量控制工具。