第3章 批次效应识别：用PCA、聚类图看批次问题

上一章咱们把数据装进了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")

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

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

我记得有一个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和聚类热图两条路，最后综合判断是否要进入批次校正。你可以在实验室的白板上画一遍，脑子里就有谱了。

批次效应识别流程 标准化表达数据 PCA 降维可视化 样本距离聚类热图 样本是否按批次分离？ PC1/PC2 颜色块明显分层 树状图是否按批次分支？ 同批次样本率先聚为一簇 综合判断：是否存在批次效应？ 是 → 进入校正环节；否 → 直接差异分析 两侧结果一致才做判断。若批次混杂条件，需额外设计验证实验。 PCA 和热图只是“肉眼观察”，最终效应大小见第4章方差分解。

一点小补充：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，结果任何批次效应都被抹平了，那就是自己骗自己。

← 上一章 📖 返回目录 下一章 →