4. 常用的批次效应检测方法介绍
说实话，批次效应检测这件事，我踩过的坑比学过的理论多。早些年做多组学整合，愣是把生物学差异和批次效应混在一起，结果结论全反了。后来我养成了一个习惯：校正之前，先花十分钟做检测。这十分钟能省下后面好几天返工的时间。
下面这几种方法，是我在项目中反复验证过的。它们各自有适用场景，我建议你按顺序组合使用——先全局鸟瞰，再局部排查。
4.1 主成分分析（PCA）—— 最直观的“第一眼”
PCA 是我个人最常用的探索性工具。把高维基因表达数据压缩到前几个主成分，然后按批次给样本点着色。如果同一批次的样本在 PC 空间里扎堆，而不同批次的样本明显分离，那就有问题了。
操作起来很简单：
# R 示例：pca 后按批次着色
pca_result <- prcomp(t(expr_matrix), scale. = TRUE)
plot(pca_result$x[,1], pca_result$x[,2],
col = factor(batch_vector), pch = 19,
xlab = "PC1", ylab = "PC2")
legend("topright", legend = levels(factor(batch_vector)),
col = 1:length(unique(batch_vector)), pch = 19)你想想看，如果前两个主成分的方差解释率超过 30%，并且批次分离明显，基本可以判定存在批次效应。但要注意：PCA 太灵敏了，有时候生物学差异也会导致分离，所以需要结合后续方法交叉验证。
个人经验：我习惯在 PCA 的基础上再加一个“批次主成分贡献”图——计算每个主成分与批次标签的相关性。如果某个 PC 与批次的相关性 >0.7，那个 PC 基本就是被批次效应主导的。这个方法比肉眼判断更客观。
4.2 层次聚类与相关性热图
层次聚类是 PCA 的补充。我不看树形图的具体分支，而是看样本是否按批次聚在一起。
# 计算样本间相关系数
cor_mat <- cor(expr_matrix, method = "spearman")
# 层次聚类
hc <- hclust(as.dist(1 - cor_mat), method = "ward.D2")
# 绘制带批次标签的聚类树
plot(hc, labels = batch_vector, cex = 0.6, main = "Sample Clustering")我记得有一次，聚类树上所有批次 1 的样本先聚成一簇，批次 2 的样本再聚成另一簇——中间几乎没有交叉。嗯，那就是典型的批次效应，而且强度很大。
避坑指南：我曾经在50个样本上直接画全相关性热图，结果颜色块密密麻麻根本看不出来。后来我改用批次内平均相关系数 vs 批次间平均相关系数的对比。如果批次内相关系数比批次间高出一大截（比如差值 >0.1），基本可以确认有批次效应。
4.3 UMAP / t-SNE —— 别忘了设置随机种子
t-SNE 和 UMAP 是近年流行的可视化方法。它们能保留数据的局部结构，往往比 PCA 更清晰地呈现批次分离。但有一个大坑：随机性。
# R 示例：UMAP 可视化批次
library(umap)
umap_config <- umap.defaults
umap_config$random_state <- 42 # 固定随机种子！
umap_result <- umap(t(expr_matrix), config = umap_config)
plot(umap_result$layout, col = factor(batch_vector), pch = 19,
xlab = "UMAP1", ylab = "UMAP2")为什么一定要固定种子？因为 UMAP 每次运行结果都可能不同。你不固定种子，今天看到的分离，明天可能就消失了。我团队里有个新人没设种子，汇报的时候领导让他重现结果，跑了三次三次都不一样——那场面相当尴尬。
另外，t-SNE 对困惑度参数敏感。我曾遇到过困惑度设太小，局部结构被放大，把本来没关系的一批样本强行拉在一起。我建议先用困惑度 = 30 试跑，如果不行再调。
4.4 方差分析（ANOVA）—— 定量检测的硬指标
可视化方法多少有点主观。我们还需要一个数字指标来量化批次效应的大小。方差分析（ANOVA）正好做这个。
对于每个基因，我们计算“批次”因素解释的方差比例。如果很多基因的批次方差占比超过 30%，那批次效应就很明显了。
# 对每个基因做 ANOVA，提取批次效应占比
pve_batch <- apply(expr_matrix, 1, function(gene) {
df <- data.frame(expression = gene, batch = batch_vector)
fit <- aov(expression ~ batch, data = df)
ss <- summary(fit)[[1]]$`Sum Sq`
ss[1] / sum(ss)
})
hist(pve_batch, breaks = 50,
main = "Distribution of Proportion Variance Explained by Batch",
xlab = "PVE_batch", col = "steelblue")
abline(v = 0.3, col = "red", lwd = 2, lty = 2)红色虚线标记 0.3。如果大多数基因的 PVE 都低于这条线，说明批次效应可控；反之就要警惕了。我个人觉得，当超过 20% 的基因 PVE > 0.3 时，必须做校正。
核心指标速查表：
| 指标 | 判断阈值 | 备注 |
|---|---|---|
| PCA 批次分离 | PC1/PC2 按批次明显分离 | 需结合生物学分组验证 |
| 批次内 vs 批次间相关系数差 | 差值 > 0.1 | 适合大规模芯片数据 |
| UMAP 批次聚集 | 同批次样本在局部区域密集 | 务必固定随机种子 |
| PVE_batch > 0.3 的基因占比 | > 20% | 定量硬指标 |
4.5 一种更系统的方法：Silhouette 宽度评估
最后补充一个我最近特别喜欢用的方法——轮廓系数（Silhouette）。它不仅能评估聚类效果，还可以用来衡量批次分离程度。
计算每个样本基于批次标签的轮廓宽度。如果整体平均轮廓宽度 > 0.5，说明批次之间的分离非常明显。如果 < 0.1，说明批次混杂得比较好。
# 基于批次标签的轮廓宽度
library(cluster)
dist_mat <- dist(t(expr_matrix))
sil_batch <- silhouette(as.numeric(factor(batch_vector)), dist_mat)
mean(sil_batch[,3]) # 输出平均轮廓宽度我在一个 8 批次的数据集上跑过，平均轮廓宽度 0.62。这个数字让我果断决定要用 ComBat 校正。校正后再跑一次，降到了 0.12——效果立竿见影。
4.6 知识结构总览
上面这些方法不是孤立的。我习惯按一个固定的流程跑：先 PCA 看全貌，再 UMAP 看局部，然后用 ANOVA 定量，最后用 Silhouette 做最终确认。下面这张图把流程串了起来。
我的个人工作流：拿到新数据后，我先 PCA + 层次聚类（5分钟），如果看出明显批次分离，立刻跑 ANOVA 计算 PVE 分布。然后 UMAP 做双重确认。最后用 Silhouette 给一个数字结论。整个流程跑下来不超过半小时，但能帮我避免 80% 以上的批次校正失误。
嗯，工具和方法就这些。下一件事，就是根据检测结果决定要不要做 ComBat 校正了。不过别急，检测这一步做扎实了，后面的校正才有意义。