【原】从基因组调查数据中推断出相关的网络

摘要

基于高通量测序的技术，如16S rRNA基因分析，有可能阐明自然微生物群落的复杂内部运作--无论是来自世界海洋还是人类肠道。探索这些数据的一个关键步骤是识别这些群落成员之间的依赖关系，这通常是通过相关分析实现的。然而，自卡尔-皮尔逊时代起，人们就知道对这种技术产生的数据类型（被称为组成数据）的分析可能会产生不可靠的结果，因为观察到的数据采取的是基因或物种的相对部分的形式，而不是它们的绝对丰度。利用人类微生物组计划的模拟和真实数据，我们表明这种组成效应可能是广泛和严重的：在一些真实的数据集中，分类群之间的许多相关性可能是人为的，真正的相关性甚至可能以相反的符号出现。此外，我们表明群落的多样性是调节这种组成效应强度的关键因素，并开发了一种新的方法，称为SparCC（可在https:///yonatanf/sparcc），它能够从组成数据中估计出相关的数值。为了说明SparCC的潜在应用，我们推断了一个丰富的生态网络，连接了人体18个部位的数百个相互作用的物种。使用SparCC网络作为参考，我们估计标准方法对每个真实的相互作用产生3个虚假的物种-物种间的相互作用，并且错过了人类微生物组数据中60%的真实相互作用，而且，正如预测的那样，大多数错误的联系是在多样性最低的样本中发现。

引言

利用高通量基因组调查，如16S rRNA基因剖析，对自然界群落的研究已经成为常规，然而适当的、经过充分验证的分析方法的开发仍在进行。第一个挑战是从16S rRNA基因序列中获得可靠的和信息丰富的计数，通过过滤虚假的读数，并将剩余的读数以有意义的方式分组。一旦获得这种计数，就需要应用适合离散调查数据的分析技术。

基因组调查的一个共同目标是确定生态群落内分类群之间的相关性。相关分析为实现这一目标提供了一条很好的途径，但我们表明，当它应用于基因组调查数据（GSD）时是无效的，并可能产生误导性结果。与GSD相关的挑战来自于这样一个事实，即它们是对群落成分丰度的相对测量，而不是绝对测量。构成这些数据的计数（如16S rRNA基因读数）是由从群落中提取的遗传物质的数量或测序深度设定的，而分析通常是通过将观察到的计数归一化为总计数来开始。由此产生的分数属于一类被称为封闭式或组成式的数据，并构成其特殊的几何和统计属性。具体来说，从GSD计算相关性的标准方法在理论上是无效的。相关性估计是有偏差的，因为它们的总和必须是，分数不是独立的，并且倾向于有一个负相关，而不考虑基础绝对丰度（称为基础丰度）之间的真正相关性。因此，相关性的估计往往反映了数据的组成性质，而不是表明基本的生物过程。事实上，在1897年，卡尔-皮尔逊警告说，不要 "试图解释分子和分母包含共同部分的比率之间的相关性"，从那时起，已经证明许多其他标准分析技术在应用于此类组成数据时是无效的，而且它们的解释是不可靠的，往往会产生误导。尽管如此，这些方法仍然是微生物生态学研究中使用的主要工具。

尽管已经开发了组合数据分析的方法，但推断成分之间的依赖关系这一基本任务仍然是一个突出的挑战。一个广泛使用的方法是Aitchison的完全子成分独立性测试，它测试是否存在任何依赖关系，但不指出哪些成分是相关的，也不指出相关的大小。Filzmoser和Hron最近开发了一种方法，用于推断经过适当数学转换的成分数据的相关性，但他们的方法没有提供将转换后的变量的相关性与基础基因或物种的相关性相关的映射。

在本文中，我们首先使用来自人类微生物组计划（HMP）的模拟和真实世界的数据来证明GSD可以被 "组成 "效应严重偏离，然后确定调节其严重程度的因素。最后，我们提出了一种新的方法，称为SparCC，并表明即使在最具挑战性的数据集中，它也能以高精确度推断出相关关系。

结果

标准的相关推理技术在GSD上表现不佳

构成性假象在多大程度上影响了现实世界的GSD？我们将标准的统计方法应用于人类微生物组计划（HMP）的16S rRNA基因调查数据，这些数据测量了在个人不同身体部位发现的微生物群落的组成。每个群落的组成是以操作分类单位（OTU）来描述的。因为只有每个OTU的相对丰度是可用的，所以这些数据被认为是成分性的，因此会受到上述潜在偏差的影响。

从标准皮尔逊相关推断出的网络在不同的身体部位显示出不同的模式，表明了生物结构（图1，左栏。所有18个HMP身体部位见图S1）。具体来说，阴道中部、耳后皱襞和口腔粘膜网络的一个突出特点是存在一个与其他多个OTU呈负相关的OTU。尽管很想把生物学意义归于这些观察结果，但从具有类似分类群丰度的随机洗牌数据推断出的相关网络，但缺乏OTU之间的任何相关性（见材料和方法），再现了这一特征（图1，中间一栏），表明它可能来自封闭（规范化）过程。

这些虚假的相关关系背后的机制是直接的。在阴道中部网络中观察到的模式是由OTU （一种乳酸杆菌）的优势造成的。这个OTU的丰度为中位的，所以它的相对丰度的波动对群落其他部分的丰度有很大的影响，这仅仅是由于所有OTU的相对丰度之和为100%的要求：当乳酸菌的丰度变化时，所有其他OTU的相对丰度都以相反的方向一致变化，造成与乳酸菌的人为负相关，以及相互之间的人为正相关。

多样性和相关密度控制着组成效应（compositional effects）的严重性

组成效应在一些数据集中很严重，但在其他数据集中却很温和。我们发现，数据集中样本的多样性（通常被称为α多样性）是一个很好的预测成分效应强度的指标，它随着多样性的增加而减弱。直观地说，组成群落的OTU越少，组成效应就越差，极端的情况是群落只由两个OTU组成，而这两个OTU看起来总是完全负相关的。此外，即使在由多个OTU组成的群落中，如果只有几个OTU在群落中占主导地位，组成效应也会很显著。这种多样性的概念可以用香农有效OTU数来量化，它既可以量化OTU的数量，也可以量化群落中的主导地位。

具有已知相关性的模拟网络（见材料和方法）说明了多样性对成分假象的影响。真正的相关性（图2A-C）只有在群落具有多样性时才能恢复（图2F）。在与HMP样品相似的多样性网络中，推断出的联系往往被与优势OTU的负相关所支配，这导致了其余OTU之间的正相关（图2D,E）。这种影响是如此强烈，以至于它消除了OTU 4与OTU 3和5之间的负相关，以及OTU 1和2之间的正相关（图2E）。更糟糕的是，随着多样性的进一步减少，OTU 4与OTU 3和5之间的负相关变成了明显的正相关（图2D）。值得注意的是，这些组成效应并不局限于皮尔逊相关，也存在于非参数相关中，如Spearman相关（图S2）。

如果基础网络具有真正的正相关，那么组成效应甚至比基于群落多样性的预期更加明显。发生这种情况是因为成分之间的强相关性降低了样本的有效多样性（即两个完全相关的OTU表现为一个OTU）。这种影响可能会混淆天真的努力，即通过比较观察到的相关性和洗牌后的网络来纠正成分效应。当数据被洗牌时，如图1的中间一列，相对于未被洗牌的数据所观察到的结构，可能会出现一些虚假的连接（如对口腔粘膜样品的观察），从而对观察到的网络产生错误的信心。因此，随机化并不足以确定观察到的相关性的重要性，也不可能通过与随机网络的比较（或 "减去"）来确定相关性。

SparCC：从GSD推断相关关系的新程序

在这里，我们描述了一种从组成数据中推断相关关系的新技术，称为SparCC（组成数据的稀疏相关关系）。SparCC估计了对数转换后的成分之间的线性皮尔逊相关关系。由于这些相关性不能准确计算（如下所述），SparCC利用一种基于以下假设的近似方法。(i) 不同成分（如OTU或基因）的数量很大，以及(ii) 真正的相关网络是 "稀疏的"（即大多数成分之间没有强烈的相关性）。随后，我们表明SparCC对违反稀疏性假设的情况具有惊人的鲁棒性。SparCC不依赖于任何特定的基础变量分布，即群落中的真实丰度可以遵循任何分布，在随后的例子中选择对数正态分布的动机仅仅是为了便于实施和经验上的适合。为了清楚起见，我们以16S rRNA基因数据为背景介绍该方法，其中成分是OTU，基础变量是它们在群落中的真实丰度，但SparCC可以应用于任何其近似值有效的成分数据。

像大多数成分数据分析技术一样，SparCC是基于对数转换的：

这种转换有几个优点。首先，新的变量含有关于OTU真实丰度的信息，因为分数的比率等于真实丰度的比率。第二，与分数本身不同，两个OTU的分数之比与分析中包括的其他OTU无关，这一特性被称为亚组合一致性。第三，这种转换在数学上是方便的，因为新的变量不再局限于单数，而是可以自由地承担任何实值。取对数消除了正性约束，并在处理变量时诱发了（反）对称性。

为了描述构成性数据集中的依赖关系，艾奇逊建议使用以下数量

通过迭代上述程序，可以实现更准确的估计。在每次迭代中，前一次迭代中确定的相关性最强的OTU对被排除在基础方差估计之外。这就加强了剩余配对中的稀疏性，产生了更好的方差和相关估计。

需要从观察到的计数中估计出OTU的分数，以应用SparCC。用样本中的总计数对每个OTU进行归一化处理（最大似然估计）对稀有OTU来说是不可靠的，因为它高估了零分数的数量。这可能会产生由测序深度的变化所引起的伪影。这些假象促使一些作者对他们的数据进行降样，使所有的样品具有相同的总计数，然而降样并不能缓解成分效应，而且需要丢弃相当一部分的可用数据。因此，我们采用了贝叶斯方法来估计成分比例（见材料和方法），这样可以评估下游分析的稳健性和分配置信值。

SparCC在模拟数据上有很高的准确性

我们使用先前描述的模拟数据集来证明SparCC在推断相关性方面的准确性，即使是在由单一OTU主导的高度问题性的组成数据中（图2G-I）。通过创建不同多样性和密度的多个模拟数据集，对SparCC进行了更系统的评估。我们用OTU之间的平均Pearson相关性来衡量密度，这样，密度大的数据集有更多强相关的OTU，挑战了SparCC使用的稀疏性假设。对于密度和多样性的每个组合，都分配了多个真正的相关网络，并对相应的数据进行了采样。使用均方根误差（RMSE）评估了由SparCC或标准相关推断的网络（图3）。标准技术只对非常多样化、稀疏的网络给出了合理的估计（Pearson RMSE ），而对于具有与HMP组中观察到的多样性相当的网络，Pearson RMSE是不可接受的，达到了与阴道中部的多样性相似的群落。Spearman相关性的表现只是稍好一些（图S3A）。相比之下，SparCC的性能与多样性无关，并且对所有的参数值都有改善，甚至对违反稀疏性假设的密集网络也有改善。事实上，SparCC达到的最差的准确性（对于不真实的密集网络），与在高度多样化的样本上使用标准相关性所达到的最佳准确性相当。此外，尽管更强的相关性可以被更可靠地估计，但使用标准方法，需要将注意力限制在特别强的相关性上，然后才能显著提高准确率，所得到的准确率最多只能与SparCC的准确率相当（图S5）。

SparCC识别HMP数据中的系统发育结构的相关性

我们用SparCC从HMP数据集（图1，右栏，图4）和它们相应的洗牌数据集（其中所有的OTU是不相关的）中推断出分类群-分类群相互作用网络。与图1所示的天真方法相比，SparCC在洗牌后的数据集中没有发现明显的相关性（数据集S1）。然而，对于真实的数据，发现了许多相关关系，与标准的Pearson相关关系有很大的不同。SparCC推断表明，平均来说，使用Pearson识别的相关OTU对是错误的，而使用Pearson识别的相关OTU对是遗漏的（见表S1按身体部位的分类）。特别值得注意的是，我们观察到OTU 3和OTU 148之间的正相关，它们都属于乳酸菌属，而Pearson网络中却没有，这可能是由于高含量的OTU 3偏向于做负相关的关系。有趣的是，使用SparCC，我们观察到系统发育相关的类群之间出现正相关的可能性更大（表S2），这一发现从表面上看似乎支持中性群落动力学的作用，因为相关的生物可能居住在类似的壁龛中，但似乎并没有通过竞争性排斥而占主导地位（尽管更复杂的情况当然是可能的）。我们预计像SparCC这样的技术将在分析这些数据以解决这一问题和其他基本的生态学问题方面发挥重要作用。

讨论

在这项研究中，我们把重点放在组合数据分析的一个突出挑战上--推断相关性。我们已经证明，在人类微生物组的16S rRNA基因调查中，组成效应是明显的，并且在这种数据特性的激励下，我们开发了一种新的程序来估计相关关系。

我们发现，物种的多样性和相互作用的密度是影响成分效应对相关性估计的严重性的两个关键因素，其中低多样性、高密度的数据是使用标准方法推断相关性的最大挑战。SparCC不依赖于高多样性，相反，它只要求相关性的稀疏性，但在实践中，即使强烈违反了稀疏性假设（所有成分对中的%30是强相关的），也是稳健的。因此，我们建议在任何具有低多样性的GSD上使用SparCC：作为一个经验法则，我们建议标准技术的有效成分数至少为50个（可能需要注意的是，如果许多OTU之间存在强正相关，有效多样性可能比估计的要低很多）。我们强调，仅仅拥有许多成分并不足以避免成分效应。例如，来自HMP的16S rRNA基因调查包括数百至数千个不同的OTU，但其有效物种数相对较低，少数物种在大多数样品中占主导地位。

GSD的一个重要子类是使用DNA微阵列、RNA-seq和ChIP-seq等技术进行的全基因组调查。这些全基因组数据也会受到组成效应的影响，然而，由于这些数据往往具有高的多样性，它们可能要轻得多或可以忽略不计。例如，通过数据库提供的微阵列实验中的平均有效基因数是针对S. cerevisiae和E. coli的。这可能解释了为什么到目前为止，生物科学中对组成效应的关注相对少于其他学科。

零值的大量出现是GSD的另一个关注领域。这些零值既可以代表群落中真正不存在的成分，也可以代表从群落中抽取的样本中偶然不存在的稀有成分。如果没有额外的知识，这些选项是无法区分的，而且根据分析的目标，研究者必须决定如何解释它们，并相应地选择分析方法。我们强调，零值的处理是一个挑战，它绝不是成分数据所特有的，而只是因为分析这些数据所采用的对数转换而凸显出来[21]。在这项研究中，我们通过增加小的假数来消除零分，详见材料和方法。补充性的方法，即对零的处理与非零值的处理不同，其挑战性大大增加，也是正在进行的研究的主题[22]。

尽管本文提出的方法允许检测社区内的相关性，但仍然存在许多挑战。首先，SparCC依赖于可靠的成分计数，正如介绍中所指出的，这并不是一件小事。第二，SparCC估计的相关性是衡量对数转换丰度之间的线性关系。推断成分之间更普遍的依赖关系的组合方法，相当于非组合数据的等级相关和相互信息，还没有被开发。第三，将群落内检测到的模式与外部因素联系起来（例如，将人类肠道微生物群落的组成与人类健康状况联系起来），以及检测群落内和群落间的时间模式需要非标准的组成方法。虽然存在一些这样的方法，但它们很少在GSD的背景下使用，也没有针对其特殊的属性进行定制。最后，GSD通常与系统发育信息（物种或基因的相关度）相关联，理想情况下，这些信息会被纳入分析中（例如，加权的UniFrac距离，它试图捕捉群落丰度和系统发育组成的差异）。我们认为，开发系统的、具有统计学意义的方法来分析GSD的组成，是了解生物群落结构、群落演化过程以及形成群落的力量道路上的一个必要步骤，因此代表了未来研究的一个重要方向。