3. 批次效应评估指标：ARI、NMI、批次熵、LISI 等定量指标

你辛辛苦苦跑完校正，t-SNE 图上一团五颜六色的细胞混在一起，看起来很美。但光靠眼睛看能算数吗？说实话，可视化只能给你一个初步印象，甚至有误导性——我早年就被漂亮的 UMAP 图骗过，后来才发现细胞类型标签根本没对齐。

所以，咱们需要数字。一套能统一量化的指标，让你知道校正到底好不好。

这一节我带你过一遍最常用的几个评估指标：ARI、NMI、批次熵和 LISI。每个我都会讲清楚它衡量什么，在代码里怎么算，以及我亲自踩过的坑。

核心原则：没有任何一个指标是万能的。好的评估策略是「聚类保真 + 批次混合」两个维度各选一个指标组合使用。

3.1 为什么需要定量指标？

我们先理清一个底层问题：可视化（t-SNE / UMAP）到底有什么缺陷？

  • 主观性强：不同的人对「混得好不好」判断不一样。
  • 参数敏感：perplexity、n_neighbors 一变，图就变了。
  • 保距性差：t-SNE 和 UMAP 都不保持全局距离，有些结构是假象。

所以，你必须用数字来回答两个问题：

  1. 生物异质性保留了吗？——同一细胞类型是否仍聚在一起。
  2. 批次效应消除了吗？——不同批次的同类细胞是否混在一起。

下面的 SVG 把评估任务和对应的指标体系映射出来了。你一眼就能看到哪个指标解决哪个问题。

批次效应评估指标体系 评估任务：批次效应校正效果 聚类保真（Cluster Fidelity） 类别混淆度（Label Confusion） 批次混杂度（Batch Mixing） 局部混合度（Local Mixing） ARI | NMI ARI 批次熵 | kBET iLISI | cLISI 记住：没有单一指标万能 组合使用「聚类保真 + 批次混合」两个维度的指标才是正确的评估方式

3.2 ARI — 看聚类标签跟真实标签对上了多少

ARI（Adjusted Rand Index）你之前大概率见过。它衡量两个聚类结果之间的相似度，跟 Rand Index 比多了一个「随机校正」。

说白了，ARI 计算的是：
把你的聚类分组 和 真实细胞类型标签 做对比，看它们有多吻合。
值域：[-1, 1]。1 表示完全吻合，0 表示随机水平，负数表示还不如随机。

公式简单提一下（别看晕）：

ARI = (RI - Expected_RI) / (max(RI) - Expected_RI)

其中 RI 是 Rand Index，Expected_RI 是随机排列下的期望 RI。

嗯，这里要注意：ARI 对标签顺序完全不敏感。你把类别 0 和 1 互换，结果不变。这是它比简单准确率强的地方。

个人习惯：在单细胞整合场景下，我通常把「校正后的聚类标签同真实细胞类型」做 ARI 对比。校正做得越好，ARI 越接近 1。不过如果校正过度（把不同细胞类型强行拉在一起），ARI 反而会掉——这也是个监警机制。

用 Python scikit-learn 算 ARI 超级简单，两行代码搞定：

from sklearn.metrics import adjusted_rand_score

# true_labels: 真实细胞类型标签 (list or array)
# pred_labels: 校正后聚类得到的标签 (list or array)
ari = adjusted_rand_score(true_labels, pred_labels)
print(f"ARI = {ari:.4f}")

你可能会问，那直接用 ARI 评价整合效果不就完了？
我建议你一定要配合同步监控「批次混合度」指标。因为：

  • 只优化 ARI → 可能聚类很准，但批次根本没混。
  • ARI 过高但批次效应还存在 → 大概率是聚类粒度太粗。—— 我见过有人用 k=2 的 KMeans 跑出 ARI=0.9，毫无意义。

3.3 NMI — 归一化互信息

NMI（Normalized Mutual Information）是另一个常用的聚类对比指标。

它跟 ARI 的区别在于：ARI 基于「对与错的计数」，NMI 基于信息论里的互信息。

为什么需要它？
简单说，ARI 偏爱平衡的簇大小。如果你的细胞类型分布严重不均（比如巨噬细胞占 80%，其余几种占 20%），ARI 会受少数类影响很大。NMI 在这方面稍微稳健一些。

计算公式取决于你用什么归一化方式（算术平均、几何平均、min 等），但结果都是 [0, 1]——1 表示完全一致，0 表示无共享信息。

from sklearn.metrics import normalized_mutual_info_score

nmi = normalized_mutual_info_score(true_labels, pred_labels)
print(f"NMI = {nmi:.4f}")
坑来了：scikit-learn 的 normalized_mutual_info_score 默认使用 average_method='arithmetic'。在某些旧版 scikit-learn 里不是这个默认值，建议显式指定。另外，当 true_labels 和 pred_labels 的簇数差异很大时，NMI 会偏高。我曾经吃过这个亏，以为 NMI=0.95 整合完美，其实只是聚类粗。

3.4 批次熵 — 看我能不能猜出你来自哪个批次

上面 ARI 和 NMI 都是在看「生物信号保留得怎么样」。那「批次混杂度」呢？

这里有一个直观指标：批次熵。

核心思想：如果校正做得好，你从任何一个小区域（比如一个细胞邻域）里随机抽一个细胞出来，应该很难猜出它来自哪个批次。换句话说，批次的分布应该是均匀的。

两个常用版本：

  • Shannon entropy 版本：在每个邻域里，看批次标签的熵有多接近理论最大值。熵越高代表混合越均匀。
  • 线性模型 R² 版本：把批次作为分类目标，训练一个分类器，看 R² 多大。R² 接近 0 则表示批次猜测不准确 → 混合好。

我用得比较多的是第二种，因为快且直观：

import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import r2_score

# X: 校正后的 embedding (n_cells × n_features)
# batch: 批次标签 (0, 1, 2, ...)
clf = LogisticRegression(multi_class='multinomial', max_iter=1000)
clf.fit(X, batch)
batch_pred = clf.predict(X)

# 注意：这里 y_true 和 y_pred 都是离散类别，用 R2 的 "score" 其实不对，
# 更严谨的是用伪 R² 或准确率。
# 但为了展示概念，我们直接用准确率作为 R² 的近似（严格讲不准确）。
# 实际项目里你可以计算分类准确率，接近 25% (4个批次) 说明批次被洗掉了。
from sklearn.metrics import accuracy_score
acc = accuracy_score(batch, batch_pred)
print(f"批次预测准确率 = {acc:.4f}")

注意：使用准确率评估批次混合度时，不要把批次标签特别稀疏的品种也算进来。如果某批次只有十几细胞，分类器很容易猜对/错，结果是虚假的。碰到这种情况，我会先均衡采样再做预测。

还有一个更标准的方法是 kBET（k-Nearest Neighbor Batch Effect Test）。它通过卡方检验来评估每个细胞局部邻域内的批次分布是否均匀。R 包 kBET，Python 有 scIB 包包装了它。不过我坦白说，kBET 计算量有点大，大数据集下我很少用，优先看 LISI。

3.5 iLISI / cLISI — 局部结构混合度

LISI（Local Inverse Simpson's Index）是我自己项目中最常用的指标。它有两个变体：

  • iLISI（Integration LISI）：衡量不同批次的细胞在局部邻域里的混合程度。值越大混得越好。
  • cLISI（Cell-type LISI）：衡量同一细胞类型在局部邻域里的纯度。值越小纯度越高（想一下：如果同一细胞类型的细胞在邻域里占比高，LISI 值会低）。

iLISI 越接近「批次数」越好；cLISI 越接近 1 越好（表示完全纯的细胞类型邻域）。

Python 里用 milopy 这个库算非常方便：

import milopy
import milopy.core as milo

# adata: AnnData 对象，包含校正后的 embedding 和批次、细胞类型信息
milopy.utils.compute_lisi(
    adata,
    batch_key='batch',
    label_key='cell_type',
    use_rep='X_pca'    # 或者 'X_umap'，看具体情况
)

# iLISI 就在 adata.obs['iLISI'] 里
# cLISI 在 adata.obs['cLISI'] 里
print(adata.obs[['iLISI', 'cLISI']].describe())
踩坑经验：iLISI 对 perplexity（或用多少邻居）敏感。我习惯用两个值都试一遍：perplexity=30 和 perplexity=100。如果两个结果矛盾，说明你的邻域结构在不同尺度下表现不一致，这时候优先看中等尺度（perplexity~50 或 30）。另外，千万不要在 UMAP 坐标上算 LISI，UMAP 会扭曲局部结构。一定要用校正后的 PCA 或 Harmony 输出空间。

3.6 怎么选指标？一张表讲清楚

最后我把几个指标放在一起对比，方便你选择：

指标 测量维度 度量对象 取值区间 对标签的要求 我推荐的使用场景
ARI 聚类保真 聚类 vs 真实标签 [-1, 1] 需要真实类别 常规整合效果评估，适合平衡数据集
NMI 聚类保真 聚类 vs 真实标签 [0, 1] 需要真实类别 类别分布极端不平衡时
批次熵 / kBET 批次混合 邻域内批次分布 视方法而定 仅需批次标签 快速检查批次是否被洗掉
iLISI 批次混合 邻域内批次多样性 [1, n_batch] 仅需批次标签 细粒度局部混合评估，我最常用
cLISI 类别纯度 邻域内细胞类型纯度 [1, n_celltypes] 需要细胞类型 监控是否过度校正

3.7 实际工作流

说了这么多指标，具体怎么用？我总结成两步：

  1. 聚类保真：计算 ARI 或 NMI（选一个，不要两个都用，否则会过拟合同一个方向）。如果类别不平衡，优先 NMI。
  2. 批次混合：计算 iLISI 和 cLISI。iLISI 越高越好，cLISI 越低越好。同时也可以看批次熵或 kBET 作为补充。

一个典型的工作流代码长这样：

from sklearn.metrics import adjusted_rand_score, normalized_mutual_info_score
import milopy
import milopy.core as milo

# 1. 聚类保真
ari = adjusted_rand_score(adata.obs['cell_type'], adata.obs['cluster'])
nmi = normalized_mutual_info_score(adata.obs['cell_type'], adata.obs['cluster'])

# 2. 批次混合 — LISI
milopy.utils.compute_lisi(adata, batch_key='batch', label_key='cell_type', use_rep='X_pca')
iLISI_median = adata.obs['iLISI'].median()
cLISI_median = adata.obs['cLISI'].median()

print(f"ARI={ari:.4f} | NMI={nmi:.4f} | iLISI_median={iLISI_median:.2f} | cLISI_median={cLISI_median:.2f}")

看到输出结果后，你就能回答三个问题：

  • ARI / NMI 高吗？—— 生物信号保留得好吗？
  • iLISI 接近批次数吗？—— 批次混得均匀吗？
  • cLISI 接近 1 吗？—— 过度校正出现没？

我个人的小建议：不要只看中位数或均值。记得绘制指标的分布图。iLISI 的直方图如果出现一个长尾（很多细胞 iLISI 很低），说明即使中位数还行，有些细胞类型的批次效应仍然很严重。这种细节会被中位数掩盖。

3.8 别忘了「稳定性」本身也是一个检验维度

最后跟各位分享一个容易被忽略的点。

你评估校正结果时，建议跑 多个随机种子。因为很多整合方法（比如 scVI、Harmony 的初始化）有随机性。同一个方法跑三次，如果 ARI 波动超过 0.05，说明算法不稳定。这时候你再怎么挑指标也没用。

我曾经在一个项目里用 MNN Correct，发现同一种参数在两次运行之间 ARI 波动达到 0.12。后来一看，原来是邻居选择阶段的随机种子没固定。所以我的习惯是：
固定种子跑一遍主要评估，然后再跑三次不同种子看方差。

← 上一章 📖 返回目录 下一章 →