微生物组β多样性PCoA分析全流程：从距离矩阵到出版级可视化

刚接触微生物组分析的人，十有八九会被β多样性这层窗户纸卡住——OTU表、距离矩阵、主坐标分析、置信椭圆，每个词好像都认识，连在一起就不知道到底在做什么。尤其是PCoA：明明样品在图上分开了，审稿人却问你“轴标签上的百分比怎么算的”；明明你用的是Bray-Curtis距离，方法学部分却写成了PCA。这篇文章我尽量把微生物组β-多样性中PCoA分析及可视化这条链路从头到尾讲透，包括数据准备、距离算法选择、主坐标计算、出版级绘图和常见坑位，让对着R文档发愁的人也能直接复现出能放进论文的图。

1. 为什么β多样性分析首选PCoA？先搞清楚它和PCA的区别

1.1 β多样性的本质：样本间“谁和谁更像”

β多样性是生态学里的核心概念，描述的是不同样点之间物种组成差异有多大。放在微生物组研究里，我们通常面对的是一张OTU/ASV丰度表：行是样本，列是微生物分类单元，单元格里是测序得到的序列数。β多样性要回答的问题就是：这些样本之间，谁的微生物群落更像，谁和谁差异巨大。

PCoA（主坐标分析）是展示β多样性最常用的降维方法之一。它的思路很朴素：先算出所有样本两两之间的距离，得到一个距离矩阵，再把这个矩阵投影到低维空间，让你能在二维平面里直观看到样本的聚集和分散。相比PCA，PCoA不直接处理原始丰度表，而是处理距离矩阵，因此可以灵活选择各种生态学距离，这也是它在微生物组分析里如此流行的根本原因。

很多初学者以为PCoA和PCA只是名字有点像，实际用起来也差不多。这个误解会在方法学描述、结果解释甚至审稿回复时都带来麻烦。简单说，PCA的输入是“样本×变量”的原始矩阵，PCoA的输入是“样本×样本”的距离矩阵。两者目标都包含降维和可视化，但底层逻辑不同，适用场景也不同。

1.2 PCA的欧氏距离局限与PCoA的解决思路

PCA本质上是基于欧氏距离的。它把原始变量线性组合成新的主成分，让投影后方差最大化，所以非常依赖数据的数值特征。可是微生物组OTU表有个很让人头疼的特点：高度稀疏、大量零值。绝大多数微生物分类单元只在少数样本里出现，直接拿这样的表去做PCA，那些零值会主导方差计算，结果往往被几个丰度特别高的物种牵着走，群落之间真正有意义的差异反而被掩盖。

PCoA则绕开了这个问题。它先从样本间距离矩阵出发，用Gower中心化把距离信息转成可以特征分解的形式，再通过特征值分解得到每个样本的主坐标。这一步的关键是，它既不要求原始数据满足正态性，也不要求数据是连续型，只要你能算出任意两个样本之间的距离，哪怕是非欧氏距离，PCoA都能帮你找到低维坐标。

所以PCoA并不是“比PCA更好”，而是“更灵活”。如果你的数据本身满足欧氏距离的使用前提，比如一些经过中心化对数比变换的数据，PCA完全可以用，甚至更合适。但在标准的16S扩增子分析里，我们通常面对的是组成型、稀疏型数据，这时候PCoA搭配Bray-Curtis或UniFrac距离，能更真实地反映群落差异。

1.3 常见距离/差异度算法：Bray-Curtis、Jaccard、UniFrac怎么选

距离矩阵是PCoA的输入，选什么距离直接决定PCoA图的形状和生物学解释。我见过太多论文，图画得挺漂亮，结果一看方法学，距离算法没写，或者写错了。这里整理一下最常用的三种：

Bray-Curtis是目前扩增子研究里最主流的选择。它直接使用丰度信息，对样本间丰度比例的差异非常敏感，而且不要求物种有无的二元判断。它的计算公式看起来很复杂，但核心思想就是“共享物种的丰度占比越高，样本越相似”。

Jaccard距离只关心物种出现与否，不关心丰度。如果你的研究重点是“这个环境下有没有某种菌”，而不是“某种菌丰度高不高”，那Jaccard比Bray-Curtis更合适。缺点是它对测序深度非常敏感，低深度样本容易丢失稀有物种，从而高估样本间差异。

UniFrac距离则额外引入了系统发育树，把物种之间的进化关系也纳入距离计算。未加权UniFrac可以捕捉到“哪些谱系是某个环境特有的”，加权UniFrac则更关注“高丰度谱系的丰度变化”。如果你的分析对象是系统发育结构差异明显的群落，UniFrac比前两种更有解释力。但要注意，用UniFrac需要保证OTU/ASV的代表序列能正确构建/匹配到系统发育树，这一步经常成为处理流程里的瓶颈。

选距离没有绝对标准，但有一条底线：研究问题决定距离算法。想谈群落组成变化，选Bray-Curtis；想谈物种有无差异，选Jaccard；想谈系统发育层面的生态分化，选UniFrac。选完之后，PCoA图、PERMANOVA、ANOSIM都要用同一个距离矩阵，否则结果串不起来。

2. 从OTU表到距离矩阵：数据预处理和距离计算全流程

2.1 输入数据的标准结构：OTU表、样本元数据、系统发育树

拿到上游分析产生的OTU/ASV丰度表后，先别急着算距离。你得确认数据格式是规范的。标准输入包括三大部分：

• OTU/ASV丰度表：行名是样本ID，列名是OTU/ASV ID，内容为序列数或相对丰度。行和列方向千万别搞反，我见过不少新手直接把Qiime2导出的feature-table转置了再读，结果后续所有分析全错。
• 样本元数据表：至少包含一列样本ID和至少一列分组信息，比如疾病组/对照组、时间点、处理方式。这个表行名最好与OTU表一致，顺序无所谓，但名称必须完全匹配。
• 系统发育树（可选）：如果打算用UniFrac距离，需要一棵包含所有OTU/ASV tip的新ick树。通常来自DADA2或deblur上游流程，构建树时要注意树尖ID必须与丰度表列名一致。

数据读入R后，我建议先执行几项基础检查：确认维度、检查是否有NA、确认行名是否重复。R里rownames(otu)不能有重复；Python的DataFrame里index也要唯一。别以为这些基础检查浪费时间，大量PCoA图异常都源于数据匹配错误。

2.2 标准化与抽平：什么时候该用哪种策略

测序深度不同是微生物组数据绕不开的问题。一个样本测了5万条序列，另一个只测了8000条，如果不做处理，直接拿原始count算Bray-Curtis距离，测序深度差异会伪装成群落差异，这在PCoA图上经常表现为样本沿着第一轴按测序深度分开，而不是按分组分开。

常用的处理策略有两种。第一种是抽平（rarefaction），让所有样本的测序深度统一到一个相同数值。R中可以用phyloseq::rarefy_even_depth，它会随机抽样一定数量的序列，相当于“大家伙都按同一个量来参会”。抽平的好处是简单、直观、被审稿人广泛接受；坏处是丢弃了部分数据，而且需要设置随机种子才能保证结果可重复。建议在做正式分析前固定一个种子，例如set.seed(20240617)。

第二种是相对丰度归一化，即每个样本的每个OTU序列数除以该样本总序列数，乘以100得到百分比。这样保留了所有数据，但组成数据的“闭合效应”会在后续距离计算里产生伪相关，所以使用相对丰度时通常还要配合适当的变换，例如平方根变换或中心化对数比变换。

到底选哪种？我的建议是：如果样本量足够，且测序深度总体不太悬殊，抽平是最省心的选择；如果样本量很小，或者已经做了严格的批次平衡，相对丰度加变换也完全可以。重要的是在论文方法部分明确写清楚，不要假装自己没处理过。

2.3 用R和Python计算β多样性距离矩阵

距离矩阵计算是PCoA前最后的步骤。R里最常用的是vegan::vegdist：

library(vegan) otu <- read.delim("otu_table.txt", row.names = 1, check.names = FALSE) # 如果已经抽平，直接用原始count; 如果没抽平，建议先归一化 otu_norm <- otu / rowSums(otu) * 100 # 计算Bray-Curtis距离 dist_bc <- vegdist(otu_norm, method = "bray")

Python则以scikit-bio为主：

import skbio.diversity.beta_diversity as beta_diversity import pandas as pd otu_df = pd.read_table("otu_table.txt", index_col=0) # 按样本行归一化 otu_norm = otu_df.div(otu_df.sum(axis=1), axis=0) * 100 bc_dm = beta_diversity("braycurtis", otu_norm.values, ids=otu_norm.index)

计算完距离矩阵，务必检查一下：矩阵是否对称，对角线是否接近0，是否有NaN。一个很常见的翻车点是OTU表里存在空行（某些样本所有OTU都是0）或空列（某些OTU在所有样本中都是0）。空样本在计算距离时会出现NaN，空列虽然不影响距离本身，但会干扰后续排序和统计。建议计算距离前先过滤掉零丰度OTU，以及检查每个样本的总丰度是否为正数。

vegdist默认会把输入表当作样本×物种矩阵，所以行名和列名一定要摆对。如果你是从Excel复制来的数据，很容易出现列名带多余空格或#，这些都需要清理干净。

3. PCoA坐标计算与坐标轴解释：轴标签上的百分比到底怎么来的

3.1 从距离矩阵到主坐标：Gower中心化与特征值分解

拿到距离矩阵之后，PCoA的计算过程可以拆成三步，虽然R或Python里已经封装好了，但理解原理能帮你少踩很多坑。

第一步，把距离矩阵D转换成一个新矩阵A，公式是A = -0.5 * D²，这里的D²指每个元素都平方。第二步，对A做Gower中心化，也就是令每个元素减去行均值、减去列均值、再加上总体均值，得到中心化矩阵G。第三步，对G做特征值分解，得到特征值和特征向量。每个样本的PCoA坐标就是特征向量乘以对应特征值的平方根。

R里最常用的基础函数是cmdscale：

pcoa_res <- cmdscale(dist_bc, k = 10, eig = TRUE) coordinates <- pcoa_res$points eigenvalues <- pcoa_res$eig

这里k=10表示取前10个主坐标用于后续可视化，eig=TRUE让我们能拿到特征值，用来计算解释度百分比。

如果使用ape::pcoa，会额外处理负特征值的问题。负特征值出现的原因是非欧氏距离矩阵无法在低维空间被完美表示，Bray-Curtis和UniFrac这类非欧氏距离都可能出现。ape::pcoa的correction参数可以选择修正方法，常见的有cailliez、lingoes等。如果你的数据负特征值很多，普通cmdscale可能给出很奇怪的坐标，此时换成ape::pcoa更稳妥。

3.2 解释度百分比的计算：别再写错轴标签了

PCoA图每个轴标签后面的百分比，是把该轴的特征值除以所有特征值之和得到的。这个值和PCA的“方差解释度”性质类似，表示这个轴在多大程度上反映原始距离矩阵中的信息。

R里可以这样手动计算：

all_eig <- pcoa_res$eig all_eig[all_eig < 0] <- 0 # 负特征值一般不计入解释度 explain <- all_eig / sum(all_eig) * 100 axis1_percent <- round(explain[1], 1) axis2_percent <- round(explain[2], 1)

然后你在ggplot绘图时，轴标签就可以写成：

xlab(paste0("PCo1 (", axis1_percent, "%)")) ylab(paste0("PCo2 (", axis2_percent, "%)"))

这是很多论文里被质疑的点。我只写“PC1”“PC2”而不写百分比，是常见的疏漏；但更严重的是有人直接把PCA的坐标都用了，轴标签还写“PCoA”，这种张冠李戴在审稿人眼里基本是硬伤。

解释度百分比还影响你对图的判断。如果PCo1+PCo2只有20%多，说明二维图只能展示整个距离结构的一小部分，样本在图上看起来聚合或分散，可能只是局部关系。面对低解释度，不要强行解读，可以考虑展示PCo3/PCo4，或者换距离算法重新评估。

3.3 如何正确解读PCoA散点图和置信椭圆

PCoA图上每个点代表一个样本，点与点之间的欧氏距离近似反映它们在目标距离矩阵中的真实距离。注意我说的是“近似”，因为降维本身会丢失信息，所以当两个点在二维图里挨得很近，它们在原距离矩阵里往往也比较接近；但平面上距离较远的点，不一定是真正的极端差异，可能只是被投影到了不同方向。

置信椭圆是PCoA图上最常见的叠加元素之一。它通常使用组内点的多元正态分布95%置信区间，描述的是“这个组样本均值的位置估计”，而不是“所有样本都落在里面”。R里ggplot2::stat_ellipse默认就是95%置信椭圆，S Size较小或组内离散度过大时，这个椭圆会显得特别大，甚至超出图边界，这时候要谨慎表述。

另外，椭圆圈住的范围不等于聚类。如果两组样本的椭圆大部分重叠，说明组间差异不显著；但如果样本量足够，组间距离差异仍然可能通过PERMANOVA检测出来。所以PCoA图要和统计检验配合使用，看图的同时必须报p值。

4. 从基础散点到出版级可视化：R/Python完整代码与参数调整

4.1 用ggplot2绘制PCoA散点图：颜色、形状、椭圆、标签一次到位

我平时最常用的PCoA可视化方式是ggplot2，因为对分组、主题和输出格式的控制非常灵活。假设我们已经从cmdscale提取好了坐标，并且合并了分组信息：

library(ggplot2) library(vegan) # 假设coordinates是样本坐标矩阵，metadata包含sampleID和group pcoa_df <- data.frame(coordinates[, 1:2]) colnames(pcoa_df) <- c("PCo1", "PCo2") pcoa_df$sampleID <- rownames(pcoa_df) pcoa_df <- merge(pcoa_df, metadata, by = "sampleID") p <- ggplot(pcoa_df, aes(x = PCo1, y = PCo2, color = group, fill = group)) + geom_point(size = 3, alpha = 0.8) + stat_ellipse(geom = "polygon", level = 0.95, alpha = 0.2, linewidth = 0.5) + scale_color_manual(values = c("#0072B5", "#BC3C29", "#20854E")) + scale_fill_manual(values = c("#0072B5", "#BC3C29", "#20854E")) + labs(x = paste0("PCo1 (", axis1_percent, "%)"), y = paste0("PCo2 (", axis2_percent, "%)")) + theme_classic(base_size = 14) + theme(legend.position = "right", legend.title = element_blank())

关于图中点的透明度，样本量少(少于20)时alpha可以设成1，样本量多时设成0.6-0.8，避免重叠点遮挡信息。点的形状也可以按批次或另一个分组设置：shape = batch，但形状不宜超过4种，否则图例非常难读。

如果你的数据每组样本特别少（比如每组只有3个），置信椭圆会非常不稳定，这时候可以改用“凸包”geom_polygon配合chull函数绘制组内样本凸包，或者干脆不画椭圆只画散点。我在给别人审稿时见过不少每组2个样本还画椭圆的，这基本等于把“样本量不足”写在了脸上，属于减分项。

4.2 用Python matplotlib绘制PCoA图

Python生态下，我推荐用scikit-bio计算距离矩阵，再用matplotlib/seaborn绘图。完整代码可以这样写：

import matplotlib.pyplot as plt import numpy as np import pandas as pd from skbio.stats.ordination import pcoa # distances已经是skbio DistanceMatrix pcoa_res = pcoa(distances) coord = pcoa_res.samples[["PCo1", "PCo2"]].copy() coord["SampleID"] = coord.index coord = coord.merge(metadata, left_index=True, right_index=True) explained = pcoa_res.proportion_explained xlab = f"PCo1 ({explained.iloc[0]*100:.1f}%)" ylab = f"PCo2 ({explained.iloc[1]*100:.1f}%)" fig, ax = plt.subplots(figsize=(6, 4)) for group, color in zip(["A", "B", "C"], ["#0072B5", "#BC3C29", "#20854E"]): sub = coord[coord["group"] == group] ax.scatter(sub["PCo1"], sub["PCo2"], label=group, color=color, s=60, alpha=0.8) ax.set_xlabel(xlab) ax.set_ylabel(ylab) ax.legend() plt.tight_layout() plt.savefig("pcoa_python.pdf", dpi=300)

相比R，Python生态在交互式可视化方面更有优势。如果你想把PCoA图做成可以鼠标悬浮查看样本名的HTML交互图，可以无缝切到plotly.express.scatter，在hover_data里填入样本ID、分组、测序深度等列，导出的HTML文件作为论文补充材料非常方便。

4.3 导出高清图与组合排版：让Figure达到期刊要求

出版级PCoA图对分辨率、字体和图例排版都有要求。R里用ggsave导出时，我习惯同时输出PDF矢量版和300dpi的TIFF版：

ggsave("pcoa.pdf", p, width = 6, height = 4.5, dpi = 300) ggsave("pcoa.tiff", p, width = 6, height = 4.5, dpi = 300, compression = "lzw")

如果你想把PCoA图和PERMANOVA结果箱线图，或者α多样性图拼在一起，推荐用patchwork包：

library(patchwork) combined <- p + boxplot_plot + plot_layout(ncol = 2, widths = c(2, 1)) ggsave("combined_figure.pdf", combined, width = 9, height = 4.5, dpi = 300)

需要注意，有些期刊对图内字体要求统一为Helvetica或Arial，所以主题设置里可以加一句+ theme(text = element_text(family = "Arial"))。还有，图例标题如果不需要就删除，避免英文组名旁边出现“group”这种多余标签。

5. 组间差异检验：PERMANOVA和ANOSIM如何配合PCoA图

5.1 PERMANOVA的假设、使用场景和R实现

PCoA图给人视觉上的“分开了”，这只是第一步。科学结论必须有统计检验支撑。目前微生物组领域最主流的组间β多样性差异检验是PERMANOVA（非参数多元方差分析），也叫Adonis。它不依赖数据正态分布，而是直接对距离矩阵进行置换检验。零假设是“不同组的样本在距离矩阵中的质心位置没有差异”。

R里最标准的实现是vegan::adonis2：

set.seed(20240617) permanova <- adonis2(dist_bc ~ group, data = metadata, permutations = 999) print(permanova)

结果里最重要的两个指标是R²和p-value。R²表示分组因素能解释的距离变异比例，数值越大说明分组差异越明显；p值表示这种差异是否显著。关于permutations数量，期刊一般要求至少999，我通常设置9999以保证结果稳定，尤其当p值接近0.05的时候。

使用PERMANOVA要特别注意两点。第一，它对组间离散度（方差）的差异敏感。如果一组的样本特别分散，另一组特别紧密，PERMANOVA可能因为离散度差异而给出显著p值，而不是真正的质心差异。所以最好同时做vegan::betadisper检验：

disper <- betadisper(dist_bc, metadata$group) permutest(disper)

如果betadisper也显著，说明组内离差确实不同，PERMANOVA的结果要谨慎解读。第二，PERMANOVA对不平衡设计也很敏感，组间样本量差异过大时，置换检验的功效会下降，尽量保持实验设计均衡。

5.2 ANOSIM与PERMANOVA如何配合使用

ANOSIM（相似性分析）也是一种基于距离矩阵的置换检验，但它比较的是秩，而不是原始距离。ANOSIM的核心输出是R值：R接近1，说明组间距离显著大于组内距离；R接近0，说明组间差异不大；R为负，则组内差异反而更大。R代码是vegan::anosim：

set.seed(20240617) anosim_res <- anosim(dist_bc, metadata$group, permutations = 999) summary(anosim_res)

PERMANOVA和ANOSIM经常一起使用，但两者的侧重点不同。PERMANOVA对质心位置差异更敏感，ANOSIM则对组间秩差异更敏感。我的习惯是：如果两者结论一致，那就放心了；如果PERMANOVA显著但ANOSIM不显著，很可能是样本量或离散度问题，需要进一步检查。下表是我在实际分析中常用的选择逻辑：

5.3 如何把统计结果标到PCoA图上

统计结果标在图上，可以极大提升信息密度。最常见的做法是在图的右上角用annotate加上一行小字：

p_value <- permanova$`Pr(>F)`[1] p_label <- paste0("PERMANOVA: p = ", format(p_value, scientific = TRUE, digits = 2)) p <- p + annotate("text", x = max(pcoa_df$PCo1) * 0.8, y = max(pcoa_df$PCo2) * 0.95, label = p_label, size = 4, hjust = 0)

如果p值特别小，比如小于0.001，我倾向于直接用科学计数法显示，例如“p = 2.5e-04”。注意不要让文字被图例遮住，必要时调整legend.position到图的底部。

除了把p值标在PCoA图上，我还会额外画一张“组间距离箱线图”：对距离矩阵按分组组合拆开，计算同一组内样本两两之间的Bray-Curtis距离，画出箱线图，这样能直观展示组内离散度。这个图通常和PCoA图作为Figure的一部分拼在一起，审稿人看了会省力很多。

6. 实战排雷：五个让PCoA结果“看起来不对”的常见原因

6.1 距离矩阵与后续分析不匹配，统计结果张冠李戴

我在实际项目里遇到过这样的情况：PCoA图用UniFrac距离画的，样本分得很开，看起来非常漂亮；但跑PERMANOVA时，同事偷懒用了之前算好的Bray-Curtis距离矩阵，结果p值也不错。虽然两个距离矩阵都有效，但论文里的PCoA和PERMANOVA对应不同距离，这在统计上是不自洽的。要解决很简单，设置好同一个距离矩阵对象，PCoA和所有统计检验都用它。

还有一个隐藏的不匹配是：绘图时用了归一化后的相对丰度矩阵，而距离矩阵用了抽平前的原始count矩阵。如果两者存在本质差异，PCoA图上的坐标就会和PERMANOVA检验的数据基础不一致。因此每一步都要记录清楚，避免“数据血统”混乱。

6.2 抽平后OTU表仍有空行和零方差样本

抽平不是万能的。尤其是当某个样本的测序深度非常低，比如只有几百条序列，抽平到几千条时它可能被抽到只剩下极少数OTU，甚至所有OTU都变成0。这样的样本在计算距离矩阵时会得到NaN或全零向量，PCoA图上可能被映射到原点，还可能扭曲整个坐标空间。

所以抽平后一定要重新过滤：去掉总丰度为0的样本，去掉在所有样本中丰度都为0的OTU。还可以检查一下每个样本的测序深度是否足够。我一般会设定一个底线：如果最低样本深度低于总体中位数的1/10，那么这个样本要么重测，要么在分析前列为可疑样本。

6.3 样本量过少却硬画置信椭圆

置信椭圆默认假设每组样本来自多元正态分布，并基于组内的方差-协方差矩阵估计。当每组只有3到4个样本时，协方差矩阵估计极不稳定，椭圆形状可能非常夸张，甚至画到图的外面去。这种情况下，要么不画椭圆，要么改用凸包（convex hull）表示样本范围，方法学里写成“Convex hulls show the distribution of samples within each group”。

如果每组样本数已经少到只有2个，那就连凸包也别画了，直接只显示散点，靠点位的聚集程度配合PERMANOVA p值来传达信息。强行加椭圆只会让读者高估结果的稳定性。

6.4 分组颜色与形状使用不当，图的区分度反而下降

颜色选择直接影响可读性。我看过太多图，三个组用了红、绿、蓝，其中红色和绿色对红绿色盲读者来说几乎无法分辨。建议优先选用色盲友好的配色方案，例如R里scale_color_brewer(palette = "Dark2")，或者手动设定#0072B5、#BC3C29、#20854E这种经过验证的颜色。

如果同图还要区分第二个分组，比如不同时间点，可以用形状。形状不宜超过4种，而且图例要确保每个形状/颜色组合都有足够大的区分度。点的大小也需要考虑，我一般设置size = 3，样本多时缩小到size = 1.5并用透明度避免重叠。

6.5 忽视批次效应和置换检验种子

微生物组测序常分多批次完成，肉眼没看到的批次差异可能混入β多样性结果。一个快速检查方法是把批次/测序run信息也放到PCoA图上，用形状或颜色标注。如果样本明显按批次聚类而不是按实验分组聚类，就要考虑是否需要做批次校正。这时候可以做PERMANOVA把批次作为第二个解释变量，看批次因素的R²是否显著。

最后，置换检验的随机性经常被忽略。同一个数据和代码，每次跑adonis2，如果没有设置种子，p值会有一点点浮动。为了让结果可复现，务必在脚本开头设置set.seed()，并在方法学中写明随机种子，这一点很多新手都会漏。

说回到我自己的项目习惯：每次完成β多样性分析，我都会用sessionInfo()固定R版本和所有包版本，把距离矩阵、坐标、统计结果全部导出成CSV存档。这样即便几个月后审稿人要求重新作图，我也能完全复现当时的版本。PCoA本身不是难点，难点在于让每一步都建立在清晰的数据和合理的参数之上。把这篇文章里提到的细节理清，你的β多样性分析至少能少走一大半弯路。