在微生物组研究中整合系统发育和功能数据

摘要

动机

经常分析微生物组功能数据，以确定微生物功能（例如基因）和感兴趣的样本组之间的关联。然而，仅考虑功能就难以区分社区功能概况变化的不同可能解释。为了帮助解决这个问题，我们开发了POMS公司这是一个实现多个系统发育感知框架的包，可以更稳健地识别丰富的功能。

结果

关键贡献是一个扩展的平衡树工作流，它包含功能和分类信息，以识别在独立分类谱系的样本组中一致丰富的功能。我们的软件包还包括一个运行系统发育回归的工作流。基于模拟数据，我们证明，与常用工具相比，这些方法能够更准确地识别具有选择性优势的基因家族。我们还表明，POMS尤其可以识别现实世界宏基因组数据集中的丰富功能，这些功能是微生物群多个成员强选择的潜在目标。

可用性和实施

这些工作流在POMS公司R包位于https://github.com/gavinmdouglas/POMS网站.

补充信息

补充数据可在生物信息学在线。

1引言

微生物组测序已被应用于表征无数环境，通常根据微生物特征的相对丰度进行分析。这些特征可能包括分类群（微生物存在）和功能（它们编码的基因和途径）。虽然这两种数据类型都被用来进行有价值的观察，但它们通常是独立分析的。然而，需要将这些数据类型联系起来，才能对微生物组数据中观察到的变化做出连贯的解释(道格拉斯和兰吉尔，2021年;Manor和Borenstein，2017年). 例如，在一组样品中富集特定微生物途径（例如，疾病或独特环境）可以代表对具有该途径的多个独立分类群的选择，这将具有重大的生物学意义。相反，由于特定微生物分类单元中存在该途径，因此也可能出现这种富集，该分类单元因兴趣组中的其他各种原因而更加丰富，使得这种富集在生物学上不那么有趣。

已经开发了两种方法，部分解决了如何链接分类和功能数据类型的挑战。鱼墨西哥卷(Manor和Borenstein，2017年)是一种通过差异丰度测试确定功能转变的分类贡献者的工具。鱼墨西哥卷然而，这是一种事后工具，在确定了重要的微生物功能后应用，而理想情况下，功能和分类数据将在测试差异功能时进行整合，以更好地确定强富集候选。

使系统进化是一种在统计测试期间显式集成功能和分类数据的方法(布拉德利等。, 2018;Bradley和Pollard，2020年). 它根据编码每个基因家族的分类群的流行率或特异性来确定功能关联。这是使用系统发育线性模型进行的，该模型解释了由于共同进化历史而产生的共生分类群的遗传相似性。利用这一框架，确定了由特定门内的一组不同分类群所贡献的重要基因家族和途径。虽然这一方法是对过去分析方法的一个关键改进，但尚未被广泛采用。这至少部分是因为系统发育要求用户将分析限制在MIDAS数据库中的基因组(Nayfach公司等。, 2016).

整合分类和功能分析的复杂性，以及微生物组数据分析中的其他挑战，源于将标准统计方法应用于原始微生物组数据的困难，因为这些数据具有组成性质。幸运的是，人们对用于分析微生物组数据的改进合成方法越来越感兴趣。具体来说，最近有人提出分析微生物组特征相对丰度的比率（而不是特征本身的丰度）来解决组成问题(格洛尔等。, 2016;莫顿等。, 2019). 然而，目前尚不清楚应使用哪些特征来稳定计算这些比率。一个建议的解决方案是比较连接不同分类群的系统发育树中每个节点两侧的分类群比率（基于等距对数比率转换）(西尔弗曼等。, 2017). 这种通用方法现在通常称为分析平衡树(莫顿等。, 2017). 虽然这是一个统计上有效的微生物组数据分析框架，但如何解释分类比率的差异通常并不清楚。

在此，我们旨在通过在平衡树框架中测试功能富集来解决功能差异丰度分析所涉及的挑战和基于平衡树方法的有限解释性。具体而言，我们的方法，称为元基因组信号的系统发育组织（POMS），重点是识别编码给定功能的多个分类单元与样本组一致关联的情况。这种情况为以下假设提供了更多的支持：功能本身赋予相关样本群中的微生物一种选择性优势，而不是该功能恰好出现在该群中更丰富的少数类群中。除了POMS工作流，我们还重新实现了系统发育，这使得这种系统发育回归框架可以应用于新的基因组。我们发现，每个工作流程都在模拟和真实宏基因组数据中精确定位有趣的功能关联，并且它们具有互补的优势和劣势。这是一个有价值的概念证明，将功能丰富性集成到平衡树分析中可以提高其可解释性，并提供新的见解。然而，重要的是，由于POMS不能将分类广度有限的功能识别为高度丰富的功能，我们主要将其视为系统发育回归和其他更敏感方法的补充工具。

2材料和方法

2.1等距对数比

2.2多项式检验

POMS中使用多项式精确检验来识别一致富集函数（CEF），如第3节所述。该测试考虑了每个函数三类功能重要节点（FSN）的计数：（i）不与平衡重要节点（BSN）相交的FSN，（ii）与BSN相交的FSNs，其中功能富集位于样本组1中相对更丰富的分类群中，以及（iii）与BSN相交的FSN，其中功能富集位于样本组2中相对更丰富的分类群中。每个类别中FSN比例的零期望基于FSN和BSN之间的质量作用交互。换句话说，空期望对应于随机分配FSN和BSN的情况。预期还假设与BSN相交的FSN在样本组1和样本组2中较高的BSN之间平均分割。此测试是使用的xmulti函数实现的X标称值R包（使用版本1.0.4测试），使用默认的对数似然比统计来计算P（P）-价值

2.3系统发育回归

我们使用来自藻门(Ho and Ané，2014年)R包，它根据每个测试功能的拷贝数回归一个载体，同时考虑到提示的系统发育相似性。所有分析均使用默认的布朗运动模型。我们在POMS公司用于运行此过程以及计算作为系统发育。后面的这些功能是根据系统发育v0.94码基(https://bitbucket.org/pbradz/系统发育)，它是在麻省理工学院的许可证下分发的。

2.4预处理功能拷贝数表

在运行所有分析之前，使用相同的过滤截止值过滤功能拷贝数表（[将宏基因组组装基因组（MAG）和其他分类群链接到基因组注释]）。首先，表格仅限于包括至少一个用于相关分析的样本中的分类群。然后，在少于五个分类群或小于0.1%的分类群中发现的任何功能都被删除。

2.5基于MAG的模拟

基于MAG的模拟基于一个大型人类荟萃分析数据集的704个对照样本(阿尔梅达等。, 2019). 这些样本都符合元数据标准，被标记为来自健康成年人而非抗生素。这些样本也来自总共至少40个样本的研究。随机抽取一个样本，创建平衡组。MAG丰度被视为先前计算的平均读取深度（即每个站点的平均映射读取数）。样本中覆盖宽度小于25%的MAG在该样本中的平均深度为0。然后，我们排除了所有未在任何样本中调用的MAG，剩下1595个MAG用于模拟。

这些模拟按照第3节所述进行。首先，将1000个重复数据集中的每个样本随机分为两组。然后，对于每个剖面，随机选择一个焦点基因家族，在一个组中，编码该基因家族的所有MAG的丰度以假计数1递增，这些丰度乘以1.5倍。这些模拟数据集被称为“焦点基因”剖面。我们还进行了平行模拟，其中随机分类群的相对丰度仅在一组中随机膨胀了相同的数量。重要的是，相同数量的分类群受到干扰，因为在每个匹配的焦点基因模拟剖面中受到影响。生成的模拟剖面称为“随机分类群”剖面。“基于分支的”剖面表示一组额外的模拟，这些模拟使用相同的选择压力，但仅适用于系统发育树中给定节点的分类群后代。测试树中693个节点中的每个节点都有一个基于分支的数据集，其中至少有五个基本分类单元（不包括根节点）。

我们对这些模拟数据集进行了两种类型的系统发育回归。首先，模型反应是分类群对第一个样本组的特异性。在第二种情况下，根据Wilcoxon检验（未修正P（P） < 0.05). 本研究中比较的其他差异丰度方法为Wilcoxon检验，采用通用单拷贝基因（USCG）丰度中位数进行归一化（使用MUSiCC工具使用的相同方法；Manor和Borenstein，2015年)，基于相对丰度的Wilcoxon检验，ALDEx2 v1.16.0(费尔南德斯等。, 2014)，DESeq2 v1.24.0(爱等。, 2014)和limma-voom v3.40.6(法学等。, 2014;里奇等。, 2015). 根据Benjamini–Hochberg截止值<0.05，对所有测试方法（包括POMS）确定了重要的基因家族。

2.6塔拉海洋数据集验证

组装的Tara Oceans宏基因组数据集取自一个预先存在的项目(德尔蒙等。, 2018). 海洋样品的环境和化学剖面以及MAG的相对丰度和注释取自出版的补充表格为这个项目。

GToTree v1.4.16(Lee和Ponty，2019年)基于共享的单拷贝基因构建系统发育树，并排除完整性低于60%、冗余度高于10%的MAG，从而保留642个MAG。如果覆盖范围大于1%，则称MAG存在于样本中。这是一个宽松的设置，是根据样本中MAG覆盖范围的经验分布选择的。

根据从KO到从KEGG网站下载的这些类别的KEGG映射，重建了高级功能（即KEGG模块和路径）(卡内希萨等。, 2016)2021年4月12日。使用PICRUS2公司脚本pathway_pipeline.py(道格拉斯等。2020)，利用MinPath（最小路径）(Ye和Doak，2009年)以及在中实现的算法人类2(弗兰佐萨等。, 2018)重建更高层次的功能丰度。对每个MAG分别进行这些重建。

由于环境数据（如盐度和营养盐浓度）是连续的，因此根据样品平衡和这些环境因素（未修正）之间的斯皮尔曼相关性，确定了显著的BSNP（P） < 0.05). 然后将这些结果离散化，以表明这些因素与每个BSN的组是正相关还是负相关。基于分类单元丰度是否显示出显著的Spearman相关性（未修正P（P） < 0.05）。还根据相对丰度本身与环境因素的函数计算了Spearman相关性。每个样本的分类丰度和每个分类群的功能丰度的交叉乘积用于产生社区范围的功能丰富度。

最后，Faith的系统发育多样性通过POMS或系统发育回归基于编码每个重要命中输出的分类群子集进行计算。这是使用皮坎特R包（v1.8.2；肯贝尔等。, 2010).

2.7病例对照枪宏基因组数据集验证

我们将某些验证分析集中在三个数据集上，这些数据集是对人类鸟枪宏基因组数据集进行的大型荟萃分析的一部分(阿尔梅达等。, 2019). 这些数据集是根据欧洲核苷酸档案中的数据集加入标识符定义的。在执行与基于MAG的模拟相同的预处理步骤后，我们使用先前生成的MAG、样本MAG丰度剖面和MAG系统发育树作为POMS的输入，但使用的是不同的样本子集。根据每个数据集病例组中每个分类单元的样本特异性进行系统发育回归。Faith的系统发育多样性计算和KEGG通路和模块的重建与Tara Oceans数据集相同。

2.8 POMS依赖性

POMS公司用R书写(R核心团队，2019年)并且依赖于以下R包（指出了本手稿中使用的版本，但不需要这些确切的版本）：猿版本5.3(天堂等。, 2004),平行v3.6.0，潘戈恩第2.5.5节(Schliep，2011年),藻门v2.6.4，字符串v1.4.0和X标称值v1.4。POMS公司v0.3.1用于所有分析。该方法的测试和开发是在运行Ubuntu v16.04.5的服务器上使用R v3.6.0和RStudio v1.2.5033进行的。

运行后，需要几个额外的R包来遵循当前的分析工作流POMS公司（同样，所示版本用于本文，但不是必需的版本）：ggtree（ggtree）(于，2020;于等。, 2017)v1.16.1，ggplot2第3.3.0版(威克姆，2016年),普利尔1.8.4版(威克姆，2011年)和重塑2第1.4.3节(威克姆，2007年). 所有显示的多面板图都是用奶牛场（v1.0.0）R包。

2.9代码可用性

POMS源代码位于：https://github.com/gavinmdouglas/POMS网站。本手稿中所有分析的代码可从以下网址获得：https://github.com/gavinmdouglas/POMS手册/.

3结果

3.1 POMS概述

POMS是一种用于分析微生物功能的平衡树框架(图1)包括基因家族和高级功能。关键输入表对应于样本的分类丰度和轴突功能丰度。还必须提供包含样本中所有分类群的系统发育树。每轴突功能丰度表对应于MAG或环境中存在的其他已知分类群的基因组注释。这种格式与功能丰度表形成对比，功能丰度表在很大程度上与特定分类群没有联系，例如通过对广泛分布的基因家族数据库的读取映射生成的功能丰度表。POMS的主要输出是一个表格，总结了每个带注释的函数的一致浓缩测试结果。

工作流首先在树的左侧和右侧使用足够的基本提示（默认情况下为10）标识树中的所有节点。所有不符合这些标准的节点都将从分析中排除。然后计算剩余节点的分类群平衡（每个样本），定义为左侧分类群相对丰度与右侧分类群相对丰富度的等距对数比。通常，在分类丰度中添加一个伪计数，以解释缺少分类单元的情况。然后，可以使用标准统计测试来确定这些系统发育平衡在样本组之间是否存在差异。简单地说，样本组之间平衡显著不同的节点表明，与另一个样本组相比，该节点左侧或右侧的分类群相对于节点另一侧的分类群具有显著更高的丰度。根据Wilcoxon秩和检验（未修正）确定样本组之间存在显著差异的余额P（P） < 默认为0.05）。或者，用户可以根据用户定义的外部测试指定哪些节点存在显著差异，从而实现更多自定义分析。

我们将平衡显著不同的每个节点称为BSN。在识别每个BSN后，POMS进一步确定节点的哪一侧在每个样本组中得到了丰富。换句话说，POMS确定在每个样本组中哪些分类群处于相对较高的水平（相对于节点另一侧的分类群）。值得注意的是，这些分类群在一组或另一组中不一定具有较高的相对丰度，而是节点两侧分类群的相对比率在样本组之间不同。

接下来，POMS测试跨节点的每个带注释函数的丰富性。更具体地说，Fisher精确测试是根据节点两侧（未修正）进行或未进行功能编码的尖端计数计算的P（P） < 默认为0.05）。重要的是，这些富集测试是在平衡树步骤中测试的所有节点上计算的：而不仅仅是在BSN上。基于此方法，我们将每个重要节点称为FSN。同样，简单地说，FSN表示测试的功能在节点两侧的分类单元中比另一侧丰富。

最后，使用多项式精确检验来测试FSN是否与BSN一致，并且一致地朝向同一样本组的方向，其频率高于偶然预期。具体来说，对于每个测试功能，该功能的一组已识别FSN（即，与另一个相比，该功能在节点的一侧丰富的节点）可以划分为三类。第一类对应于与BSN不相交的FSN（即，节点一侧的函数显著丰富，但样本平衡没有显著差异）。其他两类对应于与BSN相交的FSN：一类功能在第一个样本组中相对丰富的分类群中富集，另一类功能则在第二个样本组相对丰富的类群中富集。在我们的示例中图1这两类分别由相同和不同颜色的相交圆和正方形表示。对于每个函数，根据每个类中FSN的数量与给定FSN和BSN之间随机交集的预期比例进行多项式精确测试。我们将基于该测试的重要功能称为CEF。

虽然该测试主要用于识别在单个样本组方向富集的CEF，但该框架也可以识别向两个样本组显示混合富集信号的CEF。换言之，函数可能很重要，因为与随机期望相比，与BSN不相交的FSN会减少，而不是不断向特定样本组富集。这种情况在生物学上可能仍然很有趣，但这突出表明，需要对每类FSN的数量进行事后考虑，才能正确解释CEF。

3.2基于仿真的验证

为了验证POMS，我们首先基于包含MAG的样本生成模拟数据集。来自鸟枪宏基因组学（MGS）对照样品的MAG，以及相应的平均深度和系统发育树，是从大型人类肠道MGS meta分析中获得的。这些MAG之前用KEGG正交测井（KO）进行了注释。我们将704个对照样本分为两个大小相等的组1000次，以创建随机测试数据集。

然后，我们进行了几组仿真来评估POMS的性能。对于所有这些模拟，我们将POMS的性能与我们实现的系统发育回归进行了比较。该回归基于Wilcoxon试验确定的组间显著分类，或基于作为系统发育。尽管通过这种比较可以深入了解POMS如何与类似地集成分类和功能信息的现有方法进行比较，但更实际的问题是，与标准差异丰度工具相比，这些方法的性能如何，因为后者更常用。因此，我们还运行了一些差异丰度工具，包括ALDEx2、DESeq2、limma-voom和Wilcoxon测试（基于原始相对丰度或根据USCG丰度修正）。我们将这些标准差异丰度方法应用于缺乏分类学联系的KO丰度表。这些表中每个KO的每样本丰度是通过将编码KO的每个分类单元的分类单元丰度和KO拷贝数的乘积相加来计算的，也就是说，我们计算了分类单元丰富度和KO-拷贝数表的交叉乘积。

为了描述每个工具的基线行为，我们首先调查了显著KO的比例[基于Benjamini–Hochberg校正P（P）-在1000个随机数据集中，每个工具识别的值（BH）<0.05]。一些工具仅在少数数据集中确定了重要的KO；POMS在12个数据集中确定了显著的KO，而ALDEx2和相对丰度以及USCG校正的Wilcoxon检验方法分别确定了7个、17个和19个具有显著KO的数据集。相比之下，DESeq2、limma-voom和基于显著性和特异性的系统发育回归方法分别在543、1000、1000和1000个数据集中发现了显著的KO。尽管三种称为KO的工具在所有数据集中都是重要的，但每个数据集中重要KO的比例很低。limma voom的显著KOs的平均比例为0.163[标准差（SD）=0.038]，基于显著性的系统发育回归为0.060（SD=0.029），基于特异性的系统发育回归为0.056（SD=0.027）。

然后，我们使用三种不同的方法介绍了每个数据集中两个组之间的分类差异。在第一种方法中，对于1000个随机数据集中的每一个，我们随机选择一个由至少五个MAG编码的KO。然后我们模拟选择作用于编码该基因的基因组，方法是在这些基因组中添加一个伪计数，然后仅在一个样本组中将其丰度乘以1.5(补充图S1). 每个随机选择的基因被称为每个数据集的焦点基因，代表一个基因，该基因在两个样本组中的一个样本组对编码它的分类群具有选择性优势。这组模拟数据集被称为焦点基因图谱。在第二种方法中，我们进行了类似的模拟，但选择随机基因组（而不是编码特定焦点基因的基因组）来增加一组中的丰度。为了保持一致性，在每个图谱中，在这种方法中受到干扰的随机基因组的数量与相应的焦点基因图谱中编码焦点基因的基因组的数量相同。这组模拟数据集被称为随机分类群剖面。最后，我们还模拟了同一分支中的所有分类群，而不是随机选择的分类群受到干扰的情况。这涉及到扰乱树中每个节点下的分类单元的丰富性（对于所有至少有五个尖端且不包括根节点的节点）。这导致693个数据集，对应于树中的每个这样的节点。这些数据集被称为基于层的剖面。我们包括这些模拟，以探索POMS在多大程度上对单个分支的生长具有鲁棒性，而不是对编码类似功能的多个独立分支的生长。在生成这些模拟剖面后，我们应用POMS、系统发育回归工作流程和差异丰度方法来识别重要的KO。

我们根据焦点基因相对于给定重复中每种方法确定的所有其他重要KO的显著性排序来评估结果，以便将最高排序（即一个排序）分配给具有最小排序的基因P（P）-值。在我们基于焦点基因的模拟中，MAG编码跨剖面的焦点基因是唯一的直接选择目标，因此，焦点基因的排名应该很高（即接近1，表示较小P（P）-与其他基因相比的价值）。

比较基于模拟复制的焦点基因排名分布的方法，可以发现存在巨大差异(图2a). 最明显的区别是，所有三种系统发育方法确定的焦点基因的排名都显著高于差异丰度工具确定的重点基因。系统发育回归方法中，焦点基因的排名最高，中位数为1，其次是POMS，中位数是2，DESeq2的中位数是65.75。我们在模拟中观察到了相同的总体结果，模拟的极端选择压力较小，MAG数量也有所改变(补充图S2和S3). 此外，我们使用基于参考基因组的模拟稀疏丰度剖面观察到了相同的总体趋势，这表明该结果不是由与关注MAG相关的偏差驱动的（方法和结果；补充图S4和S5).

接下来，我们研究了在模拟复制图谱中，焦点基因排名变化的潜在因素。我们怀疑编码每个焦点基因的MAG数量会显著影响检测能力。事实上，我们发现，根据系统发育方法，除了那些由极少数MAG编码的基因外，焦点基因在总体上排名很靠前(图2b;补充图S6). 例如，POMS确定的在排名前10位的KO中具有相对排名的焦点基因的编码中值为228个MAG，而那些未列入排名前10名的KO的编码中值仅为16.5个MAG。在差异丰度测试中，这一趋势被逆转。例如，排名前10个KO的DESeq2输出中的焦点基因的编码中值为38个MAG，而不包括在前10个ko中的那些编码中值为89个MAG。

这些方法也因其称焦点基因为显著基因的频率而异。特别是，在许多情况下，POMS没有将焦点基因称为显著：我们发现，焦点基因在78.8%的POMS输出中显著，而所有其他工具将其称为显著，平均为95.1%的输出。POMS输出中焦点基因不显著的大多数情况对应于基因由极少数或几乎所有MAG编码的情况(补充图S7). 这一结果突显了POMS与其他方法之间的根本区别：POMS更有可能识别出在独立分类谱系中不断丰富的重要基因，这对于广泛编码（但不普遍存在）的基因来说是最清楚的。相反，差异丰度和某种程度上的系统发育回归方法称基因显著不同，即使这些基因只与少数类群相连。

在证明POMS成功检测到了焦点基因并对其进行了高度排序后，我们接下来将POMS与标准差异丰度测试确定的显著不同KO的总数进行比较。该分析表明，根据POMS，随机类群模拟剖面中显著不同的KO的比例明显最低（平均值=0.001；SD=0.0.017）与所有其他方法相比[例如，次低的是系统发育回归（特异性）：平均值=0.051；标准偏差=0.028；图2c]. 然而，重要的是，基于系统发育方法，随机分类群和焦点基因模拟图谱产生了显著不同的结果，特别是对于POMS，而使用差异丰度方法，这两个图谱之间的差异不太明显与随机分类群数据集相比，POMS在焦点基因数据集中检测到的显著KO的平均比例增加了026%，而基于DESeq2的相应增加（在差异丰度工具中，增加最多）为46.8%。相比之下，POMS在焦点基因数据集中确定的显著KO的平均比例仅比基于分支的分析中确定的比例高166%。这些结果表明POMS和系统发育回归(图2c)，在预期没有一致功能差异的情况下，更好地控制假阳性率，但在有真实信号的情况下仍然可以识别大量重要点击。有趣的是，基于分支的POMS结果中的异常值(图2c)表明当只有一个分支大量转移时，这种方法可能会错误地调用重要基因，至少在极少数情况下是这样。这显然是一个意想不到的结果，值得赞赏（参见讨论）。然而，更令人担忧的是，标准差异丰度测试经常将50%以上的功能确定为显著，即使在样本组差异基于随机分类扰动的情况下。无论参数设置如何，始终观察到这些总体结果(补充图S8).

最后，我们还使用了这些模拟数据集的修改版本（请参见补充方法)评估我们实现的系统发育方法的资源使用情况。我们发现POMS运行时间（在一个核心上）随基因组数量线性增加(补充图S9). 例如，在一个核心上，500个MAG的工作流需要69秒，1000个MAG需要130秒。系统发育回归工作流程显示运行时间在相同数量级上，但在输入MAG数量较大时运行速度略快于POMS。回归工作流程还消耗了POMS两倍多的内存（POMS最大为0.42GB，系统发育回归方法为0.96GB）。

3.3应用于真实数据集

接下来，我们将上述所有方法应用于几个真实的宏基因组数据集。与我们的模拟数据集不同，对于这些数据集，我们对哪些函数将被识别为CEF没有明确的预期，因此我们的评估在解释中容易受到主观偏见的影响。然而，正如我们在下面讨论的那样，鉴于文献中存在类似的发现，POMS确定的许多重要CEF是合理的。

我们首先关注作为Tara Oceans项目一部分组装的MAG。由于强大的选择压力可能会在海洋环境梯度上作用于微生物(伊巴巴兹等。, 2019;萨拉查等。, 2019)这些宏基因组样本代表了POMS的一个有吸引力的测试案例。我们的分析特别关注93个样本，其中642个MAG带有KO注释。与水样相关的环境因素也可用，包括温度、盐度、氧气水平、硝酸盐浓度、磷酸盐浓度和二氧化氮浓度。我们将POMS应用于该数据集，以独立检测与每个因素相关的微生物功能。除KO外，我们还考虑了KEGG途径和模块进行此分析。根据KO注释独立计算每个MAG的这些高级函数（见第2节）。由于环境因素数据是连续的，而不是划分成离散的组，因此我们根据每个节点的余额与因素水平之间的Spearman相关性来识别BSN。根据相关系数的符号推断每个BSN的方向。我们还根据伯克希尔哈撒韦截断值分别为0.05、0.15和0.25，将首席执行官分为严格、中等或宽松三类。

我们对所有环境因素进行了分析，共有一个CEF低于中间截止值，19个CEF小于宽松截止值(表1). 这些CEF与磷酸盐或平均盐度水平有关（分别有13个和11个BSN），事实上，一些确定的功能可能反映了这两个因素的选择性作用。例如，一条与生物膜形成有关的途径ko05111被发现与较高的磷酸盐水平有关，这与之前关于磷水平和生物膜形成之间非平凡关系的研究结果相一致(丹霍恩等。, 2004;方等。, 2009). 与较高磷酸盐水平相关的几个重要的CEF也是膜和转运相关基因，这可能与对环境应激的反应有关。还可以与与平均盐度水平相关的CEF建立合理的联系，例如与柠檬烯和蒎烯降解的正相关（ko00903）。柠檬烯和蒎烯是单萜烯，可以通过代谢改变细胞膜流动性。这可能反映出一种潜在的适应能力，这种适应能力有助于生物体在高盐度条件下生存(De Carvalho和Fernandes，2010年).

接下来，我们研究了如果使用系统发育回归或Spearman相关性来测试关联性，结果会有什么不同（见第2节）。我们发现，使用这些方法会产生与每个因素相关的更多功能。例如，根据系统发育回归，与平均盐度和磷酸盐水平相关的显著KO分别为1008和1348个（BH<0.25）。使用Spearman相关性，平均盐度和磷酸盐水平分别为4169和4974个显著KO。只有10/20和5/20个由POMS鉴定的CEF在Spearman相关性和系统发育回归结果中也具有显著性(表1)分别是。我们还根据Faith的系统发育多样性评估了编码重要功能的分类群的系统发育分布。编码POMS CEF的分类群值稍高（平均值=50.0；SD=27.2），与系统发育回归结果相比（平均值=46.9；SD=44.1），但差异不显著（Wilcoxon检验，W=18123，P（P） = 0.069;补充图S10a). 在所有情况下，POMS检测到的CEF并不是这些其他工具中排名最高的点击(补充图S11a).

然后我们研究了POMS在病例对照宏基因组数据集上的表现。我们重点研究了一个MAG数据集，该数据集是根据人类相关微生物群编译而成，并作为大规模荟萃分析的一部分发布的。我们使用该数据集的子集对应于三个疾病数据集：两个肥胖数据集和一个结直肠癌数据集。

我们分析的第一个肥胖数据集包括477名肥胖者和257名对照者的粪便样本，共含有1401个MAG。我们将POMS应用于该数据集，并确定了肥胖个体和对照个体之间不同的34个BSN。分析的第二个肥胖数据集包括251名肥胖者和159名对照者，他们的粪便微生物群中总共含有1161个MAG，其中31个是BSN。POMS分别鉴定了79、2和21个KO、途径和模块，这些都在至少一个数据集中得到了持续丰富(补充表S1). 在这两个数据集中，两个KO、无通路和两个模块被称为CEF。

重要的是，考虑到我们对人类微生物群与肥胖之间的联系的了解，这些重大发现中的许多都是合理的。例如，POMS检测到的最强CEF之一是与病例样本相关的细胞色素bd-泛醌醇氧化酶（M00153）模块。这一功能可能反映了肥胖个体通过氧化磷酸化或潜在的氧化应激适应而产生能量的广泛变化(朱夫雷等。, 2014). 由于饮食或代谢差异，其他几个CEF也符合类似的广泛选择压力。例如，参与维生素B1和B2生物合成的KO（K00941和K00794）以及参与维生素K2生物合成的模块（M00116）是病例样本中的CEF。这些维生素与肥胖之间的联系是有争议的，但在某些情况下，肥胖者的维生素水平较低(冈安蒂等。, 2014;拉维拉等。, 2021). 同样，β-内酰胺耐药性（ko01501）也与病例样本显著相关，这一点特别有趣，因为长期以来已知接触抗生素与肥胖有关(威尔金斯和赖默，2021). 因此，与这些结果一致的一个假设是，这些功能代表了肥胖个体对应激和肠道环境中资源景观的广泛适应。

最后，我们将POMS应用于75名结直肠癌患者和53名对照的粪便样本的微生物谱，其中包含1187个MAG。在这种情况下，只有14个BSN，没有重要的CEF。尽管如此，结果中唯一的异常值（唯一经过修正的特征P（P） < 0.6）是乙醛酸和二羧酸代谢（在对照组患者中富集；纠正P（P） = 0.294). 这一异常值值得注意，因为此前根据代谢物特征，它被认为是结直肠癌样本（相对于对照组）中功能耗竭最严重的异常值(阿里马等。, 2020).

与POMS结果相反，系统发育回归（基于特异性）在结直肠癌数据集（2996）中的命中率远高于肥胖数据集（1202）。此外，没有路径，只有14/79个KO和4/21个模块(补充表S1)在同一数据集中，这两种方法都称为显著。在这种情况下，编码POMS CEF的分类群显示出Faith系统发育多样性的显著较高值（平均值=114.9；标准偏差=69.9；Wilcoxon检验，W=80348；P（P） = 4.0 × 10⁻⁷)与编码系统发育回归点击数的序列相比（平均值=80.2；标准偏差=74.0；补充图S10b). 此外，根据系统发育回归或差异丰度方法，POMS CEF并不是排名最高的重要呼叫(补充图S11b–d). 这些结果突出表明，与其他系统发育方法和标准差异丰度测试相比，POMS识别出不同的显著点击集。

4讨论

在此，我们提出并验证了POMS框架：一种识别微生物组数据中CEF的新方法。通过模拟，我们已经证明POMS可以准确识别广泛编码的微生物功能，这些功能赋予了强大的选择性优势。虽然我们提出了几个不同的分析来验证和证明我们的方法，但也许最令人信服的证据是，焦点基因（即被模拟为具有选择性优势的基因）是POMS确定的最重要的功能之一，然而，标准差异丰度方法通常不是这样。另一方面，焦点基因也经常通过系统发育回归进行识别，但更难解释这种方法的输出。此外，尽管POMS在组装到MAG中的实际鸟枪宏基因组数据集上产生了合理的结果，但基于系统发育回归，这些CEF通常不是最热门的。这表明POMS和系统发育回归是互补的工具，有可能深入了解不同类型的功能关联。

例如，POMS只能识别由分类群广泛编码的功能，并且这些功能在不同谱系中存在差异；限制于单个谱系的基因，即使该谱系有许多成员，也不应被识别为CEF（尽管这种情况可能发生：请参阅下面基于谱系的模拟讨论）。这是因为POMS要求在独立的BSN中重复（并在一致的方向上）丰富函数，以确定函数的重要性。这些函数必然会在树中较深的节点上显示信号，而不是在小分支或尖端上。这种差异反映在这样一个事实上，编码重要POMS点击数的分类群与编码重要系统发育退化点击数的类群相比，往往具有明显更高的Faith系统发育多样性。因此，POMS可能会丢失许多通过系统发育回归确定的功能，并可能在树的零星尖端显示信号。然而，由POMS识别的CEF有一个明确的生物学解释：每个CEF代表一个假设，即该功能在某些环境中为编码它的分类群提供适合度优势。这不一定是通过系统发育回归确定的重大点击的情况：少数分类群可能只是因为不同的原因而开花或枯竭。

还应该认识到，POMS确定的重要CEF并不能保证只在一个方向上得到丰富。POMS多项检验评估FSN在三个类别中的分布是否偏离随机预期。两种类型的BSN都可以简单地丰富重要功能，即两个样本组中BSN的混合相对较高。这些病例在生物学上可能仍然很有趣，但与主要针对单个样本组富集的CEF相比，其解释会有所不同。因此，在解释POMS确定的任何CEF时，应考虑每类FSN的计数。

POMS方法还有其他限制。例如，确定首席执行官需要有足够的BSN来确定重大浓缩。即使一个适应性功能在系统发育上广泛分布，它也无法被POMS识别，除非有相应的BSN可以由该功能驱动。此外，POMS假定树中的所有节点都具有独立的平衡分布。这在一定程度上是无效的，因为特定的分类群更有可能因个体而异，而且即使在很长的进化距离上也可能发生分类单元共现(妈妈等。, 2020). 对于所有平衡树方法，也需要注意的是，节点平衡中可能存在相关性，这仅仅是因为节点背后的提示相同。我们通过基于枝的扰动模拟评估了这个问题的程度，这表明由于这个问题，POMS有时会产生错误的结果。然而，重要KO的比例仍远低于差异丰度工具。未来的发展重点是基于藻菌化建立类似的框架(洗脸盆等。, 2019)而不是平衡树可以帮助解决这个问题。

尽管有这些警告，但我们已经表明，与更常见的方法相比，将功能信息集成到平衡树框架中可以更好地识别可以提供选择性优势的功能。POMS只是这个通用分析方案如何实现的一个例子，未来的工作可能会结合更复杂的方法来分析跨系统发育树的平衡(西尔弗曼等。, 2017;洗脸盆等。, 2019). 尽管如此，与将差异丰度测试应用于整个社区的功能丰度的常见做法相比，目前的POMS框架代表了一种更稳健的方法。更一般地说，考虑到流行的微生物群差异丰度方法产生的广泛结果(正在接近等。, 2022)而且它们在我们的模拟实验中表现不佳，显然在这种情况下应该首选系统发育方法。

致谢

作者感谢Cecilia Noecker提供的有益反馈。

基金

G.M.D.由国家科学与工程研究委员会（NSERC）加拿大研究生奖学金（博士）资助。G.M.D.还得到了NSERC Michael Smith外国研究补充和Mitacs Globalink研究奖的支持，以实现与E.B.M.G.I.L.的合作，并获得了NSERC-Discovery Grant的支持。E.B.是特拉维夫大学Edmond J.Safra生物信息学中心的教员研究员，部分资金由ISF 2435/19号拨款提供。

利益冲突：未声明。

数据可用性

Tara Oceans MAG数据集来自FigShare（10.6084/m9.figrahare.4902923.v1）。人类病例对照数据集可在NCBI短读本档案中查阅，编号为ERP002061（肥胖#1）、ERP003612（肥胖#2）和ERP012177（结直肠癌）。我们的分析基于MAG和相关数据文件，这些数据文件是通过对这些研究和其他研究的荟萃分析获得的(阿尔梅达等。, 2019). 我们描述的用于运行模拟和实际数据示例的准备好的文件，以及应用于实际数据集的系统发育方法的输出表，可在FigShare上获得，网址为10.6084/m9。

工具书类

作者注释