批次效应识别与校正：从PCA到统计方法
做RNA-seq分析的人，迟早会撞上批次效应。我刚开始那会儿，拿到一个肝癌vs正常组织的表达矩阵，PCA一画，咦？正常样本分成了两坨，肝癌也分成了两坨。后来一查才知道，正常样本来自两个不同测序批次，癌样本也是。说白了，就是技术差异把生物学信号淹没了。
批次效应不是实验设计的锅，而是不可避免的“技术噪音”。不同试剂批次、不同操作人、不同测序仪、不同时间点，都会引入系统偏差。如果不校正，机器学习模型学到的可能全是批次差异，根本不是生物学差异。
这一节，我带你走一遍三个常用方法：PCA可视化看问题、ComBat参数校正、limma线性模型移除。每个方法我都踩过坑，会给你讲明白怎么用、何时用、别犯什么错。
1. PCA可视化：先看见问题
先别急着套校正工具。我习惯第一步是画PCA图，看看样本到底怎么聚类的。用R里自带prcomp()就能搞定，加上ggplot2染色批次信息。
# R代码示例：PCA可视化探索批次效应
pca_res <- prcomp(t(log2(counts + 1)), center = TRUE, scale. = TRUE)
pca_df <- as.data.frame(pca_res$x)
pca_df$batch <- sample_info$batch
pca_df$group <- sample_info$condition
library(ggplot2)
ggplot(pca_df, aes(x = PC1, y = PC2, color = batch, shape = group)) +
geom_point(size = 3) +
theme_minimal() +
labs(title = "PCA before batch correction")正常情况，同一种生物学条件（比如“肿瘤” vs “正常”）的样本应该先聚在一起。如果同一个条件内的样本被批次颜色分开，那就是批次效应。我见过一个项目，PCA前两个主成分解释方差不到15%，但批次却完美地把样本分成两群。这时候，不校正直接扔给模型，准确率看着高，但换个批次全翻车。
2. ComBat：经验贝叶斯校正
ComBat是我在多个项目中用得最多的方法，对付2~5个批次特别好使。它假设每个基因的批次效应是“可加的”和“乘性的”，用经验贝叶斯估计每个批次的调整参数。优点是很省心——你只需要给出表达矩阵和批次标签。
我用的是sva包的ComBat_seq函数，它专门针对counts数据（负二项分布家族）。如果你已经做了log2变换，那就用ComBat。
# 使用ComBat_seq校正counts数据
library(sva)
batch <- sample_info$batch
# 假设有2个批次：batch1, batch2
adjusted_counts <- ComBat_seq(
counts = counts,
batch = batch,
group = NULL, # 如果有生物学分组可传入，会保留组间差异
full_mod = FALSE
)
# 校正后PCA验证
pca_adj <- prcomp(t(log2(adjusted_counts + 1)), center = TRUE, scale. = TRUE)
pca_adj_df <- as.data.frame(pca_adj$x)
pca_adj_df$batch <- batch
ggplot(pca_adj_df, aes(x = PC1, y = PC2, color = batch)) + geom_point()这里有个坑：ComBat的group参数特别关键。如果你传入生物学分组（比如肿瘤vs正常），它会计算组内平均效应，不把真实生物差异当成批次去校正。我之前有一次忘了设group，结果把肿瘤和正常之间的差异也给抹平了，后续差异基因一个都没找到。嗯，这个错我犯过一次就记得死死的。
3. limma removeBatchEffect：线性模型取残差
limma包里有个removeBatchEffect函数，原理很简单：建一个线性模型，把批次效应当作协变量拟合，然后从表达量中减去批次效应的估计值。它不像ComBat那样做经验贝叶斯收缩，但计算快、灵活——你可以同时控制多个批次变量。
这个函数输入必须是log2变换后的表达矩阵（或者经过voom变换的）。用法如下：
# 使用limma的removeBatchEffect
library(limma)
# 假设有2个批次因子batch，以及生物学条件condition
design <- model.matrix(~ condition)
adjusted_expr <- removeBatchEffect(
log2_expr,
batch = batch,
design = design
)
# 然后就可以用校正后的矩阵做PCA或机器学习注意design参数是生物学分组的设计矩阵，它会保留你关心的差异。如果你不传design，函数默认所有样本属于同一组，那就会把生物学差异也一起“移除”——这基本是个灾难。我见过有人直接用removeBatchEffect(log2_expr, batch = batch)，以为批次校正会自动尊重组别，结果一画PCA，肿瘤和正常混在一起了。所以要记住：design一定要有你要保留的分组变量。
知识体系总览
下面这张图把批次效应识别的三步走概括起来：先可视化确认问题，再选择校正策略，最后验证校正效果。我用SVG画了核心流程，节点用的是我实际项目中最常用的工具和步骤。
校正后验证：别偷懒
无论用哪种方法，一定要再做一遍PCA（或者t-SNE/UMAP），看看批次是否被“混合”了。如果校正后批次仍然明显分开，说明方法不对或者参数没设好。我遇到过一种情况：ComBat校正后PCA上批次混得很均匀，但用KNN分类器预测批次标签，准确率还是高达95%——这说明批次效应并没有完全消除，只是在前两个主成分里看不出来了。所以，我建议用定量指标辅助判断：
| 指标 | 说明 | 理想值 |
|---|---|---|
| PCA批次方差贡献 | 批次主成分解释方差比例 | < 5% |
| Silhouette得分（基于批次） | 样本与同批次相似度 vs 不同批次 | 接近0 |
| KNN批次分类准确率 | 用基因表达预测样本批次 | 接近随机（1/k） |
实际项目中，我通常会跑完校正后算一个“批次分类准确率”，如果低于50%（对于2个批次），我就觉得可接受。当然，如果它仍然很高，可能就需要回归到数据本身，看看是否有个别批次存在严重质量问题——比如测序深度异常、rRNA残留等，这些不是一个统计校正能解决的。