📚 课程目录第3章批次效应检测方法

3. 批次效应检测方法：PCA、聚类热图、统计检验

好，聊完了批次效应的“长相”和“危害”，咱们得聊聊怎么把它揪出来。你手里拿着一批多组学数据，心里犯嘀咕：“这里面到底有没有批次效应？严重不严重？” 别急，检测方法就那么三板斧，我挨个给你讲讲。

我个人习惯的流程是：先用PCA（主成分分析）快速扫一眼，再用聚类热图验证一下，最后用统计检验给个量化指标。这三步走完，批次的“妖”基本就无处遁形了。

3.1 主成分分析 (PCA) — 最快的一双眼睛

PCA说白了，就是把高维数据压缩到几个主要方向上，方便我们在二维或三维空间里观察样本的分布。你想想看，如果批次效应很强，那不同批次的样本在PCA图上会自然地聚成几个“小团伙”，而不是均匀混合在一起。

我在项目里见过不少新手，拿到PCA图第一眼就盯着看有没有离群点。其实更重要的是看颜色——通常我们按批次给样本上色。如果每种颜色的点都扎堆抱团，那大概率是有问题的。

# Python 示例：用 sklearn 做 PCA 并可视化
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

# 假设 data 是基因表达矩阵（行是基因，列是样本）
# batch_labels 是批次标签
scaler = StandardScaler()
data_scaled = scaler.fit_transform(data.T)  # 转置为样本x基因

pca = PCA(n_components=2)
pcs = pca.fit_transform(data_scaled)

plt.figure(figsize=(8,6))
for batch in set(batch_labels):
    idx = batch_labels == batch
    plt.scatter(pcs[idx, 0], pcs[idx, 1], label=f'Batch {batch}', alpha=0.7)
plt.xlabel('PC1 ({:.1f}%)'.format(pca.explained_variance_ratio_[0]*100))
plt.ylabel('PC2 ({:.1f}%)'.format(pca.explained_variance_ratio_[1]*100))
plt.legend()
plt.title('PCA: 样本按批次着色')
plt.show()

我的小技巧： 如果PC1和PC2的方差解释率加起来不到30%，可以考虑多画几个PC的组合（比如PC2 vs PC3）。另外，可以尝试用不同形状标记不同的生物学分组，看看批次效应是否真的独立于生物学差异。

记得有一次，我帮一个合作实验室看他们的RNA-seq数据。PCA一看，第一批和第二批样本完美地分在了PC1的两边，而生物学条件在PC2上才有微弱差异。我跟他们说：“兄弟，你们这个批次效应比生物学效应还大，必须先校正。” 这就是PCA最直观的价值。

3.2 聚类热图 — 从全局看批次的“指纹”

PCA能让你快速看到样本的宏观分布，但它的缺点也很明显——它只展示了前几个主成分，可能会遗漏一些藏在更深维度的信息。这时候，聚类热图就派上用场了。

我通常会选一批变异最大的基因（比如前1000或前2000个），然后做样本和基因的双向聚类。热图的颜色代表表达量的高低。如果批次效应很强，你会发现样本的聚类树状图会先按批次分开，而不是按生物学条件分开。这比PCA更“实锤”。

  关键指标： 查看热图上方的样本聚类树。如果同一批次的样本在树状图上形成独立的“枝”，且枝的长度远小于不同批次之间的枝长，那就说明批次效应主导了样本的相似性。

# Python 示例：用 seaborn 绘制聚类热图
import seaborn as sns
import numpy as np

# 假设 data 是基因表达矩阵，取 top 1000 高变异基因
top_var_genes = np.argpartition(np.var(data.values, axis=1), -1000)[-1000:]
data_top = data.iloc[top_var_genes, :]

# 创建批次注释
batch_colors = [{'Batch': batch} for batch in batch_labels]

# 绘制热图
sns.clustermap(data_top, 
               z_score=0,  # 按行标准化
               col_colors=batch_colors,  # 列标签颜色
               figsize=(10, 12),
               cmap='vlag',
               xticklabels=False,
               yticklabels=False)

避坑指南： 我曾经犯过一个低级错误——没有对基因进行标准化就直接做热图。结果高表达基因的颜色把整体对比度都拉垮了，根本看不出批次模式。记得要么用z-score，要么对每行基因做最大最小值归一化。

3.3 统计检验 — 给批次效应一个“数字”

图像会骗人，但数字不会。PCA和热图可以给你直观的感受，但如果我们想说服审稿人或者自己，需要一个量化的评价指标。常见的统计检验方法有几种，我挑最常用的说说。

3.3.1 ANOVA / Kruskal-Wallis 检验

对于每个基因，我们都可以检验它在不同批次之间的表达量是否存在显著差异。如果大多数基因的p值都小于0.05，那就有理由怀疑存在广泛的批次效应。

# R 示例：对每个基因进行 ANOVA 检验
# 假设 expr 是基因表达数据框，batch 是批次向量
p_values <- apply(expr, 1, function(x) {
  summary(aov(x ~ batch))[[1]][["Pr(>F)"]][1]
})
# 或者用 Kruskal-Wallis（非参数）
p_values_kw <- apply(expr, 1, function(x) {
  kruskal.test(x ~ batch)$p.value
})
# 看看有多少基因的 p < 0.05
sum(p.adjust(p_values, method='fdr') < 0.05, na.rm=TRUE)

一点个人看法：我其实不太建议仅仅依赖ANOVA或KW检验。因为当样本量很大时，极微小的表达差异也可能被判定为“显著”，而这点差异对下游分析可能完全没有影响。所以更推荐下面这个方法。

3.3.2 PVCA (Principal Variance Component Analysis)

PVCA是我个人比较偏爱的方法。它能把表达量的总方差拆解成不同成分——有多少是批次引起的，有多少是生物学变量引起的，有多少是随机误差。这样你就能很直观地说：“批次效应占了总方差的15%”。如果这个比例超过了生物学分组所占的比例，那就需要介入校正了。

# R 示例：使用 PVCA 包（简化版本）
# 安装：BiocManager::install("pvca")
library(pvca)

# pData 是样本信息，包含 batch 和 condition
# expr 是表达矩阵
pct_threshold <- 0.6  # 过滤低变异基因

pvcaObj <- pvcaBatchAssay(expr, 
                          pData,
                          ~batch + condition,
                          pct_threshold)
# 输出结果
barplot(pvcaObj$dat,
        xlab = "方差成分",
        ylab = "方差贡献比例",
        col = c("steelblue", "coral", "grey"),
        las=2)

  我的判断标准： 如果批次效应贡献的方差比例超过10%，且大于生物学分组贡献的比例，那基本上就是“红色警报”了。

3.4 一张图看懂整个检测流程

说了这么多，我把整个检测流程整理成了一张结构图。你可以把它当作你的“排雷”路线图。

3.5 综合对比：三种方法怎么选？

我最后用一个表格总结一下这三种检测方法的差异和使用场景。方便你查漏补缺。

方法	优点	缺点	适用场景
PCA	速度快、直观、一眼看全局	只展示前几个PC，可能遗漏信息	初步筛查，快速判断有无问题
聚类热图	展示全貌，能看出样本层级关系	基因多时图很大，难以精细解读	验证PCA发现，或样本量较小时使用
ANOVA/PVCA	量化指标，可重复性好	需要一些统计基础，PVCA计算稍慢	写文章、审稿、正式报告必备

一句忠告： 不要光盯着好看的PCA图就下结论。我曾经见过一批数据，PCA图上样本混合得非常好，但聚类热图一画，树状图劈成两半——原来批次效应藏在第四主成分上了。所以，三板斧都抡一遍再拍板，不亏。

这就是批次效应检测的核心方法。你学会了吗？下一步，当你确认批次效应存在且影响较大时，就该考虑怎么消除它了。不过那部分我们后面再聊。