📚 课程目录第3章批次效应的可视化探索

3. 批次效应的可视化探索：PCA、热图、箱线图

数据预处理完了，也做了初步的差异分析。下一步要做的，就是亲眼看看这批次的“鬼”到底藏在哪里。我有个习惯：不先做可视化，绝不上ComBat校正。因为你需要知道问题有多严重，是系统性偏移还是个别样本捣乱，这直接影响你的校正策略。

可视化不是走形式，它是你跟数据之间的“坦白局”。今天咱们就聊聊三种最趁手的工具：PCA、热图和箱线图。

3.1 为什么一定要先可视化？

说白了，统计检验只能告诉你“有没有”，但可视化能告诉你“是什么样”。我有一个特别深的教训：有一次拿到一个肿瘤转录组数据集，用批次效应检测的统计方法跑出来p值巨大，以为没事。结果PCA一画，两批样本在PC1上完美分离——后来发现是批次标签搞反了。嗯，从那以后，我宁愿多花半小时画图，也绝不偷懒。

可视化要回答三个核心问题：

批次间是否有明显的分布偏移？
偏移是全局性的，还是只在特定基因/样本上？
批次效应和生物学变异，哪个在主导数据变异？

带着这些问题，我们来逐个击破。

3.2 PCA：看全局，抓主体变异

我的个人偏好： 做PCA之前，先把数据做一下中心化缩放（scale = TRUE）。虽然PCA内部会做，但提前标准化能让离群样本不那么“霸道”，不至于一个样本拖垮整个主成分。

PCA的思路很简单：把高维数据投影到几个最能解释变异的方向上。如果批次效应很强，你会发现样本点按批次簇拥在一起，即使前两个主成分只占总变异的20%，也能一眼看出问题。

为什么PCA图能让人一眼看出批次问题？因为批次效应通常是全局性的，它会影响很多基因的表达模式，所以会在前几个主成分上暴露无遗。

实战中，我通常会画两种PCA图：

按批次着色 —— 看批次聚类情况。
按生物学变量着色 —— 看生物学分组是否与批次混在一起。

来看看R语言的实现代码：

library(ggplot2)
library(ggrepel)

# 假设expr_matrix是表达矩阵，行是基因，列是样本
# batch_vector是批次标签向量
pca_res <- prcomp(t(expr_matrix), scale. = TRUE, center = TRUE)
pca_df <- data.frame(PC1 = pca_res$x[,1], PC2 = pca_res$x[,2], Batch = batch_vector)

ggplot(pca_df, aes(x = PC1, y = PC2, color = Batch)) +
  geom_point(size = 3, alpha = 0.8) +
  stat_ellipse(level = 0.95) +
  labs(title = "PCA by Batch") +
  theme_minimal()

踩坑提醒： 我曾经遇到过一个情况——PCA图按批次没有明显分离，但按另一个生物学因子（比如疾病亚型）分离得特别好。千万别高兴太早，继续往下看热图！因为如果批次和生物学因子是完全混杂（confounded）的，PCA看不出来，但热图会暴露。

如果要把PCA结果组织成一个更容易对比的表格，可以参考下面这个思路：

观察指标	无批次效应	有批次效应
PC1/PC2 散点分布	样本均匀混合，无聚类	不同批次样本成簇分离
前几个PC的解释方差	变异主要由生物学因素解释	批次因子贡献大
与生物学分组对比	生物学分组清晰可见	生物学分组被批次掩盖

3.3 热图：看细节，抓特异性模式

PCA看的是全局大势，热图看的是细节肌理。尤其是当批次效应只影响少数基因（比如管家基因）时，PCA可能完全无感，但热图上会出现明显的“竖条纹”。

你想想看：热图的行是基因，列是样本，颜色深浅代表表达量高低。如果样本按批次排列后，出现了一条横跨整个列的色带——那就是批次效应在“显形”。

我个人的习惯是：先做一次层次聚类（hierarchical clustering），看看样本是否按批次聚类，而不是按实验分组聚类。

library(pheatmap)

# 选择变异度最高的前500个基因
gene_vars <- apply(expr_matrix, 1, var)
top_genes <- names(sort(gene_vars, decreasing = TRUE))[1:500]

# 样本注释数据框
annotation_col <- data.frame(Batch = batch_vector, Group = group_vector)
rownames(annotation_col) <- colnames(expr_matrix)

pheatmap(expr_matrix[top_genes, ],
         scale = "row",
         annotation_col = annotation_col,
         show_rownames = FALSE,
         show_colnames = FALSE,
         cluster_cols = TRUE,
         main = "Top 500 Variable Genes Heatmap")

  核心判断标准： 如果列（样本）的聚类树先按批次分叉，再按生物学分组分叉——说明批次效应的强度已经超过了生物学效应。这时候就要立刻启动校正流程了。

记得有一次，我处理一个代谢组数据。PCA图看不出任何问题，PC1和PC2上样品点完全混在一起。但我坚持画了热图，结果发现：在样本列的树状图上，有两棵独立的大分支，且分支内的样本全是同一个批次的。原来这批数据的批次效应主要体现在一群代谢物上，而这些代谢物的变异性很小，在PCA里被“淹没”了。从那以后，热图成了我检验批次效应的“照妖镜”。

3.4 箱线图：看分布，量化偏移

箱线图是三种方法中最“直接”的：每个样本所有基因的表达值分布画成一个箱子，看看不同批次的箱子是不是整体偏高或偏低。

我个人认为，箱线图最适合用在校正前后对比，但在检测阶段，它有另一个妙用——找离群样本。

library(reshape2)
library(ggplot2)

# 计算每个样本表达量的中位数
sample_medians <- apply(expr_matrix, 2, median)
plot_df <- data.frame(Sample = names(sample_medians),
                      MedianExpr = sample_medians,
                      Batch = batch_vector)

ggplot(plot_df, aes(x = Batch, y = MedianExpr, fill = Batch)) +
  geom_boxplot(outlier.color = "red", outlier.shape = 19) +
  labs(y = "Median Expression Level", title = "Boxplot of Median Expression by Batch") +
  theme_minimal()

你可能会问：箱线图到底看什么？三个字——错位感。

如果Batch1的中位数线在5.2，Batch2的中位数线在4.8。这个0.4的差距，放在整个表达量尺度上，可能意味着批次间的系统性偏差。
如果某一个批次内有样本中位数显著偏离该批次的整体水平——那这个样本可能就是技术性离群点，需要考虑剔除。

我的经验： 箱线图虽然直观，但受限于它只看中位数和四分位距。如果批次效应只影响高表达的基因（比如上调基因），箱线图可能不敏感。所以我会配合PCA和热图一起看，形成“三角验证”。

3.5 本章小结：三种图如何组合使用？

我习惯给可视化做一个简单的“诊断矩阵”：

可视化方法	核心观察点	优势	局限
PCA	主成分上样本的聚类情况	快速捕捉全局模式	对低贡献变异不敏感
热图	样本聚类树、列色带模式	暴露特异性批次效应	高度依赖基因选择
箱线图	各样本表达量中位数分布	量化整体偏移、发现离群	对局部效应不敏感

这三兄弟必须一起用。别觉得麻烦——你想想看，如果漏掉了批次效应就上了ComBat，不仅校正效果打折扣，还可能引入新的假象。这种坑，我踩过，不希望你踩。

下面这张图，把本章的检测逻辑串起来：

当你通过这三种图的组合验证，确认批次效应确实存在，并且了解它的具体模式后，下一步才是真正的校正操作。但那是第4章的事了。在这之前，把可视化做到位，让你自己闭上眼都能说清楚批次效应长什么样——这才是分析师的底气。