主成分分析（PCA）原理

讲 PCA 之前，我先问你一个问题：如果给你一个 100 维的数据集，你怎么看它的整体分布？直接画图？不现实。人眼只能看到 3 维，更别说 100 维了。PCA 就是干这个的——它能把高维数据压缩到我们能看懂的低维空间，而且尽量保留最重要的结构。

我最早接触 PCA 是在做单细胞 RNA-seq 数据分析的时候。那时候一个样本就有 20000 个基因的表达量，你能想象 20000 维的空间长什么样吗？反正我想象不出来。PCA 帮我降到了 2 维，虽然丢了很多信息，但细胞亚群的分离情况一眼就能看出来。这就是它的价值。

协方差矩阵——数据集中关系的地图

想搞清楚数据在哪些方向上变化大，必须先知道变量之间的相关性。协方差矩阵就是干这个的。假设我们有 m 个样本，每个样本有 n 个特征（变量），数据用 X（m×n 矩阵）表示。协方差矩阵 C 是个 n×n 的对称矩阵，第 i 行第 j 列的元素是第 i 个特征和第 j 个特征的协方差。

公式很简单：

C = (1/(m-1)) * (X - μ)ᵀ (X - μ)
其中 μ 是每个特征的均值向量

为什么用 m-1 而不用 m？这是无偏估计，样本协方差的分母是 m-1。我自己写代码时经常忘了这一步，结果算出来的特征值总是偏小。嗯，细节决定成败。

协方差矩阵的对角线是每个特征的方差，非对角线反映特征间的线性关系。正值表示正相关，负值表示负相关，绝对值越大关系越强。

特征值分解——把矩阵拆开看

获得了协方差矩阵 C，下一步就是找到它的特征向量和特征值。数学上来说，对于矩阵 C，如果存在非零向量 v 和标量 λ 满足：

C · v = λ · v

那么 v 就是特征向量，λ 是对应的特征值。特征向量代表了数据变化的“方向”，特征值代表了在那个方向上“变化有多大”。说白了，特征值越大的方向，数据在这个方向上越分散，包含的信息越多。

我自己的理解方式：把 C 想象成一个橡皮筋网格，特征向量是它的拉伸方向，特征值是拉伸幅度。PCA 就是找到拉伸最快的那些方向。

具体到数值计算，我们可以用 Python 的 NumPy 来演示：

import numpy as np
# 假设 X 是已经标准化后的数据矩阵（m×n）
C = np.cov(X, rowvar=False)
eigvals, eigvecs = np.linalg.eigh(C)  # 对对称矩阵用 eigh 更快
# 特征值默认升序，降序排列
idx = np.argsort(eigvals)[::-1]
eigvals = eigvals[idx]
eigvecs = eigvecs[:, idx]

降维思想——从 n 维到 k 维

拿到特征值和特征向量后，降维就简单了：按特征值从大到小排序，取前 k 个特征向量组成投影矩阵 W（n×k），然后把原始数据 X 投影到新的 k 维空间：

Y = X · W

Y 就是降维后的数据，每行对应一个样本，每列是一个主成分。通常 k 取 2 或 3 用于可视化。但怎么选 k？看累计方差贡献率。计算每个特征值占总特征值之和的比例，累积到 80% 或 90% 时对应的 k 就是推荐值。

total = eigvals.sum()
ratio = eigvals / total
cumsum = np.cumsum(ratio)  # 累计贡献率
# 例如：取前 5 个主成分贡献率 0.85，则 k=5

一张图看懂 PCA 全流程

下面是 PCA 从原始数据到低维投影的完整流程。我把核心步骤画成了框图，你对着看应该能秒懂。

从图上你能看出，整个流程是流水线式的：原始数据 → 标准化 → 协方差矩阵 → 特征分解 → 排序选向量 → 投影 → 得到降维结果。其中虚线表示你可以用累计贡献率来指导 k 的选择。

回到批次效应场景

在单细胞或批量转录组数据中，批次效应往往表现为不同批次的数据在低维空间里形成“团块分离”。用 PCA 降维到 2D 后，如果不同颜色的批次点明显分成几撮，那就说明批次效应存在。但 PCA 只能告诉你“有没有”，不能帮你“消除”。它更像一个诊断工具。

我个人习惯的做法是：先跑一次 PCA 看看总体分布，如果批次混在一起（理想情况），那太好了，直接继续分析。如果各批次错位，那就需要用后面的 ComBat、Harmony 或 UMAP 局部对齐来处理。

记住，PCA 的强项在于它的数学透明性和可解释性。特征向量就是主成分加载量，你可以看哪些原始变量对主成分贡献最大——这在某些场景下能帮你找到批次效应的源头（比如某个批次特异的 PCR 扩增引物）。