第3章 批次效应的检测方法
拿到数据后，第一步不是急着校正，而是——先看看有没有批次的鬼影子。我见过太多人上来就一顿校正操作，结果把生物学差异也给抹掉了，得不偿失。检测批次效应，说白了就是回答两个问题：“批次之间真的不一样吗？” 和 “这个不一样到底有多大？” 这一章咱们就聊聊怎么用可视化和统计检验来回答这两个问题。
3.1 可视化方法：先“看”再“算”
我个人习惯，数据拿到手先画个PCA。为什么？因为PCA没有假设、不参杂主观偏见，它只是忠实地把高维数据投影到几个主成分上。如果你发现不同批次的样本在主成分图上各自抱团、互不搭理——嗯，那就是批次效应在跟你打招呼。
3.1.1 PCA：快速定位全局偏移
以RNA-seq数据为例，我一般用前5000个高变基因做PCA。代码大概长这样（R语言，大家熟悉）：
# 假设 expr 是基因表达矩阵，rows=基因，cols=样本
# batch_vec 是批次标签向量
pca_res <- prcomp(t(expr), scale. = TRUE)
pca_df <- as.data.frame(pca_res$x)
pca_df$batch <- batch_vec
library(ggplot2)
ggplot(pca_df, aes(PC1, PC2, color = batch)) +
geom_point(size = 2, alpha = 0.7) +
theme_minimal()画出来之后，如果不同批次的点像三明治一样分层（PC1就把批次分开了），那基本可以定性为“批次效应显著”。不过这里有个坑——有时候生物学差异也很大，PC1可能是处理组vs对照组，而不是批次。所以我通常会在图上把生物学分组也着色一遍：
ggplot(pca_df, aes(PC1, PC2, color = condition, shape = batch)) +
geom_point(size = 2, alpha = 0.7)3.1.2 t-SNE：局部结构更敏感
PCA是线性方法，有时候批次效应是非线性的。这时候t-SNE（t-distributed Stochastic Neighbor Embedding）就派上用场了。它擅长保留局部结构——说白了，就是让相似的点靠得更近，不相似的点拉得更远。
用Python的scikit-learn可以快速搞定：
from sklearn.manifold import TSNE
import matplotlib.pyplot as plt
# expr 是样本×基因的数组（PCA之后通常取前50个PC作为输入）
tsne = TSNE(n_components=2, perplexity=30, random_state=42)
tsne_out = tsne.fit_transform(expr_pca)
plt.scatter(tsne_out[:,0], tsne_out[:,1], c=batch_colors)
plt.xlabel('t-SNE 1')
plt.ylabel('t-SNE 2')我记得有一次项目，PCA上批次混在一起看不太清楚，我改用t-SNE，立刻就把两批样本分成了两个团——同一处理组在不同批次间距离远大于同一批次内不同处理组。这就是典型的“批次效应大于生物学效应”，需要严重警惕。
3.2 统计检验：给“眼见”一个量化证据
可视化看到分离后，我们要用数字说话。两个常见指标：ANOVA（方差分析）和 R²（方差解释比例）。
3.2.1 ANOVA：检验批次间均值是否有差异
对每个基因，我们可以做单因素ANOVA，看批次这个因子能否显著解释该基因的表达变异。P值<0.05的基因越多，批次效应越普遍。
# R语言对每个基因做ANOVA
p_values <- apply(expr, 1, function(gene) {
fit <- lm(gene ~ batch + condition, data = meta)
anova(fit)["batch", "Pr(>F)"]
})
# 有显著批次效应的基因比例
prop_sig <- mean(p_values < 0.05, na.rm = TRUE)
print(prop_sig)不过要注意：ANOVA对样本量敏感。如果某个批次只有2~3个样本，ANOVA可能非常不靠谱。我一般会把批次内的样本量少于5的剔除掉，或者改用置换检验。
3.2.2 R²：批次解释了多少方差
R²表示批次因子能解释该基因表达变异的比例。比P值更直观——R²越高，批次效应越“大”。比如R²=0.6意味着这个基因60%的变异归因于批次，那它基本不能用于生物学分析。
# 对每个基因计算批次解释的R²
r_squared <- apply(expr, 1, function(gene) {
fit_full <- lm(gene ~ batch + condition, data = meta)
fit_reduced <- lm(gene ~ condition, data = meta)
var_full <- sum(residuals(fit_full)^2)
var_reduced <- sum(residuals(fit_reduced)^2)
1 - var_full / var_reduced
})
summary(r_squared)| 检测方法 | 优点 | 缺点 | 我常用的场景 |
|---|---|---|---|
| PCA | 快速、无参数、可同时看宏观结构 | 线性假设，可能遗漏非线性效应 | 初始筛查首选 |
| t-SNE | 擅长展示局部聚类，能捕捉非线性批次效应 | 随机性大、计算慢、难解释全局意义 | 验证PCA发现的异常聚类 |
| ANOVA | 每个基因有统计量，可做全局分布分析 | 对样本量均匀性敏感 | 当我需要“确凿P值”时 |
| R² | 直接度量效应大小，不受样本量影响 | 对离群点敏感 | 我用来设定阈值（比如R²>0.2的基因太多就要校正） |
3.3 综合判断流程：一张图说清楚
你可能觉得上面这些方法各说各话，怎么整合？我曾经踩过坑：只看了PCA觉得没问题，结果ANOVA跑出来一大半基因显著相关批次，最后t-SNE又发现了暗藏的批次分离。所以后来我总结了一个标准流程：
① 先做PCA看全局，标注批次颜色
② 如果PCA有疑问（比如分离不明显），补一个t-SNE
③ 对每个基因做ANOVA+R²，画分布图：R²密度分布 + 显著比例
④ 若显著基因比例 > 20% 或 R²中位数 > 0.1 —— 启动校正
⑤ 如果上述指标正常但仍有怀疑，做置换检验（shuffle批次标签）来评估偶然性
为了让你更直观地理解这个流程，我画了张框架图：
看到这张图你应该就明白了：可视化负责“定性”，统计检验负责“定量”。两关都过了，才下结论。如果可视化看着没问题但统计检验一堆显著，别急着信视觉——我遇到过样本数足够大时，ANOVA会检出一堆微小差异，但R²只有0.01。这种情况下，我会放宽标准：只看R²>0.05的基因占比。