4. 批次效应的来源与识别
说到批次效应，我刚开始做多批次质控的时候，也犯过糊涂。那时候我总觉得，只要实验流程严格遵守SOP，数据就应该是一样的。结果呢？第一次整合两个批次的芯片数据，PCA图上分成了两个清晰的“孤岛”。那个场景，我现在还记得很清楚。说白了，批次效应就像是数据的“胎记”，不是你流程做得好就能避免的。
4.1 批次效应到底从哪来？
你想想看，一个实验从样品采集到最终数据，中间经历了多少环节？每一个环节都有可能引入变异。我把常见的来源归纳成这么几类：
- 样本本身的变化：采集时间不同、保存条件不同、冻融次数不一样。我遇到过一个项目，第一批样本冻了两年，第二批是新采集的，结果两批次的基因表达谱差异巨大。
- 试剂与耗材的波动：同一家公司的试剂，不同批号之间也有差异。别低估这个影响，有时候就换了一个批号的裂解液，信号强度都变了。
- 操作人员的影响：这一点我深有体会。同一个protocol，不同人做，甚至同一个人不同时间做，结果都会有些差别。手法上那一点点的差异，放到高维数据里就放大了。
- 仪器设备的漂移：任何检测设备都有“老化”和“漂移”的问题。特别是荧光检测类的仪器，随着使用时间增加，检测灵敏度会有轻微变化。
- 环境因素：实验室温度、湿度，甚至季节变化，都可能造成微小的系统误差。
我个人的经验是：不要等到整合数据时才发现批次效应。最好在实验设计阶段就做好规划，比如把样本随机分配到不同批次，或者用一个“桥接样本”贯穿所有批次。亡羊补牢虽然不晚，但成本高得多。
| 来源类别 | 典型例子 | 影响程度 |
|---|---|---|
| 样本本身 | 采集时间、保存时长、冻融次数 | 高 |
| 试剂耗材 | 不同批号的试剂、芯片 | 中-高 |
| 操作人员 | 不同实验员操作、手法差异 | 中 |
| 仪器设备 | 不同扫描仪、检测器漂移 | 高 |
| 环境因素 | 温度、湿度、季节变化 | 低-中 |
嗯，这里要注意。上面这个表格只是一个大致的经验。具体哪个因素影响最大，得看你的实验体系是什么。比如做蛋白组学，试剂的批号差异可能就比做转录组大得多。
4.2 怎么判断是否存在批次效应？
讲完了来源，咱们聊聊怎么“抓”这个批次效应。别急着上复杂的统计模型，先用眼睛看。可视化是最直接的手段。
4.2.1 主成分分析（PCA）
这是我个人最推荐的第一眼看的方法。把所有样本做一个PCA，然后用颜色标记出不同的批次。
# R语言示例，做个简单的PCA
pca_result <- prcomp(t(expression_matrix), scale.=TRUE)
plot(pca_result$x[,1], pca_result$x[,2],
col=as.factor(batch_info),
pch=19, xlab="PC1", ylab="PC2")
如果不同批次的样本在PC1或PC2方向上明显“分群”，那就说明批次效应非常显著。我遇到过最夸张的情况，PC1直接就把两个批次分开了，连生物学差异都看不出来。
避坑指南：我曾经只盯着PCA图看，发现批次完全混在一起，就以为没有批次效应。结果后面做差异分析，还是发现了很多假阳性。为什么？因为批次效应可能不是“全局”的，它可能只影响一部分基因。所以PCA没问题，不代表数据完全干净。还得结合其他方法一起看。
4.2.2 箱线图和密度图
看全局特征。比如看所有样本的基因表达量中位数，如果不同批次之间差异很大，那就有问题。
# 用箱线图看样本表达量分布
boxplot(expression_matrix, outline=FALSE,
col=as.factor(batch_info)[sample_names],
main="Expression Distribution by Batch")
如果某个批次的箱体整体“高”一截或者“低”一截，那就是批次效应。我个人的习惯是，先用箱线图扫一眼，心中有数，再去做PCA。
4.2.3 相关性热图
样本间相关性热图也能看出端倪。如果批次内的样本相关性明显高于批次间的相关性，那就有批次效应了。
# 计算样本间相关性矩阵
cor_matrix <- cor(expression_matrix, method="spearman")
# 用热图可视化
heatmap(cor_matrix,
ColSideColors=batch_colors,
RowSideColors=batch_colors)
同一批次的样本应该能聚在一起，但如果批次效应太强，所有批次内的样本就自己抱团了，完全无视生物学分组。这种时候你就知道，该上Rank Normalization了。
4.2.4 统计检验
可视化之后，我还是会补一个统计检验。比如用ANOVA计算每个基因在批次间的差异，看p值的分布。
# 对每个基因做ANOVA，检验批次效应
p_values <- apply(expression_matrix, 1, function(gene) {
summary(aov(gene ~ batch_info))[[1]][1, "Pr(>F)"]
})
# 看p值分布
hist(p_values, breaks=50, main="P-value Distribution for Batch Effect")
如果p值直方图在0附近有一个尖峰，那就说明确实有很多基因受批次影响显著。
我上次去一个CDMO工厂做质控辅导，他们的QA团队就问我：“我做了这么多检验，到底信哪个？” 我的回答很简单：“先看PCA和箱线图，这是定性判断。再用统计检验定量，这是证据。两个都指向有批次效应，你就别心存侥幸了。”
这张图梳理了今天讲的核心逻辑。左边是批次效应的“罪魁祸首”，右边是我们用来“破案”的方法。中间那个虚线箭头意味着：你确认了批次效应的存在，接下来就需要评估它的严重程度，然后决定要不要用Rank Normalization来“摆平”它。