微生物绝对丰度变化新方法(QMD)

文章题目：QMD: A new method to quantify microbial absolute abundance differences between groups
发表杂志：iMeta
发表时间：2023年1月5日
通讯作者：刘星吟
原文链接：https://onlinelibrary.wiley.com/doi/10.1002/imt2.78

摘要

组间差异菌识别(differential abundant taxa， DA)是揭示不同条件下微生物群变化的主要微生物组学分析方法之一。已有报道称，使用微生物相对丰度数据开展DA分析获得的结果与使用微生物绝对丰度(即单位重量或单位体积下样本包含的微生物真实丰度)进行的分析结果相比存在偏差，而这种偏差则会导致假阳性率过高以及对DA结果的错误解读。本研究提出了一种基于微生物相对丰度来估计不同条件下各类群中微生物绝对丰度变化的新方法(quantification of microbial absolute abundance differences，QMD)。在构建的QMD模型中，组间微生物绝对丰度的差异为其相对丰度的差异及偏差之和。正确性验证和仿真实验表明，QMD可以鲁棒地处理小样本量、大比例差异丰度类群以及多样化丰度变化尺度。与其他方法相比，QMD能帮助我们更好地理解组间微生物的变化。此外，我们还提供了QMD软件，以辅助开展微生物群组件变化的研究。

引言

测序技术的发展使我们能够深入了解微生物群的组成和复杂微生物群的动态。尽管微生物绝对丰度(即单位重量或单位体积下样本包含的微生物真实载量)是代表微生物群组成的最真实数据，但大量的分析依然使用相对丰度。作为微生物组学分析主要方法之一，DA有助于揭示不同条件下微生物群的变化。学界已经表达了对使用相对丰度数据或测序原始计数数据进行DA分析的关注。由于缺乏采样和扩增的标准化，原始计数数据通常无法直接用来描述微生物丰度的变化。另一方面，某个微生物绝对丰度的增加或其他微生物绝对丰度的减少均会引起相对丰度的增加，所以使用相对丰度解释菌群的组间差异也很困难。

菌群相对丰度的组间变化与绝对丰度的组间变化之间存在偏差。这一偏差可能会导致DA中的错误发现率(FDR)升高。一些DA方法已经在尝试降低偏差对DA的影响。

ANCOM（微生物群组成分析）是一种广泛使用的DA工具。ANCOM方法通过分析微生物丰度与所有剩余的其他微生物丰度的对数比来进行组间差异微生物的识别，其已集成到QIIME 2 中。丰度比的引入巧妙地避免了相对丰度数据和绝对微生物丰度之间的差距。Weiss等人报道，在Mann-Whitney U检验、DESeq 、DESeq2 、edgeR 、Voom 和宏基因组测序等方法中，ANCOM在FDR控制上表现最好。然而，ANCOM无法量化微生物绝对丰度的组间差异。

DR（差异排序）表明，微生物组间相对丰度变化的大小排序与绝对丰度变化的大小排序相同。DR通过从相对丰度数据构建多项式回归模型来量化绝对微生物丰度的变化。DR的缺点是其未提供筛选出显著变化差异菌的方法。DR建议差异菌群的分析可以重点关注相对丰度变化较大的微生物。

ANCOM-BC(带有偏差修正的微生物群组成分析)使用线性回归框架中的E-M算法来估计抽样分数（sampling fraction），即样本中微生物的预期原始计数与其在生态系统中的绝对丰度之比。当样本量很小时，ANCOM-BC的性能较弱。

在本研究中，我们提出了一种基于微生物相对丰度来估计不同条件下各类群中微生物绝对丰度变化的新方法(quantification of microbial absolute abundance differences，QMD)。同时我们构造了基于QMD开展DA的统计检验假设模型。分析表明，微生物绝对丰度组间差异可以用其相对丰度加上组间的总微生物丰度的变化来估计。我们使用一元线性规划模型对总微生物丰度变化(即偏差)进行估计。QMD方法得到了真实的微生物绝对丰度变化数据的验证。与其他方法相比，QMD在微生物差异丰度定量、DA识别、运行速度上都表现出较好的性能。

结果

微生物差异丰度的量化

前期研究发现，微生物绝对丰度组间差异与其相对丰度的组间差异不一致。例如，与对照组相比，生酮饮食组小鼠乳酸杆菌绝对丰度降低(图1A和B)，但是其相对丰度增加(图1C和D)。

图1. 微生物总丰度变化带来的组间差异菌群识别的偏差

(A)-(B)，分别为10天正常饮食组小鼠和生酮饮食组小鼠的小肠菌群绝对丰度。(C)-(D)，分别为10天正常饮食组小鼠和生酮饮食组小鼠的小肠菌群相对丰度。结果显示生酮饮食组小鼠乳杆菌的绝对丰度下降，而相对丰度上升。(E) 相对丰度组间差异与绝对丰度组间差异的关系。

是什么导致了这种偏差呢？我们研究了微生物绝对丰度组间差异和相对丰度组间差异之间的联系。微生物Taxai的载量与其绝对丰度之间的关系可表述为

其中a_{i,j}^K代表K组、样本j中Taxai的绝对丰度。m_j^K代表K组、样本j的取样质量。A_{i,j}^K代表K组、样本j中Taxai的微生物载量。K代表组，为便于说明，使用 C(对照组)和T(处理组)表示K。

因此，K组(样本j)中总微生物载量可计算如下:

其中s_j^K为K组(样本j)的抽样分数，即样本中微生物的预期原始计数与其在生态系统中的绝对丰度之比。抽样分数的概念继承自Shyamal Das Peddada团队，覆盖了采样过程、扩增和测序操作对绝对丰度的影响。根据他们提出的假设，抽样分数是常数，根据抽样分数的定义，我们有

根据原始计数数据，K组样本j中的Taxai的相对丰度o_{i,j}^K计算如下

根据1~4，相对丰度和绝对丰度之间的关系可以表示为

即

其中(Λ_j^K)/(m_j^K)正好是所有微生物绝对丰度之和。使用Ψ_j^K表示(Λ_j^K)/(m_j^K)，我们有

公式7表明，绝对丰度与相对丰度数据呈线性相关。文献已经证明了公式7的正确性。相对丰度数据乘以总微生物丰度即菌群的绝对丰度。利用这种方式推断的微生物绝对丰度与其实际的绝对丰度呈现强正相关。

我们更关注微生物丰度的组间变化。如前所述，各个实验组常由若干个样本构成。样本量从数十到数百不等。我们使用绝对丰度组间对数差异来衡量微生物绝对丰度的变化。记MiK为K组中Taxai的对数绝对丰度的期望。

因此，Taxai的绝对丰度组间差异为

注：这里使用以2为底的对数进行组间差异的计算。

因此，我们有

将∆Ψ表示为所有微生物丰度总的组间倍数变化

我们有

进一步的，∆Ψ作为公式9的常数部分，不同微生物偏差∆Ψ是相同。因此，要获得各个微生物绝对丰度的组间差异，只需要量化微生物总丰度变化∆Ψ，公式9中相对丰度组间变化则可以直接利用相对丰度数据进行计算。为了量化∆Ψ，我们假设:有相当一部分微生物丰度在组间的变化并不大。这个假设可以转化为即最小化所有微生物丰度组间变化的绝对值之和，即∑(Taxa_i)|∆M_i|。考虑到不同微生物的检出率不同，将Taxai的平均检出率(Di^T+Di^C)/2作为微生物丰度组间变化绝对值的权重。检出率指微生物检出的样本占总样本量的比率。例如，我们得到100个样本，在50个样本中检测到某微生物，则其相应的检测率为50%。由此，QMD模型通过寻找∆Ψ，使所有微生物丰度变化的加权绝对值最小，即见公式10。

其中Di^C和Di^T分别是对照组和实验组中的Taxai的检出率。

为了在公式10中找到最优点，QMD模型中，以0.01步进从-10到10进行遍历。该遍历范围涵盖了与对照组相比微生物丰度增加近1000倍或减少至千分之一的情况。我们提供的QMD软件提供了修改遍历范围和步骤的功能。将∆Ψ定量后，可根据公式9计算出各Taxai的绝对丰度组间差异∆M_i。

基于QMD的组间差异菌群识别

差异菌识别方法，例如t检验和Mann-Whitney U检验，其检验能力有限。公式9表明∆Ψ的存在使传统方法FDR过高。例如，如果∆Ψ为正，假设实验组中的总微生物丰度增加两番，那么相对丰度组间差异将小于绝对丰度的组间差异。这种情况下，如果相对丰度没有变化，微生物绝对丰度会增加两番，这意味着处理组中的微生物绝对丰度比对照组中的高四倍。负的∆Ψ则效果相反。

于是，我们构建了一个新的统计检验假设：

H0：微生物绝对丰度组间无变化
即：

公式12提出了基于QMD模型的组间差异菌群识别方法。在量化∆Ψ后，采用Mann-Whitney U检验进行统计学检验。为了便于QMD的应用，我们提供了QMD软件https://github.com/Xingyinliu-Lab/QMD/tree/master/GUI_QMD/release。我们还在https://www.youtube.com/watch?v=LPti1vgoiCo上提供了视频教程。QMD软件支持Benjamini/Hochberg FDR校正。

QMD正确性验证

近年来，一些研究开始关注微生物绝对丰度的测量，出现了采用数字PCR、spike-in 、流式细胞术、定量PCR 、hamPCR 等实验方案，支持微生物绝对丰度的测量。我们从上述文献中选择了LSI, STOOL, B1和B2四个实例，用来验证QMD方法的正确性。表S1中给出了这四个实例的详细信息。

根据实验获得的微生物绝对丰度，可以计算出微生物总丰度变化∆Ψ(表1)。显然，∆Ψ可以很好的描述相对丰度的组间差异和绝对丰度组间差异之间存在的偏移(图2A，图S1A)。该发现验证了公式9。

表1. 正确性验证的结果总结

注：MWU表示直接对相对丰度进行Mann-Whiney U检验，以鉴定差异丰度的菌群。

图2. 使用LSI数据集进行QMD验证

(A)验证相对丰度组间变化与绝对丰度组间变化之间的关系。开展实验观察到的真实总微生物丰度变化，以及QMD模型推断得到的总微生物丰度变化均可以校正偏差。(B)可视化LSI中的QMD对最优化过程，以找到总微生物丰度变化(∆Ψ) ̂。(C) 比较QMD、DR和ANCOM-BC绝对丰度组间差异定量结果。

QMD使用优化遍历对总微生物丰度的变化(∆Ψ) ̂进行推断。如图2B和S1B所示，QMD的目标函数是凸函数，没有局部最优陷阱。推断的(∆Ψ) ̂接近观测到的真实∆Ψ，尤其是在LSI和B2的情况下(表1)。推断的(∆Ψ) ̂加上相对丰度变化接近观测到的绝对丰度组间差异(图2A、图S1A)。该结果验证了公式10中给出的QMD模型。此外，我们还比较了QMD、ANCOM-BC和DR在微生物绝对丰度差异定量结果。在这四种情况下，QMD的表现都优于其他方法(图2C、图S1C、图S2A和B，表1)。

为了验证QMD在组间差异菌群识别中的作用，我们计算了四种实例下QMD、ANCOM和Mann-Whitney U检验的错误发现率（FDR）。QMD在LSI、STOOL和B2实例上的表现超过了其他方法。Mann-Whitney U检验则在B1实例中表现最佳(表1)。

接下来，我们比较了QMD、ANCOM和Mann-Whitney U检验在人肠道微生物数据集上的性能(表S2)。微生物总丰度变化表明不同疾病患者的肠道微生物丰度变化存在不同模式(图S3A)。在结直肠癌(CRC)患者中可以发现总微生物丰度下降的趋势(图S3A)。这种模式可能与不同疾病患者中不同的宿主-微生物群相互作用类型有关。抗生素鸡尾酒疗法等抗生素的使用也可能影响总微生物丰度的变化。例如，不同炎症性肠病(IBD)患者的微生物总丰度显示出相对较大的波动。这提示采用QMD开展多中心数据时应当谨慎，充分考虑人群的异质性、抗生素使用史、性别、年龄和其他微生物群相关因素。

另一个发现是对于丰度变化较小的微生物，Mann-Whitney U检验可能会给出相反识别结果(图S3B)。对于丰度变化较大的差异微生物识别，ANCOM可能产生假阴性结果。一个可能的原因是ANCOM关注的是微生物之间的丰度比。在大多数情况下，当多个微生物绝对丰度同时发生变化，这些微生物的相互丰度比值可能不会像微生物丰度本身那样发生很大的变化。

QMD在微生物差异丰度定量中表现稳健

我们构建了三个系列的数据集: H2029, Obesity, 和GP (Global Patterns)进行仿真实验(详情参见方法) ，以此研究样本量、差异微生物比例以及微生物丰度变化数值大小对QMD性能的影响。

对QMD估算的总微生物丰度变化(∆Ψ) ̂和∆Ψ真值进行线性回归。在所有三个模拟序列中(图3A、图S4A、图S5A、表2)其回归系数R-squared(adj)均大于0.97。这意味着QMD可以较好的完成∆Ψ估计。

表2. QMD定量的总微生物丰度变化及其真值间的线性回归

图3. H2029系列仿真实验结果

(A) QMD定量的微生物总丰度变化与其真值呈线性相关，表明QMD是一种合格的微生物总丰度变化估计工具。(B) ~ (C)，差异微生物比例和样本量影响QMD、DR和ANCOM-BC的MAE。(D) QMD、FDR调整的QMD、ANCOM、ANCOM- BC方法进行差异菌群识别的FNR。(E)-(F)，差异微生物比例和样本量影响各方法的FNR。(G)，QMD、FDR调整的QMD、ANCOM、ANCOM- BC方法进行差异菌群识别的FPR。(H)-(I)，差异微生物比例和样本量影响各方法的FPR。

我们比较了ANCOM-BC、DR、QMD和RAC这些方法开展组间丰度差异量化的结果。RAC方法是指直接将组间相对丰度差异记做为绝对丰度差异。QMD的平均绝对误差中位数为0.69（H2029数据集）、0.77（Obesity数据集）、0.01（GP数据集）。显然，QMD在量化丰度变化方面优于其他三种方法(见图3B和C、图S4B和C、图S5B和C、表3) 。

表3. ANCOM-BC、DR、QMD和RAC量化绝对丰度差异的MAE中位值

QMD建立在一部分微生物绝对丰度变化相对较小的假设基础上。正如预期的那样，本研究中发生丰度变化的微生物的比例影响了估计的准确性。MAE随比例增加而增大。与其他方法相比，QMD在处理大比例微生物丰度变化时更为稳健。在丰度变化的微生物比例小于90%时，QMD的中位MAE保持在1以下。当丰度变化的微生物比例低于20%时，DR达到较小的MAE中位值(低于1)。当丰度变化的微生物比例低于10%时，DR的表现优于QMD(图3B、图S4B、图S5B)。

相反的，大样本量可以提高统计推断能力。除RAC外，其余四种方法的MAE随样本量增加而略有下降。例如，对于样本量小于10的实例，QMD的中位MAE为0.98。样本量为90–100时，该值提高至0.64(图3C、图S4C、图5C)。

QMD在组间差异微生物识别中很好地平衡了FNR和FPR

我们使用FPR(假阳性率)和FNR(假阴性率)评估和比较了FDR调整的QMD、QMD、ANCOM、ANCOM-BC和MWU等方法的组间差异微生物识别能力。MWU方法将Mann-Whitney U检验直接应用于相对丰度比较。这里我们并没有比较DR方法，因为DR方法未能提供差异微生物识别功能。

仿真结果显示ANCOM-BC的FNR严重超高，MWU的FPR严重超高。这两种方法在FNR和FPR控制之间没有保持很好的平衡。QMD和FDR校正的QMD在FNR和FPR上均优于ANCOM。ANCOM中FPR和FNR的中位数较高(图3D和G，图S4D和G，图S5D和G，表4和5)。

FDR校正的QMD在FPR的表现优于单独使用QMD(图3H和I，图S4H和I，图S5H和I，表4)。在微生物丰度变化比例较小或小样本量的情况下，FDR校正的QMD其FNR更大(图3E和F，图S4E和F，图S5E和F，表5)。综上所述，当样本量大于50时，我们建议对QMD进行FDR调整。如果样本量较小，则不需要进行FDR调整。

和MWU的FPR中位值

和MWU的FNR中位值

讨论

微生物绝对丰度的差异识别可以帮助我们揭示组间微生物的变化，已经成为微生物组学分析的重要任务。相对丰度组间差异与绝对丰度组间差异之间的偏差可能会在差异微生物识别时产生误导。在本研究中，我们开发并验证了一种新的方法QMD来量化微生物的绝对丰度组间差异。

我们的分析显示，相对丰度组间差异与绝对风度组间差异之间的偏差等于总微生物绝对丰度变化∆Ψ。在推断量化这个偏差之后，QMD通过相对丰度组间差异加上该偏差来推断绝对丰度组间差异。QMD只处理一个未知变量：总的微生物丰度变化∆Ψ。这大大降低了QMD模型的复杂性，并使其分析与以前的方法(如DR，一次性估计数百个微生物绝对丰度的未知变化)相比更为稳健。本研究进行的仿真实验验证了这一结论。QMD比DR和ANCOM-BC更准确地估算了微生物绝对丰度差异。随着差异微生物所占比例的增加，DR的表现迅速恶化。如作者所述，ANCOM-BC方法则在小样本量的情况下表现较差。

由QMD估算的微生物丰度变化为组间微生物群动态变化提供了定量描述。我们提供了p值和经FDR调整的q值，以识别具有统计学意义的组间差异微生物。仿真实验结果表明，即使在小样本或差异微生物所占比例较大的情况下，QMD也能很好地控制FNR和FPR。

在我们的研究中，我们发现微生物相对丰度变化与绝对丰度变化呈线性相关，相对丰度变化数值最大的微生物与绝对丰度变化最大的微生物相同。DR得出的结论与之类似，微生物相对丰度变化排序与绝对丰度变化排序相同。在这一点上，如果研究人员仅关注丰度变化最大的微生物，那么基于相对丰度的DA识别结果与QMD相似。另一方面，对于其他的微生物，则需要考虑微生物总丰度的变化。如图S3所示，微生物总丰度的变化可能会带来相反的变化趋势。在DA鉴定时，还应考虑微生物群落丰度的零膨胀和过度分散问题。QMD去除了微生物零计数，在方法构建时不使用伪计数，避免了潜在的偏倚结果。整合其他过度分散相关模型可能是未来QMD发展的方向之一。

鉴于相对丰度数据的数学特性，如果所有微生物的丰度同时发生下降或上升，QMD可能会发生功能退化。在极端情况下，例如使用广谱抗生素治疗造成所有类群的微生物绝对丰度普遍下降。此时，公式9将退化为∆Mi=∆Ψ，组间的相对丰度不会发生变化。

方法

我们在McMurdie等人的方案基础上，构建了仿真框架 。图S6描述了基于真实数据变换策略的仿真框架。选取了Global Patterns (GP), H2029, 和Obesity, 3个数据集作为原始池。选择GP数据集是因为已经有若干项研究使用了GP数据集。另外两个数据库分别包含来自GMrepo的20 ~ 29岁健康人群和肥胖人群的肠道微生物组数据。

由于a_{i,j}K=Ψ_{jK}*o_{i,j}^K为每个样本指定均匀分布的ΨjK，以帮助生成绝对丰度数据。通过对人工绝对丰度数据池的重采样，获得了初步的仿真实例。重抽样样本量服从均匀分布U(6,100)。

对于每个仿真实例，按照均匀分布U(5%,95%)随机设定差异微生物所占比例。然后，根据比例随机选择微生物子集。实验组的微生物丰度变化与该子集相对应。微生物绝对丰度变化服从均匀分布U(-10,10)。

由此，我们获得了GP、H2029和Obesity系列共1500个仿真实例，每个实例500个。仿真条件包括样本大小U(6,100)、差异微生物所占比例U(5%,95%)和绝对丰度变化U(-10,10)的随机组合(图S7)。

在性能评估阶段，主要考察两类指标。在绝对丰度组间差异定量部分，采用平均绝对误差(Mean absolute Error, MAE)评估ANCOM-BC、DR、QMD和RAC的定量准确性。计算公式如下：

其中(∆Mi) ̂为定量的Taxai绝对丰度差异。

采用FNR和FPR对QMD、FDR调整的QMD、ANCOM、ANCOM-BC和MWU的差异丰度类群识别能力进行了评价。仿真实验重扰动的微生物被认为是差异微生物。未被扰动并被鉴定为DA的微生物被认为是假阳性。被扰动但未被鉴定为DA的微生物为假阴性。下面的决策表用来计算FNR和FPR。QMD模型的阳性、QMD模型的阴性数据是通过统计检验从QMD模型得出的，而真实数据的阳性、真实数据的阴性是从计算FNR和FPR的仿真数据集中得出的(表6)。

表6. 计算FNR和FPR的决策表

为了验证QMD模型正确性，我们选择了通过实验测量微生物绝对丰度的研究所提供的数据集。Barlow等提供了每个微生物的绝对丰度数据。根据绝对丰度数据，按照公式7计算相对丰度数据。Vandeputte等提供了通过流式细胞术计数的每个样本的总类群丰度。Vandeputte等的相对丰度来源于16s rRNA测序，微生物的绝对丰度由总丰度乘以其相对丰度得出。结合绝对和相对丰度数据，计算II型误差(FNR)和I型误差(FPR)对QMD模型进行验证。

在仿真中，ANCOM-BC(https://github.com/FrederickHuangLin/ANCOM-BC-Code-Archive/blob/master/scripts/ancom_bc.R), ANCOM(https://github.com/mortonjt/ancomP), DR( DR 算法提取自https://github.com/knightlab-analyses/reference-frames/blob/master/ipynb/simulation-benchmark.ipynb)在默认参数集下运行。在ANCOM-BC中，程序丢弃了零比例> = 0.9的微生物，并且使用Bonferroni pvalue adj方法，设置tol.EM = 1e-5，max.iterNum = 100， perNum = 1000， alpha = 0.05。在ANCOM分析中设置permutations=1000。在DR中，batch_size=3，learning_rate=1e-3，beta_scale=1。

据报道，伪计数可能会导致结果的偏倚。在QMD算法中并未使用伪计数，因为我们不能确定零计数是由于测序深度不足或只是该微生物不能定植。QMD提供了两种方法对筛选出的微生物进行处理，一种是直接将这些微生物排除在下一步分析之外，第二种方法是将这些微生物的丰度进行求和，并作为“其他”参与下一步分析。在仿真中，QMD丢弃少于5个样本中检出的微生物，使用Benjamini/Hochberg进行独立检验进行FDR校正，permutations设置为500。MWU方法将Mann-Whitney U检验直接应用于相对丰度。RAC方法是直接使用相对丰度组间差异作为绝对丰度差异。对于所有方法，用0.05来筛选有显著变化的微生物。

16S rRNA序列分析

为了验证QMD的正确性，我们从ENA(European Nucleotide Archive)获得了B1和B2的原始16S rRNA测序数据（PRJEB21504和ERP023761）。使用QIIME 2创建特征表，使用SILVA朴素贝叶斯分类器对ASV进行分类。

为了比较不同疾病间微生物总丰度的变化规律，从GMrepo中获得了PRJEB7949、PRJNA368966、PRJNA385949、PRJNA388210、PRJNA389280、PRJEB6070、PRJNA445640的属水平丰度。

对于PRJNA763023、PRJNA433269，通过vsearch v2.17.1_linux_x86_64提取特征表。

参考
https://mp.weixin.qq.com/s/W26ELwEFQzaw0XoaZzjBPQ