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

数据规范化

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

差异丰度分析

对于每个功能我，我们通过计算两个样本来比较两个处理群体中的丰度t吨统计的。具体来说，我们计算平均比例、和方差每种治疗t吨从中n个_t吨对受试者（矩阵中的列）进行采样：

然后我们计算两个样本t吨统计的：

其功能t吨超过特定阈值的统计数据可以推断为在两个处理中差异丰富（双面t吨-测试）。

评估重要性

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

这种方法不适用于所有列的可能排列数量有限的小样本大小。作为一种启发式方法，如果两种治疗中使用的受试者少于8名，我们将所有受试者集合在一起t吨将统计数据合并为一个空分布，并将p值估计为：

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

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

多假设检验修正

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

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

我们基于Storey和Tibshirani使用以下算法计算q值[21]该方法假设真零测试的p值均匀分布，这一假设适用于Metastats中使用的方法。给定p值的有序列表，第页 ₍₁₎≤第页 ₍₂₎≤…≤第页 _(米)，（其中米是功能的总数）和一系列值λ======================================================================================0，0.01，0.02，…，0.90，我们计算

接下来，我们将用三次样条曲线表示三个自由度，并让最后，我们估计了每个有序p值对应的q值。第一，.然后针对我======================================================================================m-1个,m-2个,…,1,

因此，用p值进行假设检验对应的q值为注意，此方法产生真实q值的保守估计，即。我们的软件为用户提供了使用p值或q值阈值的选项，无论数据的复杂性如何。

处理稀疏计数

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

在单独的窗口中打开

图2

检测稀疏特征的差异丰度。

2×2列联表用于Fisher精确检验稀有特征之间的差异丰度。（f）₁₁是特征的观察数我在来自治疗1的所有个体中。（f）₂₁是非特征的观察数我在治疗1的所有个体中。（f）₁₂和（f）₂₂对治疗2的定义类似。

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

注意，在微阵列文献中提出了一种处理稀疏特征的替代方法。微阵列（SAM）方法的显著性分析[28]使用修改的t吨统计的。由于我们数据的离散性，以及之前在数字基因表达背景下进行的研究表明，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为代谢子系统提供了类似的信息，这些信息可以作为制表符分隔的文件下载。

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

与人类肥胖相关的类群

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

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

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

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

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

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