第1章 批次效应概述：从一次翻车经历说起
说起批次效应，我就想起几年前的一个项目。当时我在做一个多中心的肿瘤标志物筛选，三个中心的血清样本，同一个实验方案，连试剂盒都是同一批采购的。结果PCA一画，三个中心的数据各自抱团，压根不按疾病分组来——我盯着那张图愣了五分钟。嗯，那就是批次效应的典型翻车现场。
从那以后，我养成了一个习惯：拿到任何组学数据，第一件事不是跑差异分析，而是先查批次效应。今天咱们就把这个"数据隐形杀手"彻底聊透。
1.1 什么是批次效应
批次效应（Batch Effect），说白了就是——不同批次处理的样本之间，那些跟生物学无关、纯粹由实验流程引入的系统性差异。
它跟"实验误差"不太一样。随机误差是上下乱跳，今天大明天小，你多做几次重复就能平均掉。批次效应是"定向偏移"：这一批所有样本的表达量整体偏高了，下一批整体偏低了。它不是噪声，是系统偏差。
举个直观的例子：
- 你在周一跑了一板芯片，温度22℃，试剂新鲜，操作顺手；
- 周三又跑了一板，温度升到28℃，试剂放了三天，换了个实习生加样；
- 结果两板数据放在一起，周一的样本整体表达量比周三高了30%。
这个30%的差异跟生物学完全没关系，它就是批次效应。
1.2 批次效应的产生原因
我在项目中统计过，批次效应的来源其实就那几类。下面这张表总结了最常见的"罪魁祸首"：
| 原因类别 | 具体表现 | 常见场景 |
|---|---|---|
| 实验条件波动 | 温度、湿度、pH值、孵育时间 | 芯片杂交、PCR扩增、细胞培养 |
| 操作者更换 | 加样手法、判读标准、操作习惯 | 多中心实验、长期队列 |
| 试剂批次变更 | 抗体效价、酶活性、试剂浓度 | ELISA、Western Blot、RNA-seq |
| 仪器设备漂移 | 激光强度、检测器灵敏度、管路残留 | 流式细胞仪、质谱仪、测序仪 |
| 时间跨度 | 样本采集日期、冻存时长、反复冻融 | 纵向队列研究、回溯性分析 |
1.3 批次效应对数据分析的影响
批次效应不是"有点讨厌"的问题，它是能让你整个分析报废的坑。我总结了三类最严重的后果：
- 假阳性激增 —— 你筛出来的"差异基因"，其实只是批次间的差异。两个批次的样本如果碰巧对应了不同的实验组，那几乎所有的差异都是批次效应，跟生物学无关。
- 真实信号被掩盖 —— 生物学差异可能只有20%，而批次效应有50%。后者直接把前者淹没了，你什么都看不到。
- 模型泛化能力归零 —— 用带批次效应的数据训练分类器，模型学到的其实是"区分批次"而不是"区分疾病"。换一批新样本，准确率直接跳水。
1.4 如何识别批次效应
别等到分析做完了才发现有问题。我自己的流程是：拿到数据后，三步走，十分钟就能判断批次效应严不严重。
第一步：PCA可视化
用主成分分析降维到二维，按批次着色。如果不同批次的样本各自抱团、明显分离，那就敲响警钟了。
# R语言示例：PCA快速检查批次效应
library(ggplot2)
library( factoextra)
# 读取表达矩阵（行=基因，列=样本）
expr_matrix <- read.csv("expression_data.csv", row.names = 1)
# 读取批次信息
meta_data <- read.csv("meta_data.csv")
# 转置后做PCA
pca_obj <- prcomp(t(expr_matrix), scale. = TRUE, center = TRUE)
# 可视化，按批次着色
p <- fviz_pca_ind(pca_obj,
geom.ind = "point",
col.ind = meta_data$Batch,
palette = c("#3498db", "#e74c3c", "#2ecc71"),
addEllipses = TRUE,
ellipse.level = 0.95)
print(p)第二步：层次聚类
计算样本间距离，做聚类树，用批次信息标记叶节点颜色。如果同一批次的样本在树中形成独立分支，那批次效应基本实锤了。
# 计算样本间距离
dist_mat <- dist(t(expr_matrix))
# 层次聚类
hc <- hclust(dist_mat, method = "ward.D2")
# 可视化，按批次着色
plot(hc, labels = meta_data$Batch, main = "样本聚类树（按批次标记）")第三步：统计检验
对每个基因做ANOVA，检验"批次"因素的显著性。统计显著的比例如果远高于预期的5%，说明批次效应是系统性的。
# 对每个基因做批次效应检验
p_values <- apply(expr_matrix, 1, function(gene) {
df <- data.frame(expression = gene, batch = factor(meta_data$Batch))
summary(aov(expression ~ batch, data = df))[[1]][["Pr(>F)"]][1]
})
# 看看多少基因的批次效应显著（p < 0.05）
cat("批次效应显著的基因比例:", mean(p_values < 0.05, na.rm = TRUE))上面三步做完，你对数据里的批次效应就有了底。清不清晰？严不严重？心里有数了，下一步才能决定是用ComBat还是limma做校正。