📚 课程目录第4章批次效应的诊断——PCA图

4. 批次效应的诊断——PCA图：如何通过主成分分析初步判断批次结构及其与生物学因素的混杂

拿到RNA-seq数据的第一件事，我一般不会直接跑差异分析。我会先看看PCA图。

为什么？因为PCA图能最快告诉你，你的数据里藏着什么“脏东西”。批次效应、组织差异、处理效果，全都会在PCA图上留下痕迹。你需要的不是一双慧眼，而是一套方法论。

4.1 PCA到底在做什么

说白了，PCA就是降维。把几千个基因的表达量，压缩成几个主成分（PC1、PC2...），让你能在二维或者三维空间里看到样本之间的关系。

这里有两个核心直觉：

PC1（主成分1）：解释数据变异最大的方向。通常代表你最关心的生物学因素（比如肿瘤 vs 正常）。
PC2（主成分2）：解释剩余变异中最大的方向。这里就经常藏着批次效应。

我个人判断批次效应的惯例是——先不管PC1，直接看PC2、PC3。如果样本在PC2/PC3上按批次分开，而PC1按生物学因素分开，那就比较干净。如果PC1都被批次占了，那麻烦就大了。

小技巧： 绘制PCA时，我习惯标注颜色为“批次”，形状为“实验条件”。如果颜色聚成团，那就是批次效应；如果形状和颜色混在一起，那就是理想状态。

4.2 如何从PCA图上读信息

我整理了一个表格，你可以对照着看自己的PCA图：

PCA图表现	可能含义	严重程度	下一步动作
样本按批次明显分离在PC1上	批次效应完全主导了表达变异	严重（必须处理）	使用Combat、limma removeBatchEffect
批次在PC2/PC3上分离，生物学因素在PC1上分离	批次效应存在但不完全掩盖生物学信号	中等（需要处理）	检查混杂程度，考虑回归或匹配
批次和生物学因素在同一PC上混合	混杂（confounding）	高度危险	尝试分层分析、敏感性分析
样本随机分布，批次混杂	批次效应微弱或不存在	低	继续分析，但建议验证
单个批次内的样本离群	可能是技术误差或样品质量问题	局部问题	核实该样本的测序QC指标

我曾经接手过一个项目，PCA图上PC1完美分开了“处理组”和“对照组”。大家都觉得漂亮。结果一查批次表，发现所有处理组都是批次1，所有对照组都是批次2。这就是典型的混杂。PCA图上的分离本质是批次效应，不是生物学效应。那批数据基本废了一半。所以，只看PCA图不够，你得看“批次 vs 生物学因素”的关联表。

4.3 实战：在R里跑PCA诊断

我会用 DESeq2 或者 prcomp 做PCA。下面是我常用的套路：

# 使用DESeq2的vst或者rlog变换后的数据做PCA
library(DESeq2)
library(ggplot2)

# 假设dds是DESeqDataSet对象
vsd <- vst(dds, blind = TRUE)

# 提取变换后的表达矩阵
mat <- assay(vsd)
# 转置，行是样本，列是基因
pca <- prcomp(t(mat), center = TRUE, scale. = FALSE)

# 提取前几个主成分
pca_df <- as.data.frame(pca$x)
pca_df$sample <- rownames(pca_df)
pca_df$batch <- colData(dds)[rownames(pca_df), "batch"]
pca_df$condition <- colData(dds)[rownames(pca_df), "condition"]

# 绘制PC1 vs PC2
ggplot(pca_df, aes(x = PC1, y = PC2, color = batch, shape = condition)) +
  geom_point(size = 3) +
  theme_minimal() +
  labs(title = "PCA: PC1 vs PC2")

这里要注意，blind = TRUE 能让变换不看到你指定的模型设计，否则可能会过度拟合。

避坑指南： 我曾经犯过一个错误，直接用原始count做PCA。结果PC1被高表达基因的文库大小差异主导，根本不是批次或生物学效应。所以一定要用vst、rlog或者CPM/log2标准化后的数据。

4.4 从PC variance占比中找线索

还有一个容易忽略的线索——每个PC解释的方差百分比（也就是 scree plot）。

如果前两个主成分合计占比超过70%，说明你的数据被少数几个因素强烈主导。这时候要警惕，这个“主导因素”到底是生物学因素还是批次效应。更常见的情况是前10个PC的方差分布非常平缓，每个PC占比不到10%。那意味着影响因素很多，批次效应可能只是其中之一。

# 查看每个PC的方差占比
pca_summary <- summary(pca)
pca_var <- pca_summary$importance[2, ] * 100

# 绘制Scree plot
barplot(pca_var[1:20], xlab = "Principal Component", ylab = "Variance Explained (%)",
        col = "steelblue", main = "Scree Plot")

嗯，这里有个规律我观察了很久——如果在实验设计中有意识地平衡批次（比如每个批次里都包含所有条件），那么批次效应通常只会在PC2或PC3之后显现。如果实验设计时完全没有考虑批次，那就大概率在PC1上出现。

4.5 判断混杂：PCA + heatmap组合拳

光看PCA不够，我习惯把PCA结果和批次-条件关联矩阵结合起来看。

比如你可以画一张热图，行是样本，列是批次和条件，颜色代表样本属于哪个批次或哪个条件。然后把热图和PCA的样本顺序对齐。如果看到“条件”列的颜色变化和“批次”列几乎一模一样，那就说明严重混杂。

用图形化方式表达这一章的核心逻辑，我画了一张流程图：

4.6 一个实际案例

我记得有一个公开数据集：GSEXXXXX，包含12个样本，3个批次，每个批次4个样本。我拿到数据后直接做了PCA。

PC1: 32.5% variance
PC2: 18.7% variance
PC3:  9.2% variance

PCA图上：
- PC1（x轴）：样本按批次1/批次2/批次3从左到右排列
- PC2（y轴）：混合无规律

结论：批次效应主导了PC1。

为什么会这样？因为实验设计里批次1全是case，批次2全是control，批次3全是treatment。批次和条件完全混杂。这种情况下，如果不做任何校正，直接做差异分析，结果将是“批次效应”被当成“处理效应”。

怎么办？

第一步：用limma的removeBatchEffect去掉批次效应（前提是批次与条件不完美混杂）。
第二步：再跑一次PCA，看PC1是否变为生物学因素。
第三步：如果依然无法分离，需要考虑匹配分析或者重新设计实验。

关键判断点： 如果批次和生物学因素在PCA图上完美对应（即一张图里颜色和形状分得一模一样），那就说明你无法从统计上区分它们。此时任何校正方法都会引入偏差。

4.7 总结几条实用经验

不要只看PC1/PC2：我见过PC1-3都干净的样本，却在PC4上暴露了批次效应。多看看前10个PC。
记录每次PCA的方差占比：这能帮你量化批次效应的贡献。
准备三张PCA图：一张按批次上色，一张按条件上色，一张按测序批次/文库/lane上色。对比着看，总能发现问题。
注意离群点：单个样本跑偏太远，别急着删。有可能是那个样本的批次标签搞错了——之前一个同事就是在这里栽了跟头。

看到PCA图之后，下一步就该考虑怎么量化批次效应、怎么校正、怎么验证。但那是后面的章节，这里我们先停在“诊断”阶段。记住一个原则：诊断不到位，校正就是白费功夫。