做微生物组分析的朋友十有八九都被PCoA图刷过屏。群落样本在二维平面上一颗颗散开不同分组用不同颜色圈起来哪个处理组和对照组分得开、哪个样本偏离了队伍全都一目了然。这张看着简单的散点图背后其实是β多样性分析里最核心的一步先用距离算法量化样本间的群落差异再用主坐标分析把高维距离信息压缩到可视化的低维空间。这篇文章我就从PCoA的原理讲起把数据准备、距离矩阵选择、R语言实操、ggplot2出图以及环境因子拟合整条流程完整走一遍最后再列几个我实际跑数据时踩过的坑。无论你是刚接触16S扩增子数据分析的新手还是已经能跑通流程但想把PCoA图做得更专业的老手这篇内容应该都能给你一些参考。1. 先搞清楚PCoA到底在算什么1.1 β多样性样本之间的“距离感”要说PCoA得先从β多样性说起。微生物组研究里α多样性描述的是一个样本内部的物种丰富度和均匀度比如Shannon指数、Chao1指数而β多样性描述的是样本与样本之间的群落组成差异。想象你手里有10个土壤样本的OTU表每个样本里有几千个OTU的丰度信息。你想知道这10个样本是不是按照“污染组vs对照组”自然分成两拨但OTU有几千个每个样本就是一个几千维的向量人眼根本没法直接看。β多样性的思路就是把这几千维的信息浓缩成一个数字——两个样本之间的距离或差异度。比如样本A和样本B的Bray-Curtis距离是0.35样本A和样本C的距离是0.78那就说明A和B的群落组成更接近。这一步的核心产出是距离矩阵一个n×n的方阵n是样本数矩阵里每个值代表对应两个样本之间的群落差异。拿到这个距离矩阵之后才算真正进入PCoA。1.2 PCoA与PCA的关系从线性到广义主坐标分析Principal Coordinates AnalysisPCoA和主成分分析Principal Component AnalysisPCA名字很像确实也有血缘关系但思路不太一样。PCA是对原始数据矩阵直接做特征值分解它隐含的前提是用欧氏距离来衡量样本差异找到方差最大的方向进行投影。而PCoA更灵活它先算好任意一种距离矩阵然后对这个距离矩阵做经典多维缩放Classical Multidimensional Scaling把样本点嵌入到一个新的低维坐标系里使得新坐标系中点与点之间的欧氏距离尽可能接近原始距离矩阵中的差异值。用大白话说PCA是直接对数据降维PCoA是对“样本间的距离”降维。所以PCoA的输入不是OTU表本身而是OTU表算出来的距离矩阵。这也是为什么PCoA能兼容Bray-Curtis、Jaccard、UniFrac等各种各样的距离算法而PCA只能老老实实用欧氏距离。还有一个容易混淆的方法是NMDS非度量多维尺度分析。PCoA追求的是“保留距离的数值大小”NMDS追求的是“保留样本间距离的排序关系”。如果你关注的是样本间差异的绝对大小选PCoA如果你更关注排序关系且数据中有大量非线性结构NMDS往往表现更好。实际论文里PCoA出现频率更高一个重要原因就是它能把差异的百分比解释度算出来标在坐标轴上读者更容易理解。1.3 为什么微生物组几乎默认用PCoA微生物组数据有个鲜明特点OTU表极度稀疏而且很多OTU的丰度分布呈长尾少数高丰度物种占了绝大部分序列大量稀有物种只出现在少数样本里。这种情况下直接用欧氏距离算样本差异结果很容易被高丰度物种主导低丰度物种携带的信息几乎被淹没。Bray-Curtis距离等非欧氏距离对丰度组成更鲁棒它基于物种丰度的变化计算差异对零值不像欧氏距离那么敏感因此更贴合微生物群落的生态学意义。而PCoA恰好能处理任意距离矩阵两者搭配就成了微生物组β多样性分析的主流组合。另外像UniFrac这类利用系统发育信息的距离也是先算距离矩阵再做PCoA。可以说“算距离矩阵 PCoA降维 可视化”在微生物组领域已经是事实标准。2. 分析前必须做好的数据准备2.1 从测序下机到OTU表/特征表PCoA分析的第一步输入是OTU表或者现在更流行的叫法是特征表Feature Table。无论你用QIIME2、mothur还是DADA2流程最后都会得到一张行是OTU/ASV、列是样本的丰度矩阵。有的流程还会同时输出代表序列、分类学注释结果以及每个样本的元数据信息。拿到特征表之后先别急着算距离我建议先花点时间做三件事。第一看一眼测序深度分布如果某些样本的序列数特别少比如只有几百条而其他样本有几万条那这个样本很可能就是后面PCoA图上的离群点。第二检查样本ID是否和元数据表里的行名完全一致这看起来是小事但我在实际分析中遇到过太多次因ID格式不一致导致合并时报错的情况。第三确认特征表里的数值是整数序列数还是相对丰度这直接决定了后续标准化方式。在这些基础检查做完之前就往下跑分析后面出了问题很难排查。磨刀不误砍柴工。2.2 数据过滤与标准化这一步决定了距离算得对不对特征表里的OTU很多只在一两个样本里出现一次这类稀有OTU是不是要过滤掉取决于你的研究问题。如果你关注的是稀有物种对群落差异的贡献那就保留如果只是想看核心群落结构可以按“在至少20%的样本中出现且总丰度大于某阈值”来过滤。我个人的习惯是先用一个比较宽松的阈值跑一遍再收紧阈值跑一遍看结论是否稳定。标准化问题更关键。用Bray-Curtis距离时通常需要先对样本做相对丰度转换也就是把每个样本的序列数归一化到0到1之间消除测序深度差异的影响。QIIME2里可以用qime diversity core-metrics-phylogenetic它会自动做稀疏化rarefaction如果你用Rvegan包里有个decostand(x, method total)可以按行做总和标准化。注意稀疏化和总和标准化是两种思路。稀疏化是随机抽平到同一深度会丢失一部分数据总和标准化只是缩放不改变数据本身。如果你的样本测序深度差异很大两种方式都跑一遍对比一下最稳妥。如果差异不大我建议优先用总和标准化信息损失更少。还有一种更严格的标准化方式是CLR中心对数比变换它适用于成分数据分析的场景比如你想比较不同样本的物种丰度相对变化。但在常规的Bray-Curtis距离PCoA流程中CLR不是必须的做不做看你的研究设计。2.3 距离矩阵的选择逻辑距离算法的选择是PCoA分析中最需要动脑子的地方没有绝对的对错关键看你想回答什么问题。Bray-Curtis是最常用的它基于物种丰度的变化同时考虑了物种是否存在和丰度高低对大多数微生物群落分析场景都适用。Jaccard只关心物种存在与否不考虑丰度适合你对“物种有无”而非“丰度变化”更感兴趣的场景。UniFrac则把系统发育信息引入距离计算weighted UniFrac考虑丰度unweighted UniFrac只看谱系是否相同。实践中unweighted UniFrac往往能拉开稀有物种驱动的差异而weighted UniFrac更容易被高丰度物种主导。如果你手头有系统发育树建议两个都算分别做PCoA从不同角度解释数据。选距离矩阵的一个实用建议多试几种如果结论方向一致那你的生物学发现就很稳健如果不同距离得到完全不同的分组模式这时候你反而应该高兴因为这意味着有值得深挖的故事——可能是样本测序深度问题也可能是稀有物种确实在驱动差异。3. PCoA分析与可视化的完整实操3.1 R语言路线vegan包一行算距离一行做PCoAR语言做PCoA最常用的包是vegan。假设你已经读入了一个名为otu_table的数据框行名是样本ID列名是OTU ID元数据存在metadata里其中分组列叫Group。第一步计算Bray-Curtis距离矩阵library(vegan) # 按样本总和标准化 otu_norm - decostand(otu_table, method total) # 计算Bray-Curtis距离 bc_dist - vegdist(otu_norm, method bray)第二步做PCoA。vegan里可以用cmdscale做经典PCoA也可以用wcmdscale做加权版本。cmdscale输出的是原始坐标点但为了后续画图方便我习惯把它转成数据框pcoa_res - cmdscale(bc_dist, k 5, eig TRUE)这里的k是你想保留的维度数一般取5到10就够了eig TRUE是为了拿到特征值用来计算各轴的解释度。如果你用的是phyloseq包流程更简洁library(phyloseq) # physeq 是你的phyloseq对象 pcoa_ps - ordinate(physeq, method PCoA, distance bray) plot_ordination(physeq, pcoa_ps, color Group)phyloseq的好处是它会自动根据你的phyloseq对象中的分组信息生成ggplot对象方便后续直接用ggplot2语法继续调整。3.2 用ggplot2画一张可直接投稿的PCoA图拿到PCoA坐标之后大部分人第一个问题就是怎么把图做得像论文里那样好看我的答案是直接放弃vegan自带的plot函数改用ggplot2。原因很简单ggplot2的分组映射、颜色设置、主题调整、置信椭圆等功能都是现成的出图效果和自定义空间远胜基础绘图函数。先整理坐标数据并计算各轴解释度# 提取样本坐标 pcoa_points - as.data.frame(pcoa_res$points) colnames(pcoa_points) - paste0(PCo, 1:ncol(pcoa_points)) pcoa_points$SampleID - rownames(pcoa_points) # 合并分组信息 plot_data - merge(pcoa_points, metadata, by SampleID) # 计算解释度 eig - pcoa_res$eig explain - round(eig / sum(eig[eig 0]) * 100, 1)这里有个易错点eig里可能包含负特征值尤其是使用Bray-Curtis这类非欧氏距离时可能出现几个轴的特征值为负。计算解释度时分母用“所有正特征值之和”是更常见的处理方式。然后画图library(ggplot2) p - ggplot(plot_data, aes(x PCo1, y PCo2, color Group)) geom_point(size 3, alpha 0.8) stat_ellipse(aes(fill Group), geom polygon, alpha 0.2, level 0.95) labs(x paste0(PCo1 (, explain[1], %)), y paste0(PCo2 (, explain[2], %))) theme_bw() theme(panel.grid element_blank(), legend.position right) ggsave(PCoA_Bray_Curtis.pdf, p, width 6, height 5)stat_ellipse默认画的是95%置信椭圆如果样本量太少或者组内离散度过大椭圆可能画不出来这时候可以把level调低到0.8或者用geom_polygon手工画凸包convex hull。出图格式优先PNG和PDF双份PNG用于快速查看PDF用于投稿排版。3.3 进阶如何用phyloseq快速完成整套流程如果你每次从头写vegan代码觉得繁琐可以封装一个小函数以后只输入特征表和元数据就能直接出图。我自己的做法是写了一个名为quick_pcoa的函数内部处理数据标准化、距离计算、PCoA、合并分组、ggplot出图以及自动保存PNG和PDF参数只需要传入特征表、元数据、分组列名、距离类型这几个关键信息。phyloseq的优势在于它把特征表、样本元数据、分类注释、系统发育树封装在一个对象里所有分析都在这个对象上操作不容易出现“ID对不上”“数据错位”这类低级错误。如果你经常做微生物组分析我建议把phyloseq的流程走熟至少要知道ordinate和plot_ordination这两个核心函数。不过也值得明白phyloseq底层调用的还是vegan的vegdist和cmdscale理解手动流程能让你在出问题时更容易排查。4. 可视化进阶让PCoA图会“说话”4.1 添加置信椭圆和样本连线分组趋势一眼看出一张只有散点的PCoA图读者往往要自己眯着眼睛找规律。加上置信椭圆或凸包以后分组趋势就非常直观了。置信椭圆适合样本量较多的组它表达的是“这个组在95%置信水平下的大致范围”凸包则是把所有样本点包起来适合样本量少或分布不规则的情况。R里画置信椭圆推荐stat_ellipse它基于多元t分布计算。画凸包需要先用grDevices::chull提取每组的外围点再用geom_polygon填充。有一个细节如果你的组之间有重叠填充色的透明度一定要低不然下面一组的点会被完全盖住。我一般把填充alpha设成0.1到0.2散点的alpha设成0.8左右保证信息层次分明。如果样本还包含时间序列或梯度变化比如同一地点不同月份采样可以画样本点之间的连线用箭头表示变化方向这样PCoA图就不只是静态快照而是能展示动态轨迹。4.2 环境因子拟合envfit投影与箭头PCoA只能告诉我们样本怎么分布不能直接告诉我们“是哪些环境因子驱动了这种分布”。要回答这个问题vegan的envfit函数非常常用。# 假设env_data是环境因子数据框行名和plot_data的SampleID一致 fit - envfit(pcoa_res$points, env_data, permutations 999)envfit会把每个环境因子作为向量投影到PCoA排序空间permutations用来检验因子的显著性。结果可以用plot(fit)直接绘图也可以手动提取箭头坐标后用ggplot2重新绘制这样能保证整体风格统一。需要注意envfit适合连续型环境变量比如pH、含水量、温度。如果你的环境因子是分类变量比如土壤类型可以改用ordihull或者factorfit它会计算分类变量的组中心并检验组间差异。实际画图时箭头长度代表环境因子对排序空间的贡献大小箭头方向代表该因子增大时样本的趋势方向。两个箭头夹角小说明这两个因子正相关夹角接近90度则基本不相关。这个信息在解读群落结构与环境因素的关系时非常有用。4.3 三维PCoA与交互式图表什么时候值得做二维PCoA图最多只能展示前两个主坐标轴但有时候前两轴的解释度加起来不够高比如只有30%多这种情况下第三轴甚至第四轴的信息就很值得展示。三维PCoA可以帮你看清样本在更多维度上的分离情况。R里画三维图可以用scatterplot3d或rgl。scatterplot3d出静态图适合投稿rgl可以交互旋转适合自己探查数据。如果你想把三维图放到网页或者报告里供别人交互查看可以用plotly包把PCoA坐标转成长表格后直接plot_ly出3D散点图。但我个人的建议是论文里优先用二维图因为印刷和审稿人阅读体验更友好。三维图更适合作为补充附图或者当你发现二维图上两个组完全重叠、但第三维能明显分开时才值得为这个“隐藏维度”专门做一张图。5. 常见问题与排查技巧实录5.1 散点挤成一团什么都看不出来这种情况我见了不下十次。原因通常出在数据标准化或者距离选择上。第一步检查一下是否做了相对丰度标准化如果直接拿原始序列数算Bray-Curtis测序深度高的样本会自然聚在一起导致所有样本挤成一团。第二步看看有没有极端离群样本有时候单个污染样本或者测序失败的样本会把整个坐标空间拉扁其他样本全部挤到角落。找到离群点后先确认它是不是技术错误导致的比如样本标签错误、DNA量极低、测序失败等。如果是剔除后重新分析如果不是这本身可能就是生物学信号。另一个不太容易想到的原因是OTU表里大量OTU全为0导致距离矩阵里有很多相同的距离值。这时候可以尝试过滤低丰度OTU后再算距离或者改用Jaccard距离压一压信息。5.2 解释度百分比怎么算为什么图上要标注坐标轴解释度的计算在R的cmdscale结果里并不是自动给你一个百分比你得自己从特征值里算。代码是eig - pcoa_res$eig explain - eig / sum(eig[eig 0]) * 100这里有几个细节容易踩坑。第一cmdscale默认只返回前k个特征值如果你设了eig TRUE它会返回所有非负特征值这才可以用来算解释度。第二当使用非欧氏距离时特征值里可能出现负数这就是上面代码里分母用sum(eig[eig 0])的原因。第三不同版本的R或者veganeig的结构可能略有差异务必先str(pcoa_res$eig)看一眼再写计算代码。论文里坐标轴标签一般写作“PCo1 (32.5%)”这样的格式百分比通常保留一位小数。如果你的前两轴解释度加起来只有20%多别慌微生物组数据本来就噪声大低解释度是常态。只要分组分离趋势清晰、PERMANOVA检验显著就仍然值得报告。5.3 分组颜色、标签重叠、图例处理等“丑图”问题PCoA图最常见的丑不是分析错了而是细节没调好。样本标签重叠是最令人头疼的问题之一。样本多的时候直接geom_text会把图糊成一片推荐用ggrepel包的geom_text_repel它会自动把标签推开避免重叠。颜色选择也很关键。R默认的ggplot2配色在屏幕上看着还行但打印出来对比度不够。建议使用RColorBrewer的Set1、Dark2等色板或者scale_color_manual手动指定色值。我常用的做法是用ggsci包的lancet或nejm配色色盲友好而且有期刊风格。图例的位置、字体大小、是否保留网格线这些细节看似不起眼但直接影响审稿人的第一印象。我自己出图的标准配置是去掉了灰色背景和网格线图例放到图内侧或右侧字体统一为Arial或Helvetica图片分辨率不低于300dpi。这些小调整花不了两分钟但图的专业感立刻提升一个档次。5.4 一个完整的排查流程参考最后总结一下我遇到PCoA图异常时的排查顺序供你参考先检查特征表是否标准化再看距离矩阵计算有没有报错或警告然后检查分组信息和样本ID是否对齐接着检查是否有离群样本最后再调整可视化细节。这四步走下来绝大多数PCoA图的问题都能定位。我个人的体会是PCoA分析并不难难的是每一步都理解自己在做什么。距离矩阵算的是样本间的“差异感”PCoA把这种差异感翻译成平面坐标可视化则是把这个翻译结果讲给别人听。只有三件事都做扎实了一张PCoA图才真正有说服力。最后再分享一个小技巧跑完PCoA以后顺手用adonis2也就是PERMANOVA检验一下分组差异是否显著。PCoA图负责让读者“看到”差异PERMANOVA负责用统计学告诉读者“差异可信”两个结果配套报告你的β多样性分析才算完整。
