2. PCoA 数学原理：从距离矩阵到主坐标的推导过程与几何意义
好，咱们进入第二讲。在讲PCoA之前，我想先问一个问题：为什么主成分分析（PCA）在批次效应检测中有时不好使？
很简单。PCA要求输入是原始数据矩阵，它默认数据是连续的、线性的、具备欧氏空间结构。但生物信息学的距离度量常常是非欧的——比如Bray-Curtis距离、Jaccard距离、UniFrac距离。你没法直接把一个距离矩阵塞进PCA。
所以我们需要一个更普遍的方法：PCoA（主坐标分析）。它不管你的距离是怎么算的，只要有个距离矩阵，它就能给你找出一组“主坐标”，让样本在新的低维空间中尽量保持原来的距离关系。说白了，PCoA是把任意形式的距离翻译成坐标。
我个人习惯把PCoA看作“距离矩阵的翻译官”。它不关心你是用哪种距离尺子量的，只关心距离关系本身。下面，我带你一步步拆开这个翻译过程。
2.1 核心逻辑：从相异性到坐标
PCoA的输入是一个 \( n \times n \) 的距离矩阵 \( D \)，其中 \( D_{ij} \) 是样本 i 和 j 之间的距离。输出是 \( n \times k \) 的坐标矩阵 \( X \)，k通常是2或3，用来画图。
整个过程可以浓缩为三步：
- 平方与中心化：将距离矩阵转化为内积矩阵 G。
- 特征值分解：对 G 做特征值分解，得到特征值和特征向量。
- 降维取坐标：选前 k 个最大的特征值对应的特征向量，组成主坐标。
听起来不复杂，对吧？但每一步都有坑。我以前踩过一个坑，就是忘记对距离矩阵做“双中心化”直接做分解，结果画出来的图完全对不上样本的真实距离。嗯，这里要注意。
2.2 第一步：把距离变成内积
你可能会问：我们明明有距离，为什么还要转换成内积？
因为距离矩阵不能直接做特征值分解——它不满足“正定”或“半正定”性质。而内积矩阵可以。你想想看，做特征值分解天生就是对内积矩阵设计的。
转换过程如下：
1. 对原始距离矩阵 D 的每个元素取平方
D_sq[i][j] = D[i][j]²
2. 对 D_sq 做“双中心化”
设 J = I - (1/n) * 1*1ᵀ
则 G = -0.5 * J * D_sq * J
这个矩阵 G 就是我们要的内积矩阵。它的几何意义是：在某个高维的“嵌入空间”中，样本 i 和 j 的点积。
这里有一点很关键：双中心化操作保证了 G 的行和与列和均为零。为什么？因为只有这样，后续的特征向量才能代表偏离中心的坐标，而不是把整个点云整体平移。
2.3 第二步：特征值分解
得到内积矩阵 G 之后，剩下的工作就和PCA的后半段一模一样了。
对 G 做特征值分解：
G = V * Λ * Vᵀ
其中：
Λ = diag(λ₁, λ₂, ..., λₙ)，λ₁ ≥ λ₂ ≥ ... ≥ λₙ
V 的每一列是特征向量
这里注意，特征值有可能出现负数——这在PCA里不会发生，但PCoA会。为什么？因为我们的距离矩阵可能对应一个“非欧几里得”空间。比如Bray-Curtis距离在某些情况下无法被完美嵌入到欧氏空间。
负特征值的处理方式：
- 如果负特征值的绝对值很大——你选的距离度量可能不适合做PCoA
- 如果负特征值很小（接近零）——忽略它们，只取正的特征值
- 在实践中，如果前两三个特征值占总特征值和的80%以上，负特征值基本不影响结果
2.4 第三步：提取主坐标
我们只取前 k 个最大的特征值及其特征向量：
主坐标矩阵 X 的每一行是一个样本的 k 维坐标：
X = V_k * sqrt(Λ_k)
其中：
V_k 是前 k 个特征向量组成的 n×k 矩阵
Λ_k 是前 k 个特征值组成的 k×k 对角矩阵
这里的 \( \text{sqrt}(Λ_k) \) 很关键——它把内积重新缩放回距离单位。不要省略它，否则你的坐标尺度就是内积尺度而非距离尺度，画出来的散点图会看起来“变瘦”。
我个人习惯在提取主坐标后，再算一遍各样本在新空间的欧氏距离，与原始距离矩阵做一次对比校验。如果两张图的点间距离一致性大于0.95，我就对当前结果有信心了。
2.5 几何意义：PCoA到底做了什么？
从几何角度看，PCoA本质上是在做这么一件事：
- 隐式恢复：从距离矩阵出发，恢复出样本在“某个欧氏空间”中的坐标
- 最紧拟合：在这个恢复出的空间中，找到最能解释距离信息的 k 个方向
- 正交投影：把样本投影到这 k 个方向上，得到主坐标
你可以这么理解：假设原来每个样本有100个看不见的维度，你只知道它们之间的相互距离。PCoA帮你在保持距离结构的基础上，用2维或3维把这100维的“形态”大致勾画出来。这是一个非常强的降维能力，但它也有代价——只保留距离，丢失了原始变量的可解释性。
举个具体的例子：样本A和样本B的Bray-Curtis距离是0.3，样本B和C的距离是0.8。PCoA找到的坐标空间中，A和B的欧氏距离接近0.3，B和C接近0.8——但你不能说“PCoA的PC1轴代表丰度差异”，因为轴没有变量标签。
这就是PCoA和PCA的根本不同：PCA的轴是有数学含义的（方差最大方向），PCoA的轴只保证距离。
2.6 知识结构总览
下面这张SVG图把PCoA的完整流程串起来了。从距离矩阵到主坐标，每一步的输入输出和关键操作都标出来了。
2.7 PCoA与PCA的对比总结
我经常在项目里被问到：PCoA和PCA到底什么关系？能不能画完PCoA再画PCA？这里我用一个表格把它们的区别理清楚。
| 对比维度 | PCA | PCoA |
|---|---|---|
| 输入数据 | 原始特征矩阵 (n × p) | 相异性/距离矩阵 (n × n) |
| 数学基础 | 协方差矩阵的特征分解 | 双中心化距离矩阵的特征分解 |
| 距离类型限制 | 必须是欧氏距离 | 任意距离度量（Bray-Curtis, Jaccard, UniFrac等） |
| 输出主成分含义 | 变量方向的线性组合（可解释特征贡献） | 无量纲坐标，不直接对应原始变量 |
| 批次效应图特点 | 展示基于方差的分离，批次与主变量混杂时不易分离 | 展示基于距离的分离，批次差异明显时更直观 |
| 计算复杂度 | O(n·p²) 一般取决于特征数 p | O(n³) 因为要对 n×n 矩阵分解，当 n 很大（>10⁴）时较慢 |
从这张表你能看出来——PCoA是PCA的泛化版本。如果直接把欧氏距离矩阵丢进PCoA，得到的结果和PCA几乎一样（只是坐标旋转上有差别，但相对位置关系完全保留）。
2.8 什么时候用PCoA？
我个人的判断标准很简单：
- 如果你手头是原始丰度表，且批次效应主要来自技术变异——我会先用PCA快速预览
- 如果你手头只有距离矩阵（比如从Qiime导出的beta多样性结果）——只能用PCoA
- 如果你的数据存在高比例的双零（比如微生物组数据中很多样本和物种都是零）——用Bray-Curtis距离+PCoA效果远好于PCA
- 如果你想以图示方式检验批次效应，不关心主成分的生物学含义——PCoA更自由
在实践中，PCoA的图上如果同一个批次的样本聚成一团，而不同批次的样本明显分开，这就说明需要做批次校正了。至于校正方法，我们后面的章节会讲。