📚 课程目录第3章批次效应识别

3. 批次效应识别

上一章咱们把数据装进了DESeq2的容器里，但先别急着跑差异分析。我个人的习惯是，在跑任何正式模型之前，一定、一定先用眼睛看看数据——说白了就是可视化检查。这一步能帮你快速判断数据里有没有“不对劲”的批次结构，也能让你对后续校正有没有效果心里有个底。

今天就聊两样最常用的招：主成分分析（PCA） 和 样本聚类热图。这两招用好了，批次问题基本藏不住。

为什么是PCA和聚类图？

你想想看，批次效应说白了就是一种系统性的技术偏差——同一个批次的样本更相似，不同批次的样本被强行拉开。这个“拉开”的动作，天然就会在数据的最主要变异方向（PC1、PC2）上留下痕迹。PCA正好是干这个的：把高维表达量压缩到两三个主成分上，让你一眼扫过去就能看到样本是按生物学分组聚在一起，还是按批次编号抱团。

聚类热图则是从距离角度看的：如果同一批次的样本在树状图里总是先扎堆，哪怕它们分属不同的处理组，那批次效应就是实锤了。我刚开始做RNA-seq时，就吃过一个亏：表面上看差异基因列表挺漂亮，结果审稿人回一句“请展示批次效应评估”，我补了个PCA图，发现样本完全按测序批次分开。嗯，从那以后我再也不敢跳过这步了。

小提示： PCA + 热图是“双保险”。PCA看宏观分离趋势，热图看微观的样本间相似度。两者结论一致时，你基本可以放心判断是否存在批次效应。

用PCA扫描批次痕迹

DESeq2里自带了一个方便函数：plotPCA。它直接基于vst或rlog变换后的数据算主成分，然后给你画个散点图。我一般会用vst（variance stabilizing transformation），因为它对低表达基因的处理比rlog更稳定。

标准的操作流程长这样：

# 先做方差稳定化变换
vsd <- vst(dds, blind = FALSE)

# 然后画PCA，用color标注批次
plotPCA(vsd, intgroup = "batch") +
  geom_point(size = 3) +
  ggtitle("PCA: batch distribution")

注意我用了intgroup = "batch"，这样点是按批次上色的。如果你不放心，可以再加个形状参数shape = "condition"，这样同一个批次内不同处理条件也能区分开。

怎么看结果？

好信号：不同条件的样本散布在同一片区域，不同批次的样本均匀混合。PC1和PC2上解释的方差比例各自不超过30%（批次效应强时，PC1常能拉到40%以上）。
危险信号：同一批次样本挤在一起，不同批次有明显的分离边界，哪怕PC1只解释了20%，只要颜色块边界清晰，就要警惕。
复杂情况：批次与生物学条件完全重合（比如处理组都是batch1，对照组都是batch2），这时PC图上的分离可能是真实差异也可能是批次效应，需要结合先验知识或设计实验来甄别——这个坑我踩过不少次。

注意： 当批次和实验条件完全混杂时，blind = FALSE 做vst可能会导致校正过度。这种情况下，先用blind = TRUE试试，或者干脆用原始count标准化后做PCA，只作为参考，别直接用来差异分析。

聚类热图：另一种视角

PCA看的是整体投影，聚类热图则保留了更多的距离信息。我项目里常用pheatmap包来画样本间的欧氏距离热图，配合树状图一目了然。

代码示例：

library(pheatmap)
library(RColorBrewer)

# 计算样本间距离
sampleDists <- dist(t(assay(vsd)))
sampleDistMatrix <- as.matrix(sampleDists)
rownames(sampleDistMatrix) <- vsd$sample_id
colnames(sampleDistMatrix) <- vsd$sample_id

# 准备批次注释
annotation_col <- data.frame(
  batch = vsd$batch,
  condition = vsd$condition,
  row.names = vsd$sample_id
)

colors <- colorRampPalette(rev(brewer.pal(9, "Blues")))(255)
pheatmap(sampleDistMatrix,
         clustering_distance_rows = sampleDists,
         clustering_distance_cols = sampleDists,
         annotation_col = annotation_col,
         annotation_row = annotation_col,
         col = colors,
         main = "Sample distance heatmap")

看热图时，我重点关注两件事：

批次是否先合并：观察树状图底部的分支，如果同一批次的样本总是率先聚成小簇，哪怕条件不同，就有批次问题。
对角线外的颜色块：理想情况下，不同批次、相同条件的样本之间距离应该比较近（颜色偏浅蓝），如果批次的差异导致同条件、不同批次的样本距离都很深（深蓝），那你就该考虑校正了。

我记得有一个Allo-HCT数据集，做免疫细胞RNA-seq，样本来自三个测序批次、四个供体组。第一次PCA图里PC1/PC2完全按批次劈开，热图里batch1和batch2各成一团。当时我以为数据废了，后来发现只是批次效应，用DESeq2的~ batch + condition校正后就正常了。这就是典型的“先诊断，再开药”。

避开常见误区

每个人都会犯的错，我列几个你重点记一下：

误区	真实情况	我的建议
PCA图看起来没问题就认为没有批次效应	可能PC1/PC2看不到，但高维PC上有	同时看PC3和PC4，或者用UMAP辅助
批次和条件重合时仍用PCA判断	无法区分，PCA会误导	做实验设计分离，或者在模型里加入交互项
只看样本间聚类热图，忽略基因层面	全局样本距离可能掩盖批次对单个基因的影响	补充检查top变异基因的热图
PCA图只给publization用，自己不看	后期校正后对比PCA才知道效果	务必保留校正前后的PCA截图

核心判断原则：

如果在PCA中，同一批次的样本有5个以上聚集在一个小区域，而其他批次在另外的区域，且这种分离与实验条件无关 → 强批次效应。
如果在聚类热图中，树状图的第一级（或前两级）分支按批次分开，且样本注释条颜色块能明显看到批次行列 → 强批次效应。
如果两者都只有轻微倾向，且PC1解释方差<20%，可以先不校正，继续分析。

图示：批次效应识别流程

下面这张图概括了今天讲的核心逻辑：从标准化后的数据出发，走PCA和聚类热图两条路，最后综合判断是否要进入批次校正。你可以在实验室的白板上画一遍，脑子里就有谱了。

一点小补充：R包之外的技巧

有些时候，你手头的数据不是DESeq2对象，而是已经标准化好的矩阵。没关系，PCA照样能画：

# 假设 exp_mat 是基因×样本矩阵（已log2标准化）
pca <- prcomp(t(exp_mat), center = TRUE, scale. = FALSE)
pca_df <- as.data.frame(pca$x)
pca_df$batch <- batch_info
ggplot(pca_df, aes(PC1, PC2, color = batch)) + geom_point()

还有一个小细节：scale. = FALSE，因为表达量的尺度差异就包含了信息，不需要额外缩放。我曾经见过有人在这步用了scale. = TRUE，结果任何批次效应都被抹平了，那就是自己骗自己。