刚接触微生物组分析的人十有八九会被β多样性这层窗户纸卡住——OTU表、距离矩阵、主坐标分析、置信椭圆每个词好像都认识连在一起就不知道到底在做什么。尤其是PCoA明明样品在图上分开了审稿人却问你“轴标签上的百分比怎么算的”明明你用的是Bray-Curtis距离方法学部分却写成了PCA。这篇文章我尽量把微生物组β-多样性中PCoA分析及可视化这条链路从头到尾讲透包括数据准备、距离算法选择、主坐标计算、出版级绘图和常见坑位让对着R文档发愁的人也能直接复现出能放进论文的图。1. 为什么β多样性分析首选PCoA先搞清楚它和PCA的区别1.1 β多样性的本质样本间“谁和谁更像”β多样性是生态学里的核心概念描述的是不同样点之间物种组成差异有多大。放在微生物组研究里我们通常面对的是一张OTU/ASV丰度表行是样本列是微生物分类单元单元格里是测序得到的序列数。β多样性要回答的问题就是这些样本之间谁的微生物群落更像谁和谁差异巨大。PCoA主坐标分析是展示β多样性最常用的降维方法之一。它的思路很朴素先算出所有样本两两之间的距离得到一个距离矩阵再把这个矩阵投影到低维空间让你能在二维平面里直观看到样本的聚集和分散。相比PCAPCoA不直接处理原始丰度表而是处理距离矩阵因此可以灵活选择各种生态学距离这也是它在微生物组分析里如此流行的根本原因。很多初学者以为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是否衡量样本间物种丰度组成差异16S扩增子最常用Jaccard否只看有无否关注物种存在/缺失的差异对稀有物种敏感未加权UniFrac否只看有无是结合系统发育树对稀有谱系敏感加权UniFrac是是结合系统发育和丰度对丰度梯度敏感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::vegdistlibrary(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_col0) # 按样本行归一化 otu_norm otu_df.div(otu_df.sum(axis1), axis0) * 100 bc_dm beta_diversity(braycurtis, otu_norm.values, idsotu_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里最常用的基础函数是cmdscalepcoa_res - cmdscale(dist_bc, k 10, eig TRUE) coordinates - pcoa_res$points eigenvalues - pcoa_res$eig这里k10表示取前10个主坐标用于后续可视化eigTRUE让我们能拿到特征值用来计算解释度百分比。如果使用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”这种张冠李戴在审稿人眼里基本是硬伤。解释度百分比还影响你对图的判断。如果PCo1PCo2只有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_indexTrue, right_indexTrue) explained pcoa_res.proportion_explained xlab fPCo1 ({explained.iloc[0]*100:.1f}%) ylab fPCo2 ({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], labelgroup, colorcolor, s60, alpha0.8) ax.set_xlabel(xlab) ax.set_ylabel(ylab) ax.legend() plt.tight_layout() plt.savefig(pcoa_python.pdf, dpi300)相比RPython生态在交互式可视化方面更有优势。如果你想把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::adonis2set.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::anosimset.seed(20240617) anosim_res - anosim(dist_bc, metadata$group, permutations 999) summary(anosim_res)PERMANOVA和ANOSIM经常一起使用但两者的侧重点不同。PERMANOVA对质心位置差异更敏感ANOSIM则对组间秩差异更敏感。我的习惯是如果两者结论一致那就放心了如果PERMANOVA显著但ANOSIM不显著很可能是样本量或离散度问题需要进一步检查。下表是我在实际分析中常用的选择逻辑检验方法主要假设输入最怕什么推荐场景PERMANOVA组间质心位置不同距离矩阵 分组离散度异质性多数分组比较ANOSIM组间秩差异大于组内秩差异距离矩阵 分组组内样本量过少辅助验证PERMANOVAbetadisper组间离散度相同距离矩阵 分组不平衡设计作为PERMANOVA辅助检查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本身不是难点难点在于让每一步都建立在清晰的数据和合理的参数之上。把这篇文章里提到的细节理清你的β多样性分析至少能少走一大半弯路。