批次效应检测方法：PCA可视化、层次聚类、相关系数矩阵、批次效应评分
检测批次效应，说白了就是回答一个问题：样本之间的差异到底是生物学原因造成的，还是实验批次带来的假象？我个人习惯把这四种方法分成两类来看。一类是视觉诊断——PCA、层次聚类、相关系数矩阵；另一类是定量打分——批次效应评分。两类互相补充，缺一不可。
你想想看，如果只靠肉眼看图，不同批次之间颜色分开了，你心里有数；但要是分不开呢？不一定代表没批次效应——有可能是维度压缩把信号丢掉了。反过来，只算一个评分数字，0.02还是0.15？没有视觉印证，你很难说服审稿人。所以我的建议是：先用视觉方法快速扫描，再用评分方法定量确认。
下面我一个一个拆开讲，每个方法都带代码和避坑经验。
PCA 可视化：最常用的第一道防线
PCA 的原理不用我多啰嗦，就是把高维表达矩阵压缩到前两个或三个主成分上。批次效应要是存在，不同批次的样本会在 PC 空间里各自抱团。
我在项目中遇到过最典型的情况：PC1 把处理组和对照组分开了，PC2 却把批次1和批次2分开了。这时候心里就有数了——批次效应确实存在，而且方向正交于生物学变量，后续校准大概率能救回来。
下面是一段标准的 Scanpy 代码：
import scanpy as sc
adata = sc.read('merged_batches.h5ad')
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000)
adata = adata[:, adata.var.highly_variable]
sc.pp.scale(adata, max_value=10)
sc.tl.pca(adata, n_comps=30, svd_solver='arpack')
sc.pl.pca(adata, color='batch',
title='PCA colored by batch',
show=True)解读 PCA 图的时候，我一般看三点：
- 批次样本是否形成独立团簇？如果是，批次效应基本实锤了。
- 团簇之间的重叠程度如何？部分重叠说明效应较弱，后续校准可行。
- PC1/PC2 的方差解释比例。如果前两维只解释了不到 10% 的方差，说明数据噪声大，PCA 可能不是最佳检测工具。
嗯，这里还要提一个坑。我曾经遇到过一个数据集，PCA 图上批次完全混合，看起来干净得很。但我换了 t-SNE 和 UMAP 一看，批次效应在局部结构里藏得很深。所以 PCA 没发现问题，不代表真的没问题——这也是为什么我们需要下面几种方法。
层次聚类：无监督视角的交叉验证
层次聚类的好处是不做降维，直接用样本间的距离矩阵来构建树状图。批次效应如果存在，同一批次的样本会在树的分支上聚到一起。
我曾经拿一个四批次的数据集做过测试——PCA 图上前三个批次混在一起，第四个批次稍微偏了一点。结果层次聚类树一画出来，四个批次在深层分支上清清楚楚分开了。这就是层次聚类的价值：它能看到 PCA 压缩掉的那部分结构信息。
代码示例：
import matplotlib.pyplot as plt
from scipy.cluster.hierarchy import dendrogram, linkage
from scipy.spatial.distance import pdist
# data 是 样本 × 基因 的矩阵，已经标准化
dist_matrix = pdist(data, metric='euclidean')
Z = linkage(dist_matrix, method='ward')
plt.figure(figsize=(12, 6))
dendrogram(Z, labels=batch_labels, leaf_rotation=90,
leaf_font_size=8)
plt.title('Hierarchical clustering of samples')
plt.ylabel('Distance')
plt.show()ward 还是 complete？我一般先用 ward，它对团簇大小的敏感性较低。如果树状图看起来太"平"，没有清晰的分支，再换 complete 试试。两种结果对比一下，心里更有底。
解读树状图时，我关注两个指标：
- 同一批次的样本是否在树中形成单系群？如果是，说明批次效应在全局距离结构里占主导。
- 不同批次的样本在树中是否交替出现？交替出现是好信号，说明生物学变异超过了批次变异。
相关系数矩阵：批量评估样本间的一致性
相关系数矩阵，说白了就是算任意两个样本之间的表达谱相似度。然后画成热图，批次效应会在热图上形成"方块"结构。
我记得有一次帮同事做技术重复的质量评估。两个人分别建库、上机，测同一个细胞系。我算了样本间的 Spearman 相关系数，技术重复之间的相关性在 0.98 以上，但不同技术重复之间的相关性只有 0.85——这就是典型的批次效应，而且来源很可能是建库操作上的差异。
代码：
import seaborn as sns
import pandas as pd
# data: 样本 × 基因 的标准化矩阵
corr_matrix = data.T.corr(method='spearman')
plt.figure(figsize=(10, 8))
sns.heatmap(corr_matrix,
xticklabels=batch_labels,
yticklabels=batch_labels,
cmap='RdBu_r', center=0.8,
square=True, cbar_kws={'shrink': 0.8})
plt.title('Sample correlation matrix')
plt.show()批次效应评分：用数字说话
前面三种方法都是视觉诊断，但审稿人更喜欢数字。批次效应评分就是干这个的。常用的几个评分我把它们整理成一个表格：
| 评分方法 | 核心思想 | 数值含义 | 我的建议 |
|---|---|---|---|
| kBET | 检验每个样本的 k 近邻中批次标签分布是否随机 | 0~1，接近 0 表示无批次效应，接近 1 表示强批次效应 | 单细胞数据首选，对局部结构敏感 |
| LISI | 计算局部逆辛普森指数，衡量批次混合程度 | 值越高表示混合越好，批次效应越弱 | 适合与 kBET 配合使用，一正一反互相验证 |
| ASW | 基于轮廓宽度的批次效应版本 | 0~1，接近 0 表示批次完全混合，接近 1 表示批次完全分离 | 对全局结构敏感，适合快速评估 |
| PCR | 用主成分回归，看批次标签能解释多少方差 | 百分比形式，越高表示批次效应越强 | 适合大样本、高维度数据，计算快 |
以 kBET 为例，代码实现如下：
import kbet
# adata.obs['batch'] 是批次标签
# adata.X 是标准化后的表达矩阵
result = kbet.kbet(adata.X,
adata.obs['batch'],
k0=10) # 近邻数，默认 sqrt(n_samples)
# result 包含每个样本的 kBET 拒绝率
acceptance_rate = result['acceptance_rate']
print(f'kBET acceptance rate: {acceptance_rate:.3f}')
# 如果 acceptance_rate < 0.5，表明存在明显的批次效应四种评分方法各有侧重点，我个人习惯至少算两个——一个全局的（ASW 或 PCR）和一个局部的（kBET 或 LISI）。两者结论一致，基本可以下判断了。
四种方法的关系与选择
讲了这么多，你可能要问：到底用哪个？我的回答是：不是二选一，而是组合使用。下面是这四种方法在整个检测流程中的位置关系：
这张图把检测流程串起来了。底色用浅蓝是为了让大家看清节点之间的关系——从原始数据出发，经过标准化，然后四个检测方法同时跑，最后汇总到综合评估里，得出最终结论。流程上的箭头都是双向验证的意思，不是单向的流水线。
总结一下我个人的经验：
- 初步筛查：先跑 PCA 和层次聚类，两个图一起看。PCA 给全局视角，层次聚类给距离结构。
- 定量确认：算相关系数矩阵，看批次间是否有系统性偏差。再算一个 kBET 或 ASW 评分，拿到数字证据。
- 交叉验证：如果 PCA 和层次聚类结论一致，评分也给出同方向信号，那批次效应基本实锤了。
- 不一致怎么办？我遇到过 PCA 显示没问题、但 kBET 评分很低的案例。这时候我会去看数据是不是存在"局部批次效应"——某些基因子集受批次影响大，但 PCA 被其他基因稀释了。处理方法是对高变异基因单独做一次检测。