Nature子刊：微生物组数据转换以提高宏基因组定量准确度

文章导读：宏基因组研究日益广泛，但其定量分析一直面临很多困难。这篇文章系统的总结了宏基因组流程中影响定量分析的各个方面，尤其是数据的组合性以及样品微生物负荷的变化。这篇文章的亮点在于将对定量分析的干扰细化到不同生态场景的微生物群落，甚至是不同类群，帮助我们深入了解宏基因组数据结构，正确认识下游分析中数量关系的可靠性，避免在研究中做出错误的研究结论（而这些错误结论在以往研究中可能并不罕见）。

Benchmarking microbiome transformations favors experimental quantitative approaches to address compositionality and sampling depth biases

基准微生物组转换处理组成和采样深度偏好以支持实验定量

作者：Verónica Lloréns-Rico, Sara Vieira-Silva, PedroJ. Gonçalves, Gwen Falony & Jeroen Raes

期刊：Nature Communications (IF=14)

时间：11 June 2021

文章摘要

虽然宏基因组测序已经成为研究宿主相关微生物群落的首选工具，但由于序列矩阵的稀疏性和组合性，下游分析和微生物组数据的临床解读仍然具有挑战性。本文评估了目前提出的计算和实验的方法来减轻这些突出问题的影响。本研究使用虚拟微生物群落的方法产生粪便宏基因组数据，对13种常用的数据转换方法在多样性估计、物种与物种关联的识别以及物种和宏数据相关性的识别方面的性能进行了基准测试，以应对不同微生物生态系统负荷的挑战。我们发现通过实验方法将微生物负荷纳入下游分析的定量方法，比旨在减轻数据组合性和稀疏性的计算策略表现得更好，不仅改善了对真阳性相关的识别，还能减少假阳性检测。当分析虚拟情节中由于炎症病变导致的微生物负荷失调时，校正采样深度的定量方法比未校正的方法具有更高的精度。总的来说，我们的研究结果提倡在微生物组研究中更广泛地采用实验定量方法，但在无法确定样品微生物负荷也推荐了一些特定情形的首选数据转换方法。

研究背景

宏基因组测序允许以前所未有的规模和通量剖析微生物群落。由于避免了可培养偏好，测序方法有望对复杂的微生物生态系统进行真正的随机取样。持续进行的标准化工作，则可以最小化系统偏差，例如DNA提取和扩增带来的偏差。然而，即使在最小化技术偏差之后，对当前高标准的宏基因组的分析和解释仍面临严峻挑战。

当前宏基因组实践的第一个挑战来自于测序数据的成比例性质。在宏基因组测序流程中，为了确保最大的测序成功率而对测序文库的制备方案进行了优化。例如，序列片段的提纯步骤是为了确保最大的产量，使得甚至是很低浓度的细胞浓度或DNA量，也能获得可以实现测序的片段浓度。然而，即使是推进自动化和标准化之后，这些过程仍无法确保测序数据量和原始样品细胞浓度联系起来。因此，得到的序列数据矩阵只能分析微生物特征（分类单元或功能）的相对比例，而原始样本的微生物群落无法定量。在这样一个比例数据结构中，相对丰度并不是独立的（也即数据的组合性compositionality）单个分类单元或代谢途径丰度的相对变化会伴随等量的其余成分的变化（也即假如一个物种丰度增加，其他物种的丰度就会降低相同的量，每一个物种的丰度数据都不是独立的），从而诱导产生负相关偏差（假如群落中一个物种绝对数量增加了，尽管其余物种绝对数量并没有变化，反映在相对丰度上就会是减少的，就会产生莫须有的负相关）。这种限制影响所有的下游微生物组分析。比例分析不能表征样品间微生物组组成或代谢潜能差异的数量大小和方向性。此外，当对生态系统密度(微生物负荷，表1)展现显著变化的群落进行分析时，数据的组合性妨碍识别微生物特征与环境或临床协变量之间的相关模式（微生物负荷指的就是样品中的细胞浓度，如果微生物组细胞绝对数量变化较大，就很难根据相对丰度挖掘真实的相关模式，例如所有的物种都随着氮元素的增加而增加，但是其中一个物种爆发式增长，反映在相对丰度上其他物种的丰度反而是降低的！）。

第二个挑战来自当前宏基因组测序方法中低并且多变的采样深度。在微生物组研究中，采样深度定义为测序的细胞（观察到的细胞大小）与样品中总的细胞（真实群落大小）的比例。与之表面相似但实际不同的是测序深度，指的是单个样品产生的测序数据总量（表1）。随着技术的发展与测序成本的降低，这些年测序深度逐渐增加（也即科学家们测得数据量越来越多），然而宏基因组分析仍然具有浅采样深度的特点（尽管测的越来越多但是还是远远不够）。最近一项研究对876个粪便样品进行鸟枪法宏基因组测序，测序深度为5.5-18.2 Gb，而平均的采样深度只有0.0045%（对于粪便这种高微生物负荷的样品，即使测20G深度仍然很低，对真实的微生物群落探索好比管中窥豹）。在评估细胞密度高（微生物负荷高）、多样化和不均匀分布的微生物生态系统时（如粪便），较浅的取样深度会导致产生的微生物群落特征矩阵具有稀疏性，无法区分物种或功能是真的缺失还是未检测到（矩阵的稀疏性是指观察到的物种或功能仅仅是真实群落的很小一部分）。采样深度带来的样本间差异使问题更加复杂。理想情况下，即使是在分析比例组成时，不同样品中微生物群落的抽样比例应该保持恒定。当评估的样本具有不同的微生物负荷，保持或创建一个独立于样本密度的均匀测序深度，就会导致以不均匀的采样深度产生序列矩阵（样本微生物负荷不同，测相同或随机的数据量，必然会导致采样深度不同甚至相差悬殊）。在上述讨论的876次粪便宏基因组调查中，实际采样深度的差异超过40倍。这种变异意味着，特定微生物组特征只在一部分样本子集中检测到可能是由于采样深度的不受控制的变化造成的。因此，尽管宏基因组分析不能确定某物种或功能真实缺失是常识，但很少有人认识到，他们可能还会将单纯的人为技术结果识别为存在。

表1. 术语表

由于日益认识到组合性和抽样不足对宏基因组分析的严重影响，导致潜在的环节策略广泛发展，包括计算性和实验性的策略。微生物组研究人员首先着手处理不同测序深度相关的问题（通常被错误地认为等同于采样深度，即假设采样的群落中微生物密度相同）。面对技术的变化以及多年来单样本测序数据量的总体增加，研究人员建议在数据集内和跨数据集的比较中进行稀疏化(或下采样，也即根据最小测序数据量进行抽平)从而处理不同样本测序数据深度不同的问题（这是对样本微生物负荷不了解下的一种妥协，也即我假设所有样本微生物负荷相同，对测序深度进行抽平，采样深度也就相同了）。

与直接按比例(相对丰度)标准化相比，将数据随机抽取子集从而均匀化测序深度，允许在不同样本间观察到的丰度（也即抽平后观察到的数目）进行样本比较，而不依赖于生成的原始序列数量。然而，基于测序深度的简易处理很快就受到了批评，不仅因为浪费和丢弃了低丰度分类单元的信息，还在于不适合微生物负荷差距显著的样本。作为响应，替代的计算方法被提出，包括一些稳健数据转换方法的使用，可以同时处理不同的测序深度干扰以及组和数据分析的限制。这些方法取自组和数据的Aitchison统计，以及基于RNA测序的转录组测序中发展的标准化和/或差异特征丰度测定方法。虽然有人建议这种转换能够在不需要缩小数据规模的情况下处理不均匀的测序深度，但尚不清楚它们在多大程度上依赖于所有分析的样本微生物负荷相同的假设。

虽然减轻组合性影响的计算方法的应用将限制人为结果(如由相对丰度引起的负相关)被识别为真正的相关，但是也不会重新获得那些已经丢失的联系微生物负荷和测序数据量的信息（也即这些方法只是为了避免错误结果，但并不能真正解决问题，丢失的信息无法复现）。为了解决这个问题，最近发展了一些实验方法来保存或重新获取这些与微生物密度的关联信息，包括DNA或细胞的spike-in对照（通过实验将人工合成的DNA标记序列添加入特定细胞的基因组，然后根据测序数据中这些序列的数量对微生物群落进行绝对定量）、使用定量PCR对测序数据进行平行化以及流式细胞仪计数微生物细胞。通过将比例转换为计数，这些方法剔除了下游分析中组合性数据的限制。除了实验方法上的差异，这些定量方法在将获得的微生物负荷纳入下游分析的方式上也有所不同。有两种方法可以区分。其中一种，绝对计数缩放程序将相对序列矩阵与实验确定的微生物负荷直接相乘，保留产生的所有测序信息。而另一种，则使用严格的缩减步骤，通过对测序数据随机抽取子集的方法来平均化采样深度，使其序列/负荷比达到最优水平，会丢掉一些序列信息以及采样深度不足的样本（可以看作特殊的抽平，不是均匀化测序深度，而是结合实验测定的微生物负荷均匀化采样深度）。这样做的目的是避免对低微生物负荷下低丰度物种的过度检测，这些低丰度物种可能干扰下游分析中关于多样性和关联的分析（也就是减少无用信息，只想获得微生物群落的主体结构，因为一些很稀有的微生物很可能是没有活性的，采样深度过大导致稀有物种数目很多会导致群落的多样性指数非常突兀）。还有一种可选情况，产生的非组合(不均匀)序列矩阵可以按比例缩放到绝对计数，以获得一定单位的数量(例如每克样品的细胞数)。

最近的一些研究使用真实数据和模拟数据比较了相对标准化和处理组合性的转换对微生物组数据集相关结构的影响。在检测物种之间相关性方面定量方法被认为优于相对标准化的分析，但迄今为止还没有针对组合性转换的系统基准。此外，由于目前的定量转换需要对微生物负荷进行实验测定，它们不适合所有类型的样品收集方法(例如，使用不记录重量的稳定缓冲液)，并这些方法是过度劳动密集型的任务。

本研究中，作者对广泛的可用的计算和实验转换方法的优点和局限性进行了系统的评估，这些方法已经被提出来处理序列数据分析中的组成性和采样深度变化。我们的基准测试表明，与其他计算方法相比，定量方法在报告样本丰富度、准确恢复物种和物种关联方面的性能得到了改进，同时最大限度地减少了假阳性关联的检测。

主要结果

本研究的工作流程：第一步，使用Matlab产生虚拟的微生物群落数据和宏数据（metadata, 图1a），这些数据中存在一定的关联结构，包括物种和物种（taxon-taxon）之间的关联以及物种和宏数据（taxon-metadata）之间的关联。第二步，使用计算机模拟细胞计数和宏基因组测序过程，获得微生物负荷以及测序深度数据（图1a）。第三步，使用13种常见的数据转换方法对上一步产生的数据矩阵进行处理（图1a和表2，大家注意，本文的数据转换方法看似有些拗口，实际上在methods中可以看出这些都是很常见的数据处理方法，只不过我们平时的叫法不同）。第四步，通过对比转换后的测序数据中能否检测到原始数据集中的关联结构来评估不同转换方法，也即能成功识别第一步数据集中物种和物种以及物种和宏数据之间的相关性则为真阳性，而原始数据没有后续人为产生的相关则为假阳性。

图1. 本研究工作流程

此研究亮点在于，还虚拟了三种常见的生态学场景，因为不同场景的微生物群落可能受数据转换的影响不同，这三种场景分别为连续演替（succession ）、爆发（blooming）和失调（dysbiosis）。如图1b所示，连续演替过程中不同样本细胞数目连续变化，且不同分类单元的比例变化不大，群落结构比较稳定；而爆发场景中某一个分类单元大爆发，可以参考蓝藻爆发或者某种病原菌爆发，这时候群落结构被某一物种主导，由于其比例突然增大，在基于相对丰度的分析中其他物种容易出现假阳性的负相关。失调场景类似于环境突然变化或者病人突然的炎症反应，导致某些分类单元的丰度出现突然的增长或减少，在这种场景中，变化较大的一般是对环境敏感的物种，对环境迟钝的物种则变化缓慢，这在下文也有讨论。

表2. 宏基因组数据转换基准

宏基因组分析的首要问题就是我们获取的数据能在多大程度上反应真实微生物群落中的物种多样性。因此，作者首先评估了不同转换方法对群落物种丰富度恢复的影响（基于样品的alpha多样性指数）。通过分析原始数据集和转换后数据集中相同样本丰富度的相关性可以发现，演替和失调场景相关性均很显著，不同转换方法区别不大；而爆发场景除了QMP（也即实验定量方法）外，效果均很不理想（图2a，b）。这是可以理解的，因为爆发场景除了组合性的严重干扰外，由于单一分类单元丰度过高很容易导致采样深度不足，其他物种难以检测。这给我们的经验是，面对爆发场景的样本，最好结合实验进行定量分析，或者使用Shannon或Simpson等结构指数。

图2. 不同转换方法在评估丰富度时的性能

根据以往研究结果，在产生的原始数据集中物种（分类单元，taxon）大部分都被强迫与负荷呈现正相关，这符合生态学研究实践以及我们的常识，也即微生物密度越高，大多数物种的绝对丰度也会越高。而在转换后的测序数据中，除QMP和ACS外，很多这种正相关难以成功恢复，尤其是在爆发场景中（图3a，b）。不过爆发场景中的爆发者和失调场景中的机会主义者（环境响应迅速）的相关性很容易被恢复，而失调场景中的反应迟钝者的相关性则全部丢失（图3c）。这说明，微生物群落里相对丰度变化大的结构主导者其与负荷的相关性更加稳健。

图3. 不同转换方法在恢复物种与负荷关联时的性能

最后，我们比较关注在现实中那些微生物与微生物以及微生物和环境之间的相关性能否通过宏基因组测序被检测到。结果发现，在物种和宏数据关联的恢复上，只有QMP和ACS性能不错，其余方法都比较差，尤其是在爆发群落中，检测到大量假阳性相关（图4a，b）。这警示我们，在得出物种与环境的关联结论时需要格外小心。从不同类群上来说，爆发场景中的爆发者与宏数据的相关更加稳健，可以比较容易的恢复（图4c）。而物种之间关联的恢复效果则比较乐观，假阳性较少（图4d）。

图4. 不同转换方法在恢复物种和物种、物种和宏数据关联时的性能

在实际研究中，我们要根据所面临的生态学场景以及研究目的确定合适的研究方法。作者模拟了一个患病情况下低负荷的失调场景的数据集，可以想象成肠道炎症腹泻脱水患者的粪便（图5a）。在这种情况下即使较低的测序深度也可以很好的恢复物种与疾病的关联，因此即使根据采样深度的缩减也不影响下游分析，反而会提高分析效率（图5b）。ACS在检测类群与疾病协变量的真实关联方面显示出更高的敏感性，不过ACS和QMP之间的敏感性差异随着患者和对照组之间负荷差异的增加而减小（图5c）。而QMP的总体精度更高，并随着测序深度的增加而扩大（图5d）。

图5. 用定量方法检测疾病相关分类单元