4. 主坐标分析实战：使用R语言计算PCoA
好了，前几章我们把概念和数学原理捋了一遍。现在，是时候动手了。这一章，我们进入真正的实战环节——用R语言把PCoA算出来。我个人习惯把这一步叫做“灵魂拷问”，因为数据是好是坏，一算便知。
为什么要单独用一章讲计算？说白了，你如果连数据格式都没准备好，后面的可视化全是白搭。我在项目中遇到过太多次，学员兴致勃勃地导入数据，结果报错一大片。嗯，今天我们就来把这个坑填平。
4.1 准备环境：安装与加载必要的R包
做PCoA，你不需要装一堆乱七八糟的包。核心就几个。我建议你把它们一次装好，省得后面还要反复折腾。
# 核心包：vegan 提供 adonis2 和 vegdist
# ape 提供 pcoa 函数
# tidyverse 用于数据处理和 ggplot2 绘图
# phyloseq 如果你处理的是微生物组数据，会很方便
install.packages(c("vegan", "ape", "tidyverse"))
# 如果你想用更现代的可视化，也可以装
install.packages("ggrepel") # 用于文字标签防重叠options(repos = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/") 设置一下。
4.2 数据加载与格式检查
PCoA吃进去的是距离矩阵，但我们通常手头有的是物种丰度表或者基因表达矩阵。所以第一步是把原始数据变成距离。
先看看数据长什么样。我会用一个模拟的OTU表格做演示，你自己换成真实数据就好。
# 模拟一个示例数据：10个样本，100个OTU
set.seed(42)
otu_table <- matrix(
round(runif(1000, 0, 50)),
nrow = 10,
ncol = 100,
dimnames = list(
paste0("Sample", 1:10),
paste0("OTU", 1:100)
)
)
# 查看数据的前几行和维度
head(otu_table[, 1:6])
dim(otu_table)dim() 和 head() 看清结构。
4.3 计算距离矩阵
这是最核心的一步。距离度量的选择，会直接影响你后续看到的PCoA图形。不同的距离，反应的数据特性不一样。下面我列几种最常见的。
| 距离类型 | 适用场景 | R函数 |
|---|---|---|
| Bray-Curtis | 生态学、微生物组（丰度敏感） | vegdist(method = "bray") |
| Jaccard | 存在/缺失数据（不在乎丰度） | vegdist(method = "jaccard") |
| Euclidean | 基因表达、代谢组（标准化数据） | dist(method = "euclidean") |
| Aitchison | 成分数据（如相对丰度） | 需要 compositions 包 |
我个人习惯，处理微生物组数据时，首选Bray-Curtis。它既能反映物种有无，又能体现丰度差异。你想想看，如果一个物种在样本A里很多，在B里很少，Bray-Curtis就能把这个信息体现出来。
library(vegan)
# 计算Bray-Curtis距离矩阵
dist_matrix <- vegdist(otu_table, method = "bray")
# 查看距离矩阵结构
class(dist_matrix) # 应该是 "dist"
as.matrix(dist_matrix)[1:5, 1:5] # 转化成矩阵看一下4.4 执行PCoA计算
有了距离矩阵，接下来就是用ape::pcoa()或cmdscale()做主坐标分解。这两个函数的核心思想是一样的：把距离矩阵里的信息“压缩”到少数几个轴上。
library(ape)
# 方法1：使用ape包
pcoa_result <- pcoa(dist_matrix)
# 查看结果的结构
names(pcoa_result)
# [1] "values" "vectors" "trace" "correction"
# 提取主坐标（前两轴）
pcoa_axes <- pcoa_result$vectors[, 1:2]
colnames(pcoa_axes) <- c("PCoA1", "PCoA2")
head(pcoa_axes)# 方法2：使用基础R的cmdscale (效果类似)
pcoa_cmd <- cmdscale(dist_matrix, k = 10, eig = TRUE)
# 查看解释方差比例
eigenvalues <- pcoa_cmd$eig
variance_explained <- eigenvalues / sum(eigenvalues[eigenvalues > 0])
variance_explained[1:5]
为什么会有两种方法？ape::pcoa()在处理负特征值时更稳健——嗯，距离矩阵不一定是“欧几里得”的，有时会出现负特征值。ape包会自动做校正。而cmdscale比较简洁，但遇到负特征值可能会直接忽略。
4.5 解读PCoA结果
计算完毕后，你手上应该有两个关键信息：一是样本在新坐标系里的坐标（主坐标矩阵），二是每个轴解释了多少方差。别急着画图，先看看数字。
# 查看各轴解释的方差比例
# 从pcoa_result里提取
rel_eig <- pcoa_result$values$Relative_eig
head(rel_eig)
# 前两轴累计解释比例
cumsum(rel_eig)[1:2]一个经验值：如果前两轴能解释50%以上的总方差，那这张PCoA图就很有信息量了。如果不到20%，那说明数据噪音很大，或者你选择的距离度量不太合适。
我曾经在一个宏基因组项目里，Bray-Curtis PCoA的前两轴只解释了8%。当时我差点放弃，后来换成Jaccard距离，一下子炸到40%。所以，遇到解释度低，先换距离试试，别急着骂数据。
4.6 构建样本元数据与合并
PCoA坐标只是数字，没有分组信息的话，它就是一堆点。你需要把样本的分组信息（比如处理组、对照、时间点）和坐标合并在一起，才能用来做可视化。
# 假设你的元数据长这样
metadata <- data.frame(
SampleID = paste0("Sample", 1:10),
Group = rep(c("Control", "Treatment"), each = 5),
Time = rep(c("Day0", "Day7"), 5)
)
# 将PCoA坐标转化为数据框
pcoa_df <- as.data.frame(pcoa_axes)
pcoa_df$SampleID <- rownames(pcoa_axes)
# 合并
plot_data <- merge(pcoa_df, metadata, by = "SampleID")
head(plot_data)4.7 快速诊断：批次效应初判
你可能会问：PCoA都算出来了，怎么知道有没有批次效应？很简单。把样本点画出来，用批次信息去上色。如果不同批次的样本各自聚成一团，那问题就大了。
library(ggplot2)
# 假设元数据里有一个Batch列
# 临时加一列模拟
plot_data$Batch <- rep(c("Batch1", "Batch2"), each = 5)
p <- ggplot(plot_data, aes(x = PCoA1, y = PCoA2, color = Batch)) +
geom_point(size = 3, alpha = 0.8) +
stat_ellipse(aes(fill = Batch), geom = "polygon", alpha = 0.1) +
theme_minimal() +
labs(
x = paste0("PCoA1 (", round(rel_eig[1] * 100, 2), "%)"),
y = paste0("PCoA2 (", round(rel_eig[2] * 100, 2), "%)")
)
print(p)画完图，你一眼就能看出：如果Batch1和Batch2的点交叉在一起，那恭喜你，批次效应不大。如果两个批次之间隔着一条银河，那就要小心了，后续分析很可能被批次效应污染。
4.8 知识体系总结
这张图基本就是我做PCoA的完整路线。你按照这个顺序走，每一步检查一遍，基本不会出错。最重要的就是第二步和第五步——距离选不对，后面全白干；分组没对齐，图就是废的。