4. 主成分分析（PCA）在批次效应检测中的应用
主成分分析——PCA，说实话这是我做组学数据分析时最常用的工具之一。你想想看，高维数据动辄几万个基因，肉眼根本看不出分布。PCA 能把维度压缩到两三个主成分，然后直接画个散点图，批次效应往往一眼就暴露了。
我个人习惯呢，在正式做差异分析之前，必先跑一遍 PCA。为什么？因为批次效应如果很严重，后续的统计推断全是扯淡。这一讲我就带你看看 PCA 到底是怎么帮我们“抓出”批次效应的。
4.1 PCA 的核心思想
说白了，PCA 就是寻找数据中方差最大的方向。第一个主成分（PC1）是数据投影后方差最大的方向，第二个主成分（PC2）是正交于 PC1 且方差次大的方向，依此类推。
当批次效应是数据中主要的变异来源时，它往往会被捕获到前几个主成分里。我在一个多中心 RNA-seq 项目里遇到过——不同测序中心的数据混在一起，PC1 直接把样本按中心分成了两簇，批次效应明晃晃地摆在那里。
4.2 PCA 检测批次效应的标准流程
流程不复杂，但每一步都有坑。我一般按下面几步走：
- 数据预处理： 标准化（Z-score）或 log2 转换，去除低表达/低变异基因。
- 运行 PCA： 用 prcomp() 或 sklearn.decomposition.PCA。
- 可视化： 用前两个或三个主成分着色，批次信息用颜色表示。
- 判断： 如果批次间样本分离，则存在批次效应。
- 量化： 计算各主成分与批次变量的相关程度。
嗯，这里要注意：标准化很重要。如果基因表达量差异巨大（比如有些基因是 1000，有些是 10），PCA 会被高表达基因主导，批次效应就可能被掩盖。我在做单细胞数据时吃过这个亏，后来老老实实先做 scale。
4.3 代码实战：用 R 和 Python 做 PCA 检测
R 版本（基于 DESeq2 的 rlog 或 vst 数据）
# 假设 expr 是 log2 归一化后的表达矩阵，行为基因，列为样本
# batch 是批次因子向量
pca_result <- prcomp(t(expr), center = TRUE, scale. = FALSE)
percentVar <- round(100 * (pca_result$sdev^2 / sum(pca_result$sdev^2)), 1)
# 提取前两个主成分
pca_df <- data.frame(
PC1 = pca_result$x[, 1],
PC2 = pca_result$x[, 2],
Batch = batch
)
# 用 ggplot2 画图
library(ggplot2)
ggplot(pca_df, aes(x = PC1, y = PC2, color = Batch)) +
geom_point(size = 3, alpha = 0.8) +
labs(x = paste0("PC1 (", percentVar[1], "%)"),
y = paste0("PC2 (", percentVar[2], "%)")) +
theme_minimal()Python 版本
import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt
# expr_df: 基因×样本 的 DataFrame
# 转置为 样本×基因
data = expr_df.T.values
scaler = StandardScaler()
data_scaled = scaler.fit_transform(data)
pca = PCA(n_components=2)
pc_scores = pca.fit_transform(data_scaled)
# 方差百分比
var_exp = pca.explained_variance_ratio_ * 100
plt.figure(figsize=(6, 5))
for batch_name, group in df_batch.groupby('batch'):
idx = group.index
plt.scatter(pc_scores[idx, 0], pc_scores[idx, 1],
label=batch_name, s=30, alpha=0.7)
plt.xlabel(f'PC1 ({var_exp[0]:.1f}%)')
plt.ylabel(f'PC2 ({var_exp[1]:.1f}%)')
plt.legend()
plt.tight_layout()
plt.show()4.4 如何解读 PCA 图？
光会画图还不行，得会读图。我总结了三种典型情况：
| 图样 | 解读 | 行动 |
|---|---|---|
| 样本按批次完全分离 | 强批次效应，批次差异 > 生物学变异 | 必须进行批次校正 |
| 样本按批次部分重叠但有偏移 | 中度批次效应，部分批次间存在差异 | 可以考虑校正或保留为协变量 |
| 样本混合均匀，批次间无区分 | 无显著批次效应 | 放心往下做差异分析 |
你想想看，如果 PC1 和 PC2 加起来解释了 60% 以上的方差，而颜色正好按批次分开，那基本上可以确定批次效应是主要干扰因素。
4.5 用 PCA 量化批次效应的指标
光看散点图还是有点主观。我习惯再用一个量化指标——PC vs Batch 的 ANOVA。
具体做法：对每个主成分（比如 PC1 到 PC10），用批次做分组进行方差分析，计算 p 值。如果某个主成分的 p 值极小（比如 < 0.001），说明该主成分与批次显著相关。
# R 代码：计算每个主成分与批次的关联
pca_vars <- apply(pca_result$x[, 1:10], 2, function(pc) {
summary(aov(pc ~ batch))[[1]]$`Pr(>F)`[1]
})
names(pca_vars) <- paste0("PC", 1:10)
print(pca_vars)我在实际项目中，通常要求前 5 个主成分中至少有 3 个与批次显著相关，才确认批次效应需要校正。当然这个阈值你可以自己调。
4.6 本章知识体系图
4.7 常见误区与我的经验
- 只看 PC1 和 PC2： 批次效应有时藏在 PC3 或 PC4 里。我建议至少看前 5 个主成分的散点图矩阵。
- 过度解读小样本： 如果每组只有 2-3 个样本，PCA 分离可能是偶然，要结合统计检验。
- 忽略 outliers： 单个异常样本可能拉偏 PCA 方向。我曾经遇到一个污染样本，PC1 方差 80% 都是它贡献的，赶紧剔除重画。