广告

  • 加载指标

开放式访问

同行评审

研究文章

检测临床宏基因组样本差异丰富特征的统计方法

摘要

目前正在进行大量研究，以确定居住在我们世界上的微生物群落的特征。这些研究旨在极大地扩展我们对微生物生物圈的理解，更重要的是，希望揭示我们与共生菌群之间复杂共生关系的秘密。这些发现的一个重要前提是计算工具能够快速准确地比较复杂细菌群落生成的大型数据集，以识别区分它们的特征。

我们提出了一种统计方法，用于根据计数数据（例如通过测序获得的数据）比较两个治疗人群的临床宏基因组样本，以检测差异丰富的特征。我们的方法，Metastats公司，使用错误发现率来提高高复杂度环境中的特异性，并使用Fisher精确检验分别处理稀疏样本特征。在各种模拟下，我们表明Metastats公司与以前使用的方法相比，性能良好，对于稀疏计数的特征，显著优于其他方法。我们在多个数据集上证明了我们的方法的实用性，包括对肥胖和瘦削人类肠道微生物组的16S rRNA调查，婴儿和成熟肠道微生物组COG功能谱，以及从85个宏基因组的随机序列推断出的细菌和病毒代谢子系统数据。将我们的方法应用于肥胖数据集揭示了原始研究中未报告的肥胖和苗条受试者之间的差异。对于COG和子系统数据集，我们首次对这些人群之间的差异进行了严格的统计评估。本文中描述的方法是第一个解决由多个受试者样本组成的临床宏基因组数据集的方法。我们的方法在不同复杂度和采样水平的数据集上是稳健的。虽然设计用于宏基因组应用，但我们的软件也可以应用于数字基因表达研究（例如SAGE）。我们的方法和免费源代码的web服务器实现可以在网址：http://metats.cbcb.umd.edu/.

作者摘要

宏基因组学的新兴领域旨在仅通过DNA分析来了解微生物群落的结构和功能。目前对社区进行的宏基因组学研究类似于对来自两个普通人群（如患病和健康人群）的多个受试者进行大规模临床试验。为了改进这类实验数据的分析，我们开发了一种统计方法，用于检测微生物群落之间差异丰富的特征，即在一个种群与另一个种群中丰富或贫化的特征。我们表明我们的方法适用于从分类信息到功能注释的各种宏基因组数据。我们还评估了瘦人和肥胖人之间肠道微生物群的分类差异，以及成熟和婴儿肠道微生物群以及微生物和病毒宏基因组功能能力之间的差异。我们的方法是第一个在多个受试者的比较宏基因组学研究中从统计学上解决差异丰度的方法，我们希望能让研究人员更全面地了解两种环境到底是如何不同的。

引用：White JR、Nagarajan N、Pop M（2009）《检测临床宏基因组样本中不同丰富特征的统计方法》。公共科学图书馆计算生物学5（4）：e1000352。https://doi.org/10.1371/journal.pcbi.1000352

编辑：Christos A.Ouzounis，英国伦敦国王学院

收到：2008年10月9日；认可的：2009年3月9日；出版：2009年4月10日

版权所有：©2009 White等人。这是一篇根据知识共享署名许可证条款分发的开放存取文章，允许在任何媒体上不受限制地使用、分发和复制，前提是原创作者和来源得到了认可。

基金：部分资金由比尔和梅琳达·盖茨基金会和亨利·杰克逊基金会提供。资助者在研究设计、数据收集和分析、决定出版或编写手稿方面没有任何作用。

竞争利益：提交人声明，不存在相互竞争的利益。

介绍

高通量、廉价测序技术的日益普及，催生了一个新的科学领域，即宏基因组学，包括微生物群落的大规模分析。细菌种群的广泛测序让我们第一眼看到许多无法通过传统方法进行分析的微生物（用当前方法只能分离和独立培养所有细菌的～1%[1]). 环境样品的研究最初侧重于单个基因的靶向测序，特别是核糖体RNA的16S亚基[2]–[5]尽管最近的研究利用了高通量鸟枪测序方法，不仅评估了微生物群落的分类组成，还评估了其功能能力[6]–[8].

近年来开发了几种软件工具，用于根据序列数据比较不同的环境。DOTUR公司[9]、利舒夫[10]，б-libshuff[11]、儿子[12]、梅根[13]、UniFrac[14]和TreeClimer[15]所有这些都集中在这样一个分析的不同方面。DOTUR将序列聚类为操作分类单元（OTU），并提供微生物种群多样性的估计值，从而为比较不同群落提供了粗略的衡量标准。SON扩展了DOTUR，并提供了一个统计数据，用于估计两个环境之间的相似性，特别是两个社区之间共享的OTU的比例。Libshuff和Ş-Libshuff提供了一个假设检验（Cramer von Mises统计）来决定两个群落是否不同，TreeClimber和UniFrac将这个问题置于系统发育背景下。注意，这些方法旨在评估是否，而不是怎样这两个社区不同。当我们开始分析微生物组对人类健康的贡献时，后一个问题尤其重要。临床试验中的宏基因组分析需要在个体分类水平上的信息来指导未来的实验和治疗。例如，我们想确定是否存在会导致人类疾病的细菌，并开发抗生素或益生菌治疗方法。这个问题首先由罗德里格斯-布里托提出等。 [16]，他们使用bootstrapping来估计与生物子系统丰度差异相关的p值。最近，Huson的MEGAN软件等。 [13]提供了一个图形界面，允许用户比较不同环境的分类组成。注意，MEGAN是上述程序中唯一一个可以应用于16S rRNA调查以外的数据的程序。

这些工具有一个共同的局限性——它们都是为精确比较两个样本而设计的——因此在临床环境中的适用性有限，因为临床环境的目标是比较两个（或多个）治疗群体，每个群体包含多个样本。在本文中，我们描述了一种严格的统计方法，用于检测临床宏基因组数据集之间差异丰富的特征（分类群、通路、子系统等）。该方法适用于高通量宏基因组数据和16S rRNA调查。我们的方法扩展了最初为微阵列分析开发的统计方法。具体来说，我们将这些方法应用于离散计数数据并校正稀疏计数。我们的研究受到了宏基因组项目对临床应用日益关注的推动（例如人类微生物组项目[17]).

注意，在数字基因表达研究（例如SAGE）的背景下，已经解决了类似的问题[18]). 卢等。 [19]采用过分散对数线性模型和Robinson和Smyth[20]在分析多个SAGE库时使用负二项分布。这两种方法都可以应用于宏基因组数据集。我们通过综合模拟将我们的工具与这些先前的方法进行了比较，并通过分析公开的数据集证明了我们方法的性能，这些数据集包括16S人类肠道微生物群调查和基于随机序列的婴儿和成熟肠道微生物群、微生物和病毒宏基因组功能调查。本文中描述的方法已作为web服务器实现，也可以从中作为免费源代码（在R中）使用http://metastats.cbcb.umd.edu.

材料和方法

我们的方法基于以下假设：（i）我们获得了与两个治疗人群（例如，患病和健康的人类肠道社区，或接触不同治疗的个人）相对应的数据，每个群体由多个个体（或样本）组成；（ii）对于每个样品，我们都获得了代表特定物质相对丰度的计数数据特征在每个样本中，例如分配给特定分类单元的16S rRNA克隆的数量，或映射到特定生物途径或子系统的鸟枪读数的数量（参见下文如何使用当前可用的软件包生成此类信息）。我们的目标是确定此类数据集中区分两个种群的个别特征，即两个种群中丰度不同的特征。此外，我们还开发了一种统计方法来衡量观察到的差异的可信度。

我们方法的输入可以表示为特征丰度矩阵其行对应于特定特征，其列对应于单个宏基因组样本。中的单元格我第个行和j个第个列是特征的观察总数我在样品中j个(图1). 每个不同的观察值在矩阵中只表示一次，即不允许重叠特征（行对应于序列集的一个分区）。

缩略图
下载：
图1。特征丰度矩阵的格式。

每一行代表一个特定的分类单元，而每一列代表一个主题或复制。The frequency of the我第个中的功能j个第个主题(c（i，j）)记录在矩阵的相应单元格中。如果有克在第一个群体中，他们由第一个群体代表克矩阵的列，而其余列表示第二个总体的主题。

https://doi.org/10.1371/journal.pcbi.1000352.g001

数据规范化

为了解释多个个体的不同采样水平，我们将原始丰度度量转换为代表每个特征对每个个体的相对贡献的分数。这导致上述矩阵的标准化版本，其中我第个行和j个第个列（我们将用它来表示（f）ij公司)是分类单元的比例我在个体中观察到j个。我们选择这种简单的归一化程序是因为它提供了计数数据的自然表示，作为相对丰度度量，然而，可以使用其他归一化方法来确保观测到的计数在样本之间具有可比性，我们目前正在评估几种此类方法。

差异丰度分析

对于每个功能我，我们通过计算两个样本来比较两个处理群体中的丰度吨统计的。具体来说，我们计算平均比例、和方差每种治疗吨其中n个吨对受试者（矩阵中的列）进行采样：然后我们计算两个样本吨统计的：其功能吨超过特定阈值的统计数据可以推断为在两个处理中差异丰富（双面吨-测试）。

评估重要性

的阈值吨选择统计数据是为了最小化误报的数量（错误地确定差异丰富的特征）。具体来说，我们试图控制p值——观察给定值的可能性吨随机统计。传统分析使用吨具有适当数量的自由度的分布。然而，这个过程的一个隐含假设是，基本分布是正态的。我们不做这个假设，而是估计吨我使用Storey和Tibshirani中描述的置换方法进行非参数化[21]。此过程也称为非参数吨-试验表明，当基本分布为非正态分布时，可以准确估计显著性[22],[23]具体而言，我们随机排列丰度矩阵的列的处理标签，并重新计算吨统计数据。注意，排列保持如下n个1重复治疗1和n个2重复治疗2。对重复此过程B类试验，我们获得B类套吨统计数据：吨10亿, …,吨M（M）0磅,b条 = 1, …,B类，其中M（M）是矩阵中的行数。对于每一行（特征），与观测值关联的p值吨统计量计算为使用吨统计值大于或等于观察值吨我:这种方法不适用于所有列的可能排列数量有限的小样本大小。作为一种启发式方法，如果两种治疗中使用的受试者少于8名，我们将所有受试者集合在一起吨将统计数据合并为一个空分布，并将p值估计为：

请注意，在我们的方法实施过程中，选择8作为截止值只是基于实验的启发式方法。我们的方法专门针对包含多个主题的数据集，适用于Rodriguez-Brito等人提出的小型数据集方法。[16]可能更合适。

除非明确说明，否则下面描述的所有实验都使用了1000个排列。一般来说，排列的数量应选择为实验中使用的显著性阈值的函数。具体来说，使用B类排列只能估计低至1的p值/B类（在我们的案例中为10−3). 在包含许多特征的数据集中，需要更多的排列来解释多个假设检验问题（下文将讨论此情况的进一步更正）。通过增加用于近似零分布的排列数，p值计算的精度明显提高，但代价是增加了计算时间。对于某些分布，可以使用称为重要性抽样的技术有效地估计较小的p值。具体地说，置换测试的目标是估计分布的尾部，导致所需的置换数量减少了95%[24],[25]。我们打算在未来版本的软件中实施这种方法。

多假设检验修正

对于复杂的环境（许多特征/分类群/子系统）吨所描述的统计数据可能会导致大量假阳性。例如，在包含1000个生物体的数据集中，选择0.05的p值阈值将导致50个假阳性。直观校正包括减少与所执行测试数量成比例的p值截止值（Bonferroni校正），从而减少误报数量。然而，当进行大量测试时，这种方法可能过于保守[21].

另一种方法是控制错误发现率（FDR），它被定义为预测集合中的误报比例[26]，而假阳性率定义为整个测试中假阳性的比例。在这种情况下，测试的重要性通过q值来衡量，q值是对每个测试的FDR的单独测量。

我们基于Storey和Tibshirani使用以下算法计算q值[21]该方法假设真零测试的p值均匀分布，这一假设适用于Metastats中使用的方法。给定p值的有序列表，第页(1)≤第页(2)≤…≤第页(米)，（其中米是功能的总数）和一系列值λ = 0，0.01，0.02，…，0.90，我们计算接下来，我们适应用三次样条曲线表示三个自由度，并让最后，我们估计了每个有序p值对应的q值。第一，.然后针对我 = m-1个,m-2个,…,1,因此，用p值进行假设检验对应的q值为注意，此方法产生真实q值的保守估计，即。。我们的软件为用户提供了使用p值或q值阈值的选项，无论数据的复杂性如何。

处理稀疏计数

对于低频特征，例如低丰度分类群，非参数吨–上述测试不准确[27]。我们进行了几次模拟（未显示数据），以确定非参数的局限性吨-测试稀疏采样特征。相应地，我们的软件仅适用于任一总体中某个特征的观察总数大于总体中受试者总数的情况（即给定特征的跨受试者平均观察数大于一）。我们使用Fisher精确检验比较了稀疏样本（罕见）特征的差异丰度。Fisher的精确测试根据超几何分布（无替换抽样）对抽样过程进行建模。将丰度矩阵中稀疏特征的频率合并，以创建2×2列联表(图2)，用作双尾测试的输入。使用中的符号图2，观测2×2列联表的零超几何概率为：

缩略图
下载：
图2。检测稀疏特征的差异丰度。

Fisher精确检验中使用了2×2列联表来区分稀有特征之间的丰度差异。（f）11是特征的观察数我在治疗1的所有个体中。（f）21是不是特征的观测值的数量我在治疗1的所有个体中。（f）12和（f）22对治疗2的定义类似。

https://doi.org/10.1371/journal.pcbi.1000352.g002

通过计算给定表格以及所有比观察到的表格更极端的表格的概率，可以计算出在假设零假设（即无差异丰度）为真的情况下，偶然获得原始表格的准确概率[27].

注意，在微阵列文献中提出了一种处理稀疏特征的替代方法。微阵列（SAM）方法的显著性分析[28]使用修改的吨统计的。由于我们数据的离散性，以及之前在数字基因表达背景下进行的研究表明，Fisher检验对检测差异丰度有效，因此我们选择使用Fisher精确检验[29].

创建特征丰度矩阵

我们方法的输入，即特征丰度矩阵，可以使用可用的软件包从16S rRNA和随机霰弹枪数据轻松构建。特别是对于16S分类分析，RDP Bayesian分类器等工具[30]和Greengenes SimRank[31]输出关于样本中每个分类单元丰度的易于解析的信息。作为一种补充的无监督方法，16S序列可以与DOTUR进行聚类[9]操作分类单元（OTU）。丰度数据可以很容易地从“*.list”文件中提取，该文件详细说明了哪些序列是同一OTU的成员。使用MEGAN可以对霰弹枪数据进行功能或分类分类[13]、CARMA[32]或MG-RAST[33].MEGAN和CARMA都能够输出分配给分类法或功能组的序列列表。MG-RAST为代谢子系统提供了类似的信息，这些信息可以作为tab分隔文件下载。

上述所有数据类型都可以很容易地转换为适合作为我们方法输入的特征丰度矩阵。未来，我们还计划为通用分析工具生成的数据提供转换器。

本文使用的数据

按照Eckburg的描述制备了人类肠道16S rRNA序列等。和莱伊等。（2006）和可在GenBank中获得，登录号：DQ793220-DQ802819、DQ803048、DQ802139-DQ810181、DQ823640-DQ825343、AY974810-AY986384。在我们的实验中，我们使用核糖体数据库项目II（RDP）目前使用的朴素贝叶斯分类器将所有16S序列分配给分类群[30]从Kurokawa的补充材料中获得了13种人类肠道微生物的COG谱等。 [34]。我们从Dinsdale的在线补充材料中获得了85个宏基因组的代谢功能图谱等。(2008) (http://www.these.org/Dinsdale补充材料/).

结果

与其他统计方法的比较

如介绍，用于分析SAGE数据的统计软件包也适用于宏基因组数据集。为了验证我们的方法，我们首先设计了模拟，并将转移的结果与Student的结果进行了比较吨-SAGE数据的检验（合并方差）和两种方法：Lu的对数线性模型（log-t）等。 [19]以及Robinson和Smyth开发的负二项（NB）模型[20].

我们设计了一项宏基因组模拟研究，从两组中抽取10名受试者，每个受试者的取样深度通过从200到1000的均匀分布中随机取样确定（这些深度对于宏基因组研究来说是合理的）。给定人口平均比例第页和分散值φ，我们从β二项分布中抽取序列Β(α,β)，其中α = 第页（1）/φ−1）和β = (1−第页)(1/φ−1). 注意，此抽样程序的数据符合Lu的假设等。以及罗宾逊和斯迈思，因此我们希望他们在这些条件下表现出色。卢等。为SAGE数据设计了一个类似的研究，然而，对于每个模拟，两个种群都使用了固定的离散度，离散度估计值非常小(φ = 0、8e-06、2e-05、4.3e-05）。尽管这些值对于SAGE数据可能是合理的，但我们发现它们并不能准确地模拟宏基因组数据。图3显示以下检查的宏基因组数据集的所有特征在每个群体中的估计离散度。离散度估计范围从1e-07到0.17，这两个种群很少有共同的离散度。因此，我们设计了模拟，以便φ从1e-08和0.05之间的均匀分布中随机选择每个种群，考虑到种群分布之间的潜在显著差异。对于每组参数，我们模拟了1000个特征计数，其中500个是在第页1 = 第页2，其余部分差异丰富，其中a*p公司1 = 第页2，并使用接收器操作特征（ROC）曲线比较了每种方法的性能。图4显示一系列值的ROC结果第页和一。对于每组参数，Metastats使用5000个置换来计算p值。Metastats的性能与其他方法一样好，在某些情况下更可取。我们还发现，在大多数情况下，我们的方法比负二项模型更为敏感，后者在高丰度特征方面表现不佳。

缩略图
下载：
图3。离散度估计(φ)用于本研究中使用的三个宏基因组数据集。

这些图比较了（A）肥胖和瘦肉人类肠道分类数据、（B）婴儿和成熟人类肠道COG分配以及（C）微生物和病毒子系统注释之间的离散值。我们在这一数据中发现了广泛的可能分散性，并且两个群体之间的分散性存在显著差异。

https://doi.org/10.1371/journal.pcbi.1000352.g003

缩略图
下载：
图4。模拟研究中比较统计方法的ROC曲线。

序列选自具有可变分散度和群平均比例的β二项式分布第页1和第页2。对于每一组参数，我们模拟了1000个试验，其中500个是在零假设下生成的(第页1 = 第页2)，其余的则有差异地丰富，其中a*p公司1 = 第页2例如，p=0.2和a=2表示包含20%种群的特征，两个感兴趣种群之间的丰度相差两倍。的参数值第页1和一显示在每个图的上方。

https://doi.org/10.1371/journal.pcbi.1000352.g004

我们的下一个模拟旨在检查每种方法在极端稀疏采样下的准确性。如下面的数据集所示，通常情况下，一个特征在一个群体中可能没有任何观察结果，因此必须使用一种统计方法来解决宏基因组数据的这种频繁特征。在与上述模拟相同的假设下，我们进行了测试一 = 0和0.01，从而显著减少了对其中一个种群特征的观察。ROC曲线如所示图5揭示了Metastats在极端稀疏性方面优于其他统计方法。保持假阳性率（x轴）不变，Metastats显示出比所有其他方法更高的灵敏度。值得注意的是，Log-t的性能较差，因为它是为同样可能稀疏的SAGE数据设计的。进一步的调查表明，如果在一个群体中没有观察结果，Log-t方法会导致高度膨胀的离散值，从而降低测试的估计重要性。

缩略图
下载：
图5。极端稀疏采样模拟研究中比较统计方法的ROC曲线。

序列从具有可变离散度和组平均比例的β二项分布中选择第页1和第页2。对于每一组参数，我们模拟了1000个试验，其中500个是在零假设下生成的(第页1 = 第页2)，其余的则有差异地丰富，其中a*p公司1 = 第页2例如，p=0.2和a=2表示包含20%种群的特征，两个感兴趣种群之间的丰度相差两倍。的参数值第页1和一显示在每个图的上方。

https://doi.org/10.1371/journal.pcbi.1000352.g005

最后，我们选择了Dinsdale的一个子集等。 [6]宏基因组子系统数据（如下所述），并将每个受试者随机分配到两个群体中的一个（每个群体20名受试者）。所有受试者实际上来自相同的人群（微生物宏基因组），因此，对于每个测试的特征，零假设都是正确的（没有特征差异丰富）。我们对这些数据运行了每种方法，记录了每个特征的计算p值。重复这个过程200次，我们模拟了5200个空特性的测试。表1显示给定不同p值阈值的每种方法产生的假阳性数量。结果表明，与其他方法相比，负二项模型导致了异常高的误报率。学生的吨-test和Metastats在估计这些空特性的重要性方面表现同样好，而Log-t表现稍好。

缩略图
下载：
表1。比较不同方法发现的假阳性。使用真实的宏基因组数据，我们通过将受试者从单个群体随机分为两个亚群体来模拟无差异丰度的特征。

https://doi.org/10.1371/journal.pcbi.1000352.t001

这些研究表明，Metastats在深度采样特征方面的表现与所有其他适用方法一样一致，并且在稀疏数据方面优于这些方法。下面我们进一步评估Metastats在几个真实的宏基因组数据集上的性能。

与人类肥胖相关的类群

在最近的一项研究中，Ley等。 [35]确定了与人类肥胖相关的肠道微生物，并得出结论，肥胖具有微生物成分，特别是厚壁菌和拟杆菌是瘦人和肥胖人之间差异丰富的细菌分支。肥胖受试者的厚壁菌相对丰度显著高于瘦肉受试者，而细菌类相对丰度则较低。此外，肥胖受试者被置于限制热量的饮食中一年，之后，受试者的肠道微生物群与瘦人更为相似。

我们获得了在Ley中测序的20609个16S rRNA基因等。并将其划分为不同分辨率的分类群（请注意，16S序列中的2261个来自之前的研究[36]). 我们最初试图用我们的方法重新建立本文的主要结果。表2说明我们的方法与原始研究的结果一致：肥胖症受试者中厚壁菌显著增多（P=0.003），瘦肉人群中拟杆菌显著增多（P<0.001）。此外，我们的方法还检测到放线菌差异丰富，这一结果在最初的研究中没有报道。肥胖受试者中约5%的样本由放线菌组成，而瘦肉受试者的放线菌频率明显较低（P=0.004）。科林塞拉和埃格特拉是观察到的最常见的放线菌属，这两种放线菌在肥胖受试者中都过多。众所周知，这些生物体会将糖发酵成各种脂肪酸[37]，进一步加强了与肥胖的可能联系。注意，最初的研究使用了学生吨-测试，导致放线菌内观察到的差异的p值为0.037，比我们的计算结果大9倍。这突出了我们方法的敏感性，并解释了为什么最初没有检测到这种差异。

缩略图
下载：
表2。瘦肉型和肥胖型人类肠道微生物区系的不同丰富分类群。

https://doi.org/10.1371/journal.pcbi.1000352.t002

为了探索我们是否可以完善初步研究的广泛结论，我们在更详细的分类水平上重新分析了数据。我们确定了三类差异丰富的生物体：梭状芽孢杆菌（P=0.005）、拟杆菌（P<0.001）和放线菌（P=0.003）。这三个是相应门（厚壁菌属、拟杆菌属、放线菌属）的主要成员，其分布与在较粗水平上观察到的分布相同。Metastats还检测到9个差异丰富的属，占两个种群16S序列的25%以上（P≤0.01）。互营球菌属,瘤胃球菌、和科林氏菌属在肥胖受试者中都得到了强化，而拟杆菌类平均而言，瘦身受试者的营养水平是正常人的八倍。

对于每个主题都有多个观察结果的分类群，我们发现我们的结果（p值估计）与用大多数其他方法获得的结果之间有很好的一致性(表2). 令人惊讶的是，我们发现Robinson和Smyth的负二项模型未能在这些数据集中检测到几个显著差异丰富的特征（例如，厚壁菌的假设检验结果的p值为0.87）。这可能部分是由于难以为我们的数据集估计其模型的参数，并进一步加强了针对宏基因组数据特征设计方法的案例。对于特定分类单元在一个种群中没有观察到的情况（例如。特拉萨基亚菌)，针对SAGE数据提出的方法似乎表现不佳。

成熟和婴儿人类肠道微生物组之间不同丰富的COG

16S rRNA的靶向测序只能提供微生物群落内多样性的概述，但不能提供有关该群落成员功能作用的任何信息。环境的随机鸟枪测序可以让我们一窥环境中生物体基因编码的功能复杂性。定义环境功能容量的一种方法是将猎枪序列映射到具有已知功能的同源序列。这个策略是黑川方明使用的等。 [34]以确定13个人（包括4名未断奶婴儿）肠道微生物群中的直向同源群（COG）集群。我们在所有受试者中检查了本研究确定的COG，并使用Metastats发现婴儿和成熟（>1岁）肠道微生物组之间存在差异丰富的COG。这是这两个群体的首次直接比较，因为最初的研究只是将每个群体与参考数据库进行比较，以找到丰富的基因集。由于该数据集测试的特征数量较多（3868个COG），且可用婴儿受试者数量有限，我们的方法使用池选项计算p值（我们选择了100个排列），然后计算每个特征的q值。使用Q≤0.05的阈值（将错误发现率控制在5%），我们检测到192个COG，在这两个种群之间差异丰富（参见表3获得成熟和婴儿微生物组中最丰富的COG许可证。完整结果作为补充材料显示在表S1).

缩略图
下载：
表3。使用q值阈值0.05，确定婴儿和成人肠道微生物组之间的COG差异丰富。

https://doi.org/10.1371/journal.pcbi.1000352.t003

成熟受试者中最丰富的COG包括信号转导组氨酸激酶（COG0642）、外膜受体蛋白，如铁转运（COG1629）和β-半乳糖苷酶/β-葡萄糖醛酸酶（COG3250）。这些COG在婴儿中也相当丰富，但相对于成熟受试者而言耗尽。婴儿保持与糖转运系统（COG1129）和转录调控（COG1475）相关的丰富COG。在最初的研究中也发现了糖转运功能的过度丰富，这加强了这样一个假设，即未经驯化的婴儿肠道微生物群是专门为消化母乳中的单糖而设计的。同样，婴儿体内铁转运蛋白的耗竭可能与母乳中相对于牛奶的低铁浓度有关[38]。尽管铁的浓度很低，但婴儿从母乳中吸收的铁非常高，并且在婴儿断奶后会变得更差，这表明了摄入这种矿物质的另一种机制。检测到一种铁氧还蛋白样蛋白（COG2440），该蛋白在婴儿体内的含量是成熟受试者的11倍，而铁氧还毒素（COG1145）在成年受试者体内的含量显著增加，这证明了可能存在不同的机制。

微生物和病毒宏基因组中不同丰富的代谢子系统

丁斯代尔的最新研究等。使用SEED平台分析了87个不同的宏基因组鸟枪样本（~1500万序列）(http://www.theseed.org网站)[6]看看生物地球化学条件是否与宏基因组特征相关。我们从本研究中分析的45个微生物和40个病毒宏基因组中获得了功能图谱。在丁斯代尔分析的26个子系统（抽象功能角色）中等。研究发现，微生物和病毒样本之间有13个显著差异（P≤0.05）(表4). RNA和DNA代谢的子系统在病毒宏基因组中显著丰富，而氮代谢、膜转运和碳水化合物都在微生物群落中富集。病毒宏基因组中高水平的RNA和DNA代谢说明它们需要一个自给自足的核苷酸来源。虽然原始研究所描述的差异不包括显著性估计，但我们的结果基本上与作者的定性结论一致。然而，由于自首次发布以来SEED数据库中的注释不断更新，我们发现我们的结果与最初报告的结果之间存在一些差异。特别是，我们发现毒力子系统的总体丰度低于之前的报道，并且在微生物和病毒宏基因组之间的丰度没有发现任何显著差异。

缩略图
下载：
表4。微生物和病毒宏基因组之间差异丰富的代谢子系统（平均百分比±标准差，p值≤0.05）。

https://doi.org/10.1371/journal.pcbi.1000352.t004

讨论

我们提出了一种处理频率数据的统计方法，以检测两个种群之间差异丰富的特征。该方法可用于分析通过分子方法生成的任何计数数据，包括环境样品的随机鸟枪测序、元基因组样品中特定基因的定向测序、数字基因表达调查（例如SAGE[29])，甚至是全基因组鸟枪数据（例如，比较组装基因的测序覆盖深度）。在模拟数据集和真实数据集上的比较表明，当应用于采样良好的数据集时，我们的软件的性能与其他统计方法相当，并且在稀疏数据上优于这些方法。

我们的方法也可以通过替换吨-用单因素方差分析测试。此外，如果每个处理中只有一个样本可用，则可以使用二次方检验代替吨-测试。[27].

未来几年，宏基因组研究将越来越多地应用于临床，需要开发新的算法和软件工具，以利用数百到数千名患者的数据。上述方法是朝着这一方向迈出的第一步，提供了一种稳健而严格的统计方法，用于识别不同丰度与疾病相关的生物体和其他特征。这些方法、相关的源代码和我们工具的web界面都可以在http://metastats.cbcb.umd.edu.

支持信息

表S1。

成熟和婴儿微生物组COG组成的比较。每个COG的相对丰度与显著性值（丰度差异仅由偶然性引起的可能性）一起表示，显著性值是使用以下几种统计程序计算得出的：亚统计（p和q值）、学生t检验、对数线性模型和负二项模型。

https://doi.org/10.1371/journal.pcbi.1000352.s001

（0.95 MB XLS）

致谢

我们非常感谢Jeff Gordon、Peter Turnbaugh和Ruth Ley在数据方面的投入和帮助。我们感谢Frank Siewerdt、Aleksey Zimin、Radu Balan、Kunmi Ayanbule和Kenneth Ryals的有益讨论。最后，我们要感谢匿名评论员的建设性意见。他们的贡献大大改进了我们的手稿。

作者贡献

构思并设计实验：JRW NN MP。执行实验：JRV。分析数据：JRW NN MP。撰写论文：JRW-NN MP。

工具书类

  1. 1Schloss PD，Handelsman J（2005）《研究不可培养微生物的宏基因组学：解决难题》。基因组生物学6:229。
  2. 2Bik EM、Eckburg PB、Gill SR、Nelson KE、Purdom EA等（2006），《人体胃中细菌微生物群的分子分析》。美国国家科学院院刊103:732-737。
  3. 三。Batzoglou S、Jaffe DB、Stanley K、Butler J、Gnere S等（2002）《阿拉伯：一个全基因组霰弹枪组装器》。基因组研究12:177–189。
  4. 4Palmer C、Bik EM、Digiulio DB、Relman DA、Brown PO（2007）《人类婴儿肠道微生物群的发展》。《公共科学图书馆·生物》5:e177。
  5. 5Sogin ML、Morrison HG、Huber JA、Welch DM、Huse SM等（2006）深海和未充分开发的“稀有生物圈”中的微生物多样性。美国国家科学院院刊103:12115–12120。
  6. 6Dinsdale EA、Edwards RA、Hall D、Angly F、Breitbart M等（2008）九个生物群落的功能宏基因组分析。自然452:629–632。
  7. 7Gill SR、Pop M、Deboy RT、Eckburg PB、Turnbaugh PJ等（2006）人类远端肠道微生物组的元基因组分析。科学312:1355-1359。
  8. 8Turnbaugh PJ、Ley RE、Mahowald MA、Magrini V、Mardis ER等（2006），一种能量获取能力增强的肥胖相关肠道微生物群。自然444:1027-1031。
  9. 9Schloss PD，Handelsman J（2005）介绍DOTUR，一个用于定义操作分类单位和估算物种丰富度的计算机程序。应用环境微生物71:1501-1506。
  10. 10Singleton DR，Furlong MA，Rathbun SL，Whitman WB（2001）《环境样品16S rRNA基因序列库的定量比较》。应用环境微生物67:4374–4376。
  11. 11Schloss PD，Larget BR，Handelsman J（2004）微生物生态学与统计学的整合：比较基因库的测试。应用环境微生物70:5485–5492。
  12. 12Schloss PD，Handelsman J（2006）《介绍SONS》，这是一种基于操作分类单元的微生物群落成员和结构比较工具。应用环境微生物72:6773–6779。
  13. 13Huson DH、Auch AF、Qi J、Schuster SC（2007）《宏基因组数据的MEGAN分析》。基因组研究17:377–386。
  14. 14Lozupone C，Knight R（2005）UniFrac：一种用于比较微生物群落的新系统发育方法。应用环境微生物71:8228–8235。
  15. 15Schloss PD，Handelsman J（2006）《引入TreeClimber》，一种比较微生物群落结构的测试。应用环境微生物72:2379–2384。
  16. 16Rodriguez-Brito B，Rohwer F，Edwards RA（2006）《统计学在比较宏基因组学中的应用》。BMC生物信息学7:162。
  17. 17Turnbaugh PJ、Ley RE、Hamady M、Fraser-Ligett CM、Knight R等（2007）人类微生物组项目。自然449:804-810。
  18. 18Velculescu VE，Zhang L，Vogelstein B，Kinzler KW（1995）基因表达的系列分析。科学270:484–487。
  19. 19Lu J，Tomfohr JK，Kepler TB（2005）《识别多个SAGE库中的差异表达：过度分散的对数线性模型方法》。BMC生物信息学6:165。
  20. 20Robinson MD，Smyth GK（2007），评估标签丰度差异的中等统计检验。生物信息学23:2881–2887。
  21. 21Storey JD，Tibshirani R（2003）全基因组研究的统计意义。美国国家科学院院刊100:9440–9445。
  22. 22Troyanskaya OG，Garber ME，Brown PO，Botstein D，Altman RB（2002），识别微阵列数据中差异表达基因的非参数方法。生物信息学18:1454-1461。
  23. 23Efron B，Tibshirani R（1993）《引导程序简介》。纽约：查普曼和霍尔出版社。
    • 24Hesterberg T（1996）有效引导模拟的控制变量和重要性抽样。统计与计算6:147–157。
    • 25Johns MV（1988）《Bootstrap置信区间的重要性抽样》。美国统计协会期刊83:709–714。
    • 26Benjamini Y，Hochberg Y（1995）《控制错误发现率——一种实用且强大的多重测试方法》。英国皇家统计学会期刊B辑方法学57:289–300。
    • 27Zar JH（1999）生物统计分析。新泽西州上马鞍河：普伦蒂斯·霍尔。
      • 28Tusher VG，Tibshirani R，Chu G（2001）电离辐射反应微阵列的显著性分析。美国国家科学院院刊98:5116–5121。
      • 29Ruijter JM、Van Kampen AH、Baas F（2002）《SAGE库的统计评估：实验设计的后果》。《生理基因组学》11:37–44。
      • 30Wang Q，Garrity GM，Tiedje JM，Cole JR（2007）Naive Bayesian分类器，用于将rRNA序列快速分配到新的细菌分类学中。应用环境微生物73:5261–5267。
      • 31DeSantis TZ，Hugenholtz P，Larsen N，Rojas M，Brodie EL，et al.（2006）Greengenes，一个与ARB兼容的嵌合体16S rRNA基因数据库和工作台。应用环境微生物72:5069–5072。
      • 32Krause L，Diaz NN，Goesmann A，Kelley S，Nattkemper TW，et al.（2008）短环境DNA片段的系统发育分类。核酸研究36:2230–2239。
      • 33Meyer F、Paarmann D、D’Souza M、Olson RD、Glass EM等（2008）宏基因组RAST服务器-用于宏基因组自动系统发育和功能分析的公共资源。BMC生物信息学9:386。
      • 34Kurokawa K、Itoh T、Kuwahara T、Oshima K、Toh H等（2007）比较宏基因组学揭示了人类肠道微生物组中普遍富集的基因集。DNA研究14:169–181。
      • 35Ley RE、Turnbaugh PJ、Klein S、Gordon JI（2006）《微生物生态学：与肥胖相关的人体肠道微生物》。自然444:1022–1023。
      • 36Eckburg PB、Bik EM、Bernstein CN、Purdom E、Dethlefsen L等（2005）《人类肠道微生物菌群的多样性》。科学308:1635-1638。
      • 37Tuomanen E（2005）《人类的微生物群落——它们的生态学及其在健康和疾病中的作用》。科学308:635-635。
      • 38Picciano MF，Deering RH（1980）婴儿期喂养方案对铁状态的影响。美国临床营养学杂志33:746–753。