第二讲：批次效应的可视化探索

说起批次效应，我个人的习惯是：先别急着建模，先拿眼睛看看数据长什么样。

你想想，数据里有没有明显的批次分组？不同批次的样本是不是天然分开了？这些用肉眼看比跑统计模型要直观得多。这章我们就来聊聊五种可视化方法——箱线图、PCA、热图、密度图、MDS。每一种都是我的老伙计，在不同项目里帮我揪出过不少批次问题。

2.1 为什么先看可视化？

说白了，可视化是所有分析的“体检报告”。

我遇到过不少项目，分析师吭哧吭哧跑了一大堆统计检验，结果后来发现数据本身就有系统性的批次偏移。那感觉就像你花了一个月研究汽车的发动机，结果发现轮胎是方的。

2.2 箱线图：最直接的表达

箱线图是入门级的批次效应查看工具。它能告诉你：不同批次的样本，表达量分布是不是一致的。

假设我们有一个表达矩阵 expr，加上一个批次标签向量 batch。R的基础绘图就能完成：

# 随机选一个基因，比如第100个
gene_idx <- 100
boxplot(expr[gene_idx, ] ~ batch,
        main = "Gene 100 表达量分布 - 按批次",
        xlab = "批次", ylab = "表达量",
        col = c("#3498db", "#e74c3c", "#2ecc71"))

如果不同批次的中位数明显不一样，或者四分位距差别很大——嗯，那就要留个心眼了。

怎么看？

• 中位数差异大：批次效应很可能是真实存在的。
• 离散度不一致：有的批次波动大，有的批次波动小，这也是问题。
• 异常点集中在某批次：这个批次可能有技术问题。

2.3 PCA：降维看全局

PCA是我最常用的方法。它能把高维数据压缩到两三个维度，然后看看批次是不是“聚成一团”。

用R的 prcomp 函数做PCA很简单：

# 转置：行为样本，列为基因
pca_result <- prcomp(t(expr), scale. = TRUE)

# 提取前两个主成分
pca_data <- as.data.frame(pca_result$x[, 1:2])
pca_data$batch <- batch

# 绘图
library(ggplot2)
ggplot(pca_data, aes(x = PC1, y = PC2, color = batch)) +
  geom_point(size = 3, alpha = 0.8) +
  theme_minimal() +
  labs(title = "PCA：样本按批次着色")

关键判断标准：

• 如果相同批次的样本在PCA图上紧挨在一起，不同批次明显分开——批次效应就相当明显了。
• 如果同批次样本也散得很开，不同批次混在一起——那批次效应可能不严重。
• PC1和PC2解释方差的比例要留意。如果头两个主成分只解释了不到20%，那这张图的可信度就要打点折扣。

我在一个单细胞RNA-seq项目里遇到过，PCA图上批次完全分成了两群。当时第一个反应就是：“完了，得好好校正了。”

2.4 密度图：分布形状的对比

密度图可以看成是箱线图的连续版本。它能更精细地显示表达量的分布形状。

# 随机选一个样本，在每个批次里各抽取一个
sample_ids <- tapply(1:ncol(expr), batch, function(x) x[1])

# 提取这三个样本的表达量
library(ggplot2)
df <- data.frame(
  expression = expr[1, sample_ids],
  batch = factor(batch[sample_ids])
)

ggplot(df, aes(x = expression, fill = batch)) +
  geom_density(alpha = 0.5) +
  theme_minimal() +
  labs(title = "密度图：同一基因在不同批次中的分布")

密度图能告诉你什么？

• 峰的位置偏移：不同批次的分布中心不同。
• 峰的形状不同：有的批次是双峰，有的是单峰，这通常说明批次间存在技术差异。
• 尾部差异：有的批次高表达基因多，这也是批次效应的信号。

2.5 热图：批次效应的“全貌”

热图结合了聚类和颜色映射，能看到整体模式。我用它来观察批次是否在聚类树中形成独立分支。

# 计算样本间的相关性矩阵
cor_matrix <- cor(expr, method = "spearman")

# 用pheatmap绘制
library(pheatmap)
annotation <- data.frame(batch = factor(batch))
rownames(annotation) <- colnames(expr)

pheatmap(cor_matrix,
         annotation_col = annotation,
         annotation_row = annotation,
         main = "样本相关性热图 - 按批次标注",
         color = colorRampPalette(c("#3498db", "white", "#e74c3c"))(50))

怎么解读？

• 同一批次的样本在热图里形成明显的色块——批次效应存在。
• 聚类树甩出明显的批次分组——很可能是批次效应主导了聚类。
• 如果批次标签和聚类结果基本无关——那你可以松一口气。

2.6 MDS：PCA的坚强替补

MDS（多维尺度分析）和PCA的思路类似，都是降维。区别在于MDS直接基于样本间的距离矩阵，不要求线性假设。

# 计算欧氏距离矩阵
dist_matrix <- dist(t(expr))

# 执行经典MDS
mds_result <- cmdscale(dist_matrix, k = 2)
mds_data <- as.data.frame(mds_result)
mds_data$batch <- batch

ggplot(mds_data, aes(x = V1, y = V2, color = batch)) +
  geom_point(size = 3, alpha = 0.8) +
  theme_minimal() +
  labs(title = "MDS：基于距离矩阵的批次效应可视化")

什么时候用MDS？

• 数据不符合线性假设时，MDS往往比PCA更靠谱。
• 想要使用非欧氏距离（比如曼哈顿距离、相关系数距离）时，MDS更灵活。
• PCA和MDS的结论相互印证，会让你的分析更有说服力。

2.7 五种方法的对比表

方法	侧重点	优点	局限
箱线图	单个基因的分布差异	简单快速，一目了然	只能看单个基因，多基因需逐个检查
PCA	全局样本的聚类模式	降维效果好，方差解释度高	对异常值敏感，需要适当标准化
密度图	分布形状的精细对比	能发现双峰等复杂分布	一次只能看一个基因/一个样本
热图	整体相关性和聚类结构	直观展示批次相关模式	样本量大时视觉过载
MDS	基于任意距离的降维	灵活性高，不限于线性假设	对距离度量敏感，需要选择合适距离

2.8 核心流程：一张图说清楚

下面这张SVG图汇总了这五种可视化方法的使用流程。说白了，就是从原始数据出发，选择合适的可视化方式，最终判断批次效应是否存在。

2.9 避坑指南

我曾经犯过的几个毛病：

• 只看一种图就下结论：PCA说没问题，但热图可能反映出细微的模式。多角度交叉验证是王道。
• 忽略了标准化：不做标准化直接跑PCA，结果全是表达量高的基因在主导。那场面，就像你在菜市场里用扩音器喊话——根本听不到正常声音。
• 样本量太少：每个批次只有两三个样本时，任何可视化都很难给出可靠判断。这时宁可多收集数据，也不要强行分析。
• 忘记检查离群样本：一个离群样本足以把PCA图扭曲得面目全非。先过滤再可视化，这是个好习惯。