一种基于树的模型，用于解决微生物组组成计数数据中的稀疏性和分类单元协方差问题

《Statistics in Medicine》：A Tree-Based Model for Addressing Sparsity and Taxa Covariance in Microbiome Compositional Count Data

　　**摘要**

微生物组的组成数据通常是高维的、稀疏的，并且表现出普遍的跨样本异质性。我们提出了“逻辑树正态”（LTN）模型，这是一种生成模型，它允许微生物组分类单元之间的灵活协方差，能够实现可扩展的计算，并有效地捕捉微生物组组成数据的其他关键特征，如零值的丰度。LTN结合了基于树的分解方法来有效地聚合稀疏的分类单元计数，并使用（多变量）逻辑正态分布 jointly 对树的分支处的相对丰度进行建模。潜在的高斯结构允许广泛的多变量分析和建模工具应用于高维数据——例如那些对协方差结构施加稀疏性或低秩假设的工具——能够被轻松地集成到模型中。作为一个通用的、完全生成的模型，LTN可以在各种场景中应用，同时，通过pólya-gamma增强实现共轭块吉布斯抽样（conjugate blocked Gibbs sampling），也提供了高效的贝叶斯推断方法。我们通过数值实验和对DIABIMMUNE研究中的婴儿队列的重新分析，展示了LTN在差异丰度分析的组成混合效应模型中的用途。我们通过数值实验和案例研究解释并展示了LTN如何通过适当考虑跨样本异质性，能够在不引入显式的零膨胀成分的情况下生成适当的零值比例。这证实了最近的一种观点，即基于计数的测序数据中的“零膨胀”通常是由于未考虑到的跨样本变异造成的。

**1 引言**

人类微生物组是居住在人体内或人体表面的所有微生物的遗传信息的集合。高通量测序技术的发展使得能够以成本高效的方式对微生物组分类单元的组成进行谱型分析，无论是通过鸟枪法宏基因组测序还是针对目标基因（例如16S rRNA基因）的扩增子测序。由此产生的数据通常是每个样本的组成计数向量，每个样本中识别出的微生物的总计数由研究的测序深度决定，只能按比例解释[1]。通过对微生物组组成数据进行监督和非监督分析的统一方法是通过构建生成模型来实现的，这允许直接整合这些数据的关键特征，包括计数性质、组成性、零值的丰度以及复杂的跨样本异质性。对于独立同分布（i.i.d.）或可交换的微生物组组成样本，生成模型可以作为构建更复杂统计模型的基础，以适应数据的额外特征。这种方法已被许多人有效用于差异丰度分析[2-7]和探索性分析[8-15]。这些先前的工作大多特别关注了适应微生物组数据中的组成性和零值普遍性的需求，提出了各种不同的方法，例如直接对零膨胀进行建模[7, 16, 17]，以及使用基于树的聚合来适应稀疏性[18-20]。然而，除了少数例外[21, 22]，这些方法通常对潜在的跨样本变异采用限制性的假设，并且往往在应用时没有明确警告使用者这些限制。一些最近的研究发现，下游分析对这些隐含的假设可能非常敏感，无论是在差异丰度分析[22]还是探索性分析[12, 23, 24]中都是如此。两种流行的微生物组组成生成模型是逻辑比率正态（LN）模型[25]和狄利克雷多项式（DM）模型[4, 10]。LN模型可以捕捉OTUs之间的丰富跨样本协方差结构；然而，当LN用于建模产生观察到的计数向量的未观察到的相对丰度时，如果OTU的数量哪怕只是中等数量，由于缺乏与通常用于捕捉微生物组组成数据计数性质的多项式采样模型的共轭性，推断在计算上将具有挑战性。因此，在实践中，LN最常用于建模实际观察到的计数比例（即OTUs上的计数除以每个样本的总计数），基本上忽略了测序读数的计数性质。另一方面，DM模型在计数上保持了多项式采样模型，但在未观察到的相对丰度向量上采用了狄利克雷分布。狄利克雷与多项式采样模型是共轭的，因此在计算上效率较高；然而，狄利克雷下引入的协方差结构由单一的标量参数（即狄利克雷伪计数参数的总和） caratterized，因此对于表征典型微生物组组成中的跨样本变异性来说过于限制性。DM模型的一种较新的泛化称为狄利克雷树多项式（DTM）模型[26, 27]，它利用了一个底层树结构（通常是系统发育树或分类树）来关联分类单元，这在一定程度上减轻了这一限制，但仍然非常有限制性。它的跨样本协方差虽然比狄利克雷的更宽松，但由于只使用参数来表征协方差结构，其中是分类单元的总数，因此仍然非常有限。此外，DTM施加的协方差结构的适当性在很大程度上取决于所施加的树结构（如系统发育树）在给定背景下准确反映分类单元之间功能关系的程度，这使得推断结果特别容易受到树规定的影响，这是采用基于树的方法进行微生物组组成分析时的一个主要问题[28, 29]。为了解决这些挑战，我们提出了一个新的生成模型，称为“逻辑树正态”（LTN），它结合了DTM和LN模型的关键特征，以继承它们各自的期望属性。LTN利用了DTM中的基于树的多项式采样分解方法，该方法利用了树聚合稀疏计数的优势[18, 19, 27, 30]，但与DTM不同的是，它在树的分支概率上采用了潜在的多变量LN分布，具有通用的协方差结构。这种更灵活的协方差结构使得推断在树的选择上比DTM更加稳健。同时，通过利用基于树的二项式分解多项式模型，并借助pólya-gamma（PG）数据增强[31]，LTN恢复了与吉布斯抽样的完全共轭性，从而避免了基于LN的模型所遇到的计算困难。LTN的完全概率、生成性质使其既可以作为独立模型使用，也可以嵌入到更复杂的模型中，以实现监督和非监督的目标。潜在的多变量高斯公式允许将常见的建模和计算技术（如基于协方差结构的稀疏性和低秩假设的各种类型正则化的技术）轻松地集成到LTN中，同时保持了典型微生物组分析数据大小（数千个分类单元和数百到数千个样本）的计算便利性。本文的其余部分组织如下。第2节在简要回顾LN和DTM后介绍了LTN，检查了其分布属性，包括其生成稀疏数据的能力，并提出了一个适用于差异丰度分析的基于LTN的混合效应模型。第3节我们通过几个数值实验研究了所提出方法在差异丰度分析中的性能，并将其与两种最先进的方法进行了比较。第4节我们对DIABIMMUNE研究中的T1D队列进行了案例研究[32]。第5节我们进行了一些讨论。

**2 方法**

**2.1 微生物组组成计数数据的结构**

已经开发了各种生物信息学预处理流程，如MetaPhLAN [33]和DADA2 [34]，用于“计算”每个样本中的微生物，并以扩增子序列变异（ASV）或操作分类单元（OTU）丰度的形式报告结果。OTU和ASV都可以作为微生物组组成下游分析的单位，近年来，由于多种原因，ASV已成为标准方法。然而，遵循微生物学研究的更广泛传统，在本文的其余部分，我们将使用更经典的术语“OTU”来指代它们中的任何一个。典型的微生物组数据集包括一个OTU（计数）表，其中是样本数量，是经过一些质量控制筛选后识别出的OTU数量，一个分类表，总结了在多个分类水平上将识别出的OTU与已知或未知分类单元的分类识别，以及可能通过医疗监测、问卷调查或其他方式测量的样本的一组协变量。人们通常还可以以系统发育树或分类树的形式构建研究中OTU的进化关系概要。有关更多细节和示例，请参见[1]。在建模和分析此类数据时需要考虑的微生物组组成数据的一些常见特征包括：（i）高维性——很大；（ii）稀疏性——OTU表中的许多计数是零值；（iii）普遍的跨样本异质性——OTU计数向量通常在样本之间表现出很大的变异性，并且OTUs之间存在复杂的协方差。我们注意到，微生物组数据中的稀疏性可以有两种形式：（1）样本内的稀疏性——每个微生物组计数样本中只有一小部分识别出的OTU具有非零计数；（2）OTU内的稀疏性——大多数OTU只在少数样本中被观察到。这两种类型的稀疏性通常都存在于微生物组数据集中。虽然我们提出的方法的主要动机在于高维性和复杂的跨样本协方差，但我们在后续的数值实验和案例研究中将提供示例，展示我们的模型如何有效地考虑这两种类型的稀疏性。

**2.2 两种流行的微生物组组成数据生成模型**

我们首先简要回顾两种最广泛使用的OTU计数生成模型，即LN模型和DM模型。它们启发了我们稍后介绍的新模型。我们首先介绍一些将在整篇论文中使用的基本符号。假设有个样本和总共个OTU，用表示。OTU表是一个矩阵，其第行表示样本中的OTU计数，第个元素表示该样本中OTU的计数。设是一个有根的全二叉系统发育树，覆盖这个OTU，其中表示内部节点、叶子和边的集合。形式上，我们可以通过其后代OTU的集合来表示系统发育树中的每个节点。具体来说，对于只有一个OTU的叶节点，我们让。然后可以从叶子到根迭代定义内部节点。也就是说，对于有两个子节点和的，我们简单地有。图1显示了一个包含四个OTU的节点的示例。图1展示了四个OTU的系统发育树示例。（a）节点的表示法。（b）从到的基于树的转换。假设是样本的OTU计数，是OTU计数的总数。给定总计数的OTU计数的自然采样模型是多项式模型

（1）

其中是底层的OTU（相对）丰度向量，位于一个-单纯形中。也就是说。相对丰度的单纯形约束可能会带来建模上的不便，传统的策略是应用所谓的逻辑比率转换将相对丰度向量映射到欧几里得空间[25]。三种最流行的逻辑比率转换选择是加性逻辑比率（alr）、中心逻辑比率（clr）[25]和等距逻辑比率（ilr）[35, 36]，这些都已经应用于微生物组组成数据。对于一个组成，clr和alr转换由给出

（1）

其中和是的几何平均数。另一方面，ilr转换使用二叉分治树结构来定义逻辑比率，它将转换为与树的内部节点相关联的“平衡”。与树的内部节点相关的平衡定义为

其中和分别是节点左右子树中的OTU数量，和分别是节点左右子树中的叶子概率组成的几何平均数。虽然在许多应用中找到合适的二叉树可能并不容易，但在微生物组背景下，系统发育树或分类树是一个常见的选择[37]。LN模型假设从这些变换计算出的对数比率是多变量高斯的。这些模型已成功用于表征微生物动态[38, 39]以及将协变量与微生物组组成联系起来[21, 40]。当OTU（ Operational Taxa Units）的数量增加时，使用LN模型进行推断可能会遇到计算上的挑战，这是因为多项式似然与多变量正态分布之间的共轭性不足。最近已经有一些努力通过近似计算来克服这一计算挑战。有兴趣的读者可以在[41]中找到关于LN模型所带来的计算挑战以及应对这些挑战的详细讨论。狄利克雷模型（Dirichlet models）是另一种经典的比例向量模型，它与多项式抽样模型共轭。因此，在微生物组建模文献中，它们通常被称为DM模型[4, 10]。DM模型在计算上比LN模型更高效，并且也能够生成非常稀疏的计数向量。然而，这类模型只提供了一个标量参数，即狄利克雷伪计数参数的总和，来表征OTU之间的协方差结构，这基本上假设所有OTU在组成约束所造成的人为依赖关系之外是相互独立的。对于典型的微生物组数据集来说，这显然是过于限制性的，因为这些数据集中的分类群展现出复杂的协方差结构[11, 42-44]。为了缓解狄利克雷模型的这一限制同时保持其计算可行性，引入了DTM（Density-Based Tree Modeling）模型[26]。例如，给定一个OTU之间的二叉树（如系统发育树或分类树），DTM利用了多项式抽样等同于沿着二元划分进行连续二项式抽样的事实。具体来说，方程(1)中的多项式抽样模型（在似然函数相同的意义上，只是有一个归一化常数不同）等同于在树的内部节点上定义的独立二项式抽样模型。

基于树的聚合方法已经在多种不同的微生物组数据分析方法中被采用，通常是因为系统发育树或分类树的容易获得性。一些例子包括上面提到的基于ilr（Invertible Learning Ratio）的LN方法，以及其他一些非生成式的树聚合方法[18, 19, 30, 45-47]。DTM的一个独特特点是，它不仅像ilr那样使用了基于树的转换来处理组成丰度向量，还将抽样模型（即多项式似然）分解为二项式模型，这是通过beta-二项式共轭性实现其计算可行性的关键。然而，DTM在解决跨样本变异方面的限制并不令人满意。需要注意的是，DTM中只有参数‘s’来表征OTU之间的协方差。实际上，DTM仍然假设这些参数‘s’是相互独立的，因此总共有 only 个标量参数，这些参数连同底层的系统发育树一起表征OTU之间的协方差。这也使得使用DTM进行推断时对系统发育树特别敏感。相比之下，LN模型，包括基于树的ilr模型，允许灵活指定OTU之间的协方差结构。这使得ilr对于二叉树的错误指定不那么敏感。此外，尽管DTM是完全生成式的，但实际上将其嵌入到更复杂的层次模型中并不容易。这是因为对于beta分布没有已知的（超）先验，因此如果存在额外的建模组件，例如对其实值的协变量效应和对未知浓度参数的先验，DTM将会遇到计算挑战。

2.3 LTN模型

我们引入了一个生成模型，它结合了DTM下的多项式抽样模型的基于树的分解和LN下的对数比率变换，从而同时享有两种模型的关键优势，并解决了它们上述的所有限制。特别是，我们的模型在表征分类群协方差和计算可行性方面都表现出灵活性，同时保持了适应稀疏性的能力。具体来说，我们对每个内部节点上的二项式分支概率应用了对数比率变换（在这种情况下只是一个逻辑变换）[49]。也就是说，我们将内部节点上的对数概率共同建模为多变量高斯分布。更正式地，对于一个组成概率向量，我们定义了基于树的对数比率（tlr）变换，其中。像其他对数比率变换一样，tlr将一个组成向量映射到一个对数概率向量。最后，将tlr建模为MVN（Multivariate Gaussian），我们就得到了LTN（Log-Ratio-based Tree Model）的完整公式：

(2)
我们将这个模型称为LTN。图2a展示了具有可交换样本的数据集的LTN模型的图形表示。（a）原始模型，（b）添加了Pólya-Gamma增强：引入了Pólya-Gamma增强以恢复共轭性。乍一看，虽然LTN提供了灵活的均值和协方差结构，但它像基于对数比率变换的现有LN模型（如alr、clr和ilr）一样缺乏共轭性。幸运的是，LTN下的似然二项式分解允许一种称为PG增强（PG augmentation）的数据增强技术来恢复共轭性。实际上，这种数据增强策略确保了即使有额外的建模层次结构（如协变量和多个方差分量），LTN模型也能保持完全共轭性，前提是采用任何共轭的（超）先验在高斯参数上，使得LTN比DTM更加通用且在计算上更易处理。图2b展示了具有可交换样本的数据集的LTN加上PG增强的图形表示。LTN的PG增强详细信息见支持信息部分A.1。

2.4 LTN的分布特性和零值的生成

在考虑将LTN应用于推理任务之前，我们首先考察了LTN的两个重要的分布特性，这些特性与它适应微生物组组成数据的常见特征的能力相关。特别是，尽管LTN的设计旨在丰富分类群的协方差并保持对维度的可扩展性，但仍然存在两个重要问题：(i) LTN在生成稀疏数据方面的能力如何——表现出OTU内部和样本内部的稀疏性；(ii) 结果分析对于二叉树的选择有多敏感？我们在这里提供了一些推理，并将在后续的数值实验和案例研究中通过实证证据来验证我们的论点。表征零值的丰度：微生物组组成数据的一个关键特征是在深层分类水平上通常观察到大量零计数。文献中已经指出，标准的生成模型不足以适应这种普遍的稀疏性，许多最新的发展都集中在引入一个零膨胀组件来适应这一数据方面[7, 16, 17]。我们承认在数据中适应零值的重要性，但与其他一些研究者[50, 51]持有相同的观点，即在生成模型中引入人为的零膨胀组件，除非有具体的生物学机制或实验协议来证明其必要性，否则可能会扭曲“非零”组分的实际含义的可识别性和解释性，从而导致下游推断的偏差。因此，除非在某些特殊情况下，我们不会在基于LTN的模型中引入显式的零膨胀组件，尽管从技术上讲构建它们并不困难。接下来我们认为，LTN的设计允许在没有显式零膨胀组件的情况下生成高度稀疏的计数，对于OTU内部和样本内部的稀疏性，原因各不相同。首先，我们注意到OTU内部的稀疏性——即，一个OTU的许多样本计数为零——是因为该OTU在大多数样本中的实际相对丰度较低。这种情况可以在LTN下轻松捕捉到，当包含所考虑OTU的树的子分支的分支概率较小时。具体来说，当OTU感兴趣的任何祖先节点上的均值参数（|μ|）较大时，就会出现OTU内部的稀疏性。另一方面，样本内部的稀疏性——即，一个计数向量有许多零，并且具有零计数的OTU在样本间变化很大——则是由于内部节点上的对数概率（log-odds）的边际方差非常大导致的。重要的是要注意，LTN因此可以在不改变分类群之间潜在的相关性的情况下适应OTU内部和样本内部的稀疏性，这对于确保结果分类群相关性的可识别性至关重要。虽然我们将在模拟和案例研究部分提供更广泛的数值证据，但在这里我们提供了一个简单的数值示例，来展示改变LTN内部节点上的对数概率的均值和方差对其影响，以及它们对OTU内部和样本内部稀疏性的影响，以支持上述论点。我们从DIABIUMME数据（稍后将描述）[32]中估计的参数和出发，生成了一个简单的LTN模型样本，假设所有样本都是可交换的。然后我们通过应用不同的缩放来修改和的值，从而调整和的绝对值以及的边际方差，分别导致图3中展示的不同水平的OTU内部和样本内部的稀疏性。这个模拟是简化的，仅用于演示LTN生成OTU内部和样本内部稀疏性的通用能力。我们将在为DIABIMMUNE项目的实际研究设计构建更复杂的基于LTN的混合效应模型时，检验LTN模型对DIABIMMUNE研究中观察到的零模式的拟合优度。（见图16。）

图3展示了按OTU（顶部）和样本（底部）划分的零值比例。顶部：样本是从具有固定均值和相关性 serta变化边际方差的LTN()生成的。底部：样本是从具有固定协方差和变化均值的LTN(生成的)生成的。基准参数值和是从DIABIMMUNE数据中的个体估计的和偏相关矩阵。这种现象是一般性的，并不特定于某个个体。对树的选择的鲁棒性：系统发育树和分类树通常可用于微生物组数据，因为它们是根据基因测序读数构建的，但在建模或分析微生物组组成时使用它们会引发一些自然的问题。毕竟，算法重建的进化树可能不一定能很好地反映给定背景下OTU的功能关系。这可能以两种方式发生。首先，即使进化关系重建得非常好，OTU之间的进化关系也可能与它们的功能关系不完全对应。例如，由于选择压力，进化上相距较远的OTU实际上可能在功能上相似或密切相关。理想情况下，应构建一个针对特定表型的功能树用于建模和分析，但这在大多数应用中通常是困难的。其次，即使OTU之间的进化关系确实与其功能关系接近，基于测序读数重建的底层进化树也常常存在相当大的不确定性，不同的生物信息学算法构建的树也可能在非平凡的方式上有所不同。因此，对于所有基于树的方法来分析微生物组数据来说，一个重要的问题是——如果所选的树对给定数据集中分类单元的功能关系描述不佳该怎么办？在构建和应用这些方法时，应该检查对树选择的敏感性。首先，我们提供一些关于LTN相比其他基于树的方法（如DTM）在树选择方面改进的鲁棒性的概念性推理，然后我们将提供数值证据来支持这一推理。树结构的选择可能会通过增加的不确定性和/或偏差来影响基于树的模型的后续推断。首先，请注意，树的选择会导致多项式均值的不同参数化。有效地将多项式均值重新参数化为对数几率（log-odds）应该涉及具有较大“结构变异性”的少数参数。（这里的“结构变异性”对应于感兴趣的变化，例如差异丰度分析中的组间差异和聚类中的簇间差异。）可以推测，在实践中，LTN表现出与其他基于树的变换（如Haar小波）观察到的“能量集中”现象类似的特点。为了理解这一点，考虑以下理想化的场景：两个分类单元（可以是任何分类级别）在功能上是相同的，因此它们会一起变化并且具有完美的相关性。在这种情况下，将它们放在一个共享父节点的兄弟节点上的树在所有样本中都会具有相同的分裂概率。特别是，如果有人想对两种条件进行差异丰度分析，那么在那个父节点上就不会有差异分裂。同时，由于基于树的聚合将两个相关分类单元的计数合并到一个分支中，两个分类单元在祖父节点上的信噪比将会很大，这与Haar小波变换的方式非常相似。这个例子说明了当分类单元在功能上相似时，树结构如何能够在更高级别的节点上集中统计信号，而不是将其分散到许多单独的分类单元上。最能有效地描述功能相似性的树将会将这样的功能上相同的分类单元完美地分成兄弟节点，而不会首先将它们分解成更小的子分类单元，从而导致“能量稀释”。实际上，所采用的系统发育树或分类树永远不会是这样一个完美的功能树，也不存在这样完全相关的分类单元。另一方面，错误的树选择将会产生相反的效果，通过在树的许多不同分割上分散功能上相似的OTU（Operational Taxa Units）的计数，从而导致普遍的“能量稀释”。在贝叶斯推断的背景下，这通常表现为后验不确定性的增加；而从抽样的角度来看，则表现为决策规则方差的增加。现在让我们转向第二个问题，即偏差膨胀的问题，我们认为LTN通常对此具有鲁棒性，而一些现有的基于树的方法（如DTM）则会受到显著的影响。错误指定的树可能会引入大量的偏差到最终的推断中，如果它对分类单元分布的潜在均值和协方差施加了强烈且不合理的约束。在DTM的情况下尤其如此，因为它假设树的分裂概率本质上是相互独立的，除了组成性约束之外。换句话说，这种相互独立性只有在系统发育树本质上是完美的树，并且任何两个分类单元之间的相关性都能通过它们在树中的距离准确地描述时才会发生，而这在几乎所有情况下都是生物学上不可信的。相比之下，LTN并不基于树来约束潜在的协方差结构（尽管可以通过对协方差的先验来纳入这样的约束）。因此，即使树的结构被错误指定，LTN也支持所有可能的协方差。在LTN中，错误的树选择通常仍然会导致效率的损失（后验不确定性和标准误差的增加），但通常不会导致较大的偏差（当然，除非LTN的先验基于树施加了强烈的约束）。我们进行了一个相对简单的模拟实验来证明这些点，并在支持信息部分C中报告了细节。具体来说，我们根据给定的树生成LTN数据，但使用严重错误指定的树来估计OTU的诱导均值和协方差。然后我们检查在正确和错误指定的树下对对数几率的均值和协方差的推断结果。我们的敏感性分析表明，即使树被严重错误指定，基于LTN的均值和协方差结构的推断通常是鲁棒的。我们将在第3.2节报告的模拟研究中提供更多关于LTN对树错误指定鲁棒性的证据，在这些研究中，我们生成了分类单元和感兴趣的结果之间的合成关联，而不考虑LTN采用的潜在树结构，然而LTN在识别差异丰度方面表现出了竞争力的性能。

3 在差异丰度分析中的应用

3.1 微生物组组成的混合效应回归模型

到目前为止，我们主要关注了LTN作为一个通用生成模型在微生物组组成计数方面的采样特性。接下来，我们将讨论其在实践中的应用。虽然原则上LTN可以作为单个组成丰度向量的先验，但我们认为在大多数应用中，例如典型的微生物组研究中，人们会使用LTN作为多个丰度向量样本的采样模型，这将是我们接下来讨论的重点。除了计算上的考虑之外，使用LTN进行推断原则上可以从贝叶斯或非贝叶斯策略（例如通过最大似然估计）来进行参数和学习。在这里，我们将关注贝叶斯方法，通过构建层次模型，将LTN作为（条件）可交换的微生物组组成样本的采样模块。构建此类层次模型的方法有很多。我们通过构建一个混合效应模型来说明这一点，鉴于此类模型在微生物组研究中的广泛应用，并重点关注检测两种条件下微生物分类单元的差异丰度这一常见任务（例如，病例与对照组等）。混合效应模型[21, 22]可以方便地纳入常见的设计特征，如协变量、批次效应、多个时间点和不同级别的重复实验。我们考虑一个常见的研究设计，其中数据集涉及可以划分为两个对比组的样本（例如，病例与对照组、处理组与安慰剂组等），并且研究人员感兴趣的是测试组间微生物组组成的差异。对于每个微生物组组成样本，让表示组指标，表示协变量的向量。我们考虑涉及测量分组结构的随机效应（例如，个体）。具体来说，我们假设每个样本都有一个与观察到的分组相关的随机效应。例如，在纵向研究中，可以表示收集第个样本的个体。我们通过以下线性模型将协变量和随机效应与微生物组组成联系起来：

（3）

其中表示潜在的组间差异，是协变量的未知固定效应系数的矩阵，是来自个体的随机效应，我们将其可交换地建模为。关于剩余的“噪声”，它表征了超出测量协变量的微生物组组成的额外变异。一个简单的选择是假设树节点上的独立噪声，并且可以采取逆伽马（inverse-Gamma）先验。有人可能会质疑独立噪声的假设，因为微生物组中的额外扰动通常也涉及OTU（Operational Taxa Units）的相关性。在采用树节点上的独立噪声时，我们做出了简化假设，即所有剩余的OTU协方差都包含在系统发育树结构中。可以通过在中的先验来进一步放宽这一假设，以纳入更灵活的OTU协方差。一种策略是采用潜在因子模型，该模型假设对噪声协方差进行低秩假设。实际上，这也可以通过修改模型（3）来轻松实现

（4）

其中表示样本的潜在因子，是相应的载荷。还可以通过使用时空过程来建模潜在因子来适应观测中的时空模式。为了阐述的简洁性，我们不再进一步讨论这个潜在因子模型的细节。模型的另一个需要指定的部分是对协方差矩阵的先验。当OTU的数量很大时，必须强制执行适当的正则化约束，以确保模型仍然是可识别的并且可以可靠地推断。现有的方法通常会对协方差或其逆矩阵施加低秩假设或稀疏性约束。对精度矩阵的常见先验是强制其向稀疏结构收缩。在高斯图上的此类模型的著名例子包括图形Lasso（gLasso）[52]和图形马蹄形（gHS）[53]先验。这两种模型都可以很容易地在LTN下使用。为了说明，在我们接下来的数值示例中，我们采用了gLasso先验，其形式为

（5）

其中DE是双指数密度函数，它本质上施加了正则化，EXP是指数密度函数，是正定矩阵的空间，是不起作用的归一化常数。在这个模型上进行后验采样可以使用[52]中提出的分块Gibbs采样方案。有了PG增强，模型的所有部分都可以从共轭完全条件分布中抽取。采样器的细节在支持信息部分A.2中提供。我们注意到，像gLasso和gHS这样的稀疏或收缩先验只是许多选择中的一种。另一个自然的替代方法是再次采用潜在因子模型，这些模型对协方差矩阵施加低秩约束，也可以很容易地在这里纳入，而不会复杂化推断算法。我们选择gLasso先验有两个原因。首先，虽然已经使用非生成方法（如SPIEC-EASI [54]）将稀疏逆协方差估计应用于微生物组网络推断，但在完全生成的微生物组模型中引入贝叶斯稀疏性先验的研究较少。LTN的生成性质允许将gLasso先验自然地集成到层次贝叶斯建模框架中，使得可以与协变量、随机效应和其他建模组件一起进行联合推断。其次，更重要的是，鉴于基于树的分解与Haar小波变换的类比，人们可能自然地期望分解具有类似的“美白效应”，就像在小波变换中对相关高斯观测值所观察到的那样，通过这种方式“波尔塔系数”的（条件）依赖性被减弱，从而导致精确矩阵的稀疏化。差异丰度分析通常涉及对两组差异的假设检验，并识别差异丰度的分类单元。在上述模型下，这可以通过在内部节点上测试一系列假设来实现：

为了进行这些局部测试，我们采用了一种贝叶斯变量选择策略，并在的条目上放置以下尖峰和平板先验：

其中是在零处的点质量，是预先指定的超参数。特别是，控制了节点上备择假设的先验边际概率，因此在控制多重检验方面起着关键作用。通常有两种选择的方法，要么直接设置，这对应于备择假设为真的节点的预期比例，要么设置，以便在期望的水平上诱导全局零（和备择）的先验概率。第一种方法不太严格，更适合于在控制误发现率（FDR）的同时识别差异丰度的分类单元。第二种方法更适合用于检验全局零假设（即两组之间没有差异的假设）。我们对协变量的系数采用连续先验分布。现在我们已经完全定义了模型，接下来我们转向决策理论的方法来识别差异丰度。在Gibbs采样从后验分布中生成样本后，我们可以计算Pr的MCMC估计值，这被称为每个变量的后验边际替代概率（PMAP），以及Pr的联合替代概率（PJAP）。当PJAP大于某个阈值时，可以使用它作为决策规则，以拒绝“两组之间任何微生物分类群都没有差异”的全局零假设。另一方面，PMAP可以用来确定哪些微生物分类群在不同组之间存在显著差异。PMAP和PJAP都可以直接从后验样本中估计得到。

3.2 数值实验

3.2.1 差异丰度分析
在进行案例研究之前，我们进行了模拟实验，这些实验的底层信号是已知的，目的是评估基于LTN的混合效应模型的表现，并将其性能与其他流行的差异丰度分析方法进行比较。为了避免指定微生物组组成的实际生成机制，我们通过向DIABIMMUNE项目[32]中的16S rRNA测序数据中注入信号来生成合成数据，该项目收集了来自芬兰和爱沙尼亚的婴儿的纵向微生物组样本。完整的数据集包含了33个受试者的777个样本中2240个OTU的计数。其中有25个受试者的样本数量超过20个，7个受试者的样本数量在某个范围内，只有1个受试者的样本数量少于10个。为了说明问题，我们关注相对丰度最高的100个OTU。这100个OTU的系统发育树被用于我们的基于LTN的混合效应模型中。我们改变了以下三个参数来控制模拟设置：样本大小、每个OTU的效应大小以及差异丰度的OTU数量。具体来说，对于每一种组合，在单次模拟运行中，我们按如下方式生成合成数据集：对于DIABIMMUNE队列中的每个受试者，我们随机抽取（不放回）该受试者的样本数量。如果数据集中该受试者的样本数量少于指定的数量，我们就抽取所有可用的样本。通过将样本随机分成两个相等大小的组（即组0和组1）来创建一个零假设数据集。另一个替代数据集是通过从相对丰度最高的20个OTU中随机选择OTU（不放回）并将这些OTU在组1中的计数增加某个值来生成的。对于每种组合，我们分别生成500个零假设数据集和替代假设数据集。我们考虑了单OTU差异和多OTU差异的情况。

值得强调的是，在我们的所有实验中，相关的单个或多个OTU都是随机选择的，不考虑底层的系统发育树结构，因此系统发育树与我们感兴趣的底层结果无关，相当于底层二叉树被严重错误指定的情况。因此，数值结果提供了LTN性能的保守估计，并证明了我们的基于LTN的模型对树结构错误指定的鲁棒性。我们将我们的基于LTN的混合效应模型在检验组间差异存在性方面的性能与MaAsLin2[55]和DirFactor[22]的性能进行了比较。MaAsLin2一次测试一个OTU的差异丰度，而与我们的方法类似，DirFactor基于OTU组成的联合模型，但不允许识别出差异丰度的单个OTU。我们在相同的模拟数据集上应用了这三种方法来检验两组之间的微生物组组成是否存在差异，同时调整了国家、采集时的年龄和样本的病例/对照状态，并考虑了受试者层面的随机效应。我们模型的性能、DirFactor和MaAsLin2是基于ROC曲线（图4）来评估的，用于检验两组之间微生物组组成的差异存在性。如3.1节所定义的，我们模型的（全局）检验统计量是后验联合替代概率（PJAP）。对于MaAsLin2，我们使用组合组标签（0/1）和OTU之间关联的最小Benjamini–Hochberg（BH）q值作为检验统计量；对于DirFactor，我们使用[22]中建议的组标签系数的范数作为检验统计量。

3.2.2 协方差估计
我们模拟了包含样本和OTU的数据集，其中OTU计数是根据DIABIMMUNE数据[32]中的系统发育树，使用LN和DTM模型生成的。模拟设置如下：

LN：OTU计数来自以下模型：

其中是基于系统发育树的ilr变换的逆变换[37]。的组成部分是从N中独立抽取的样本。我们采用了三种模型——Hub、Block和Sparse，这些模型被[42]认为是合理的微生物组组成协方差模型。Hub：100个OTU被随机划分为3个枢纽和97个非枢纽点，其中每个非零条目被设置为0.3（概率为0.5）或0.3（概率为0.5）。对角线条目的设置得足够大，以确保是正定的。Block：OTU被平均分成10个块。同一块中的每一对点以0.5的概率连接，而不同块之间的点以0.2的概率连接。每个非零条目被设置为0.3（概率为0.5）或0.3（概率为0.5）。对角线条目的设置得足够大，以确保是正定的。Sparse：设为，其中是一个对称矩阵，其下三角条目定义为，其中，和是独立的，并且足够大以确保是正定的。在这种模拟设置下，每个数据集中的系统发育树的内部节点都有一个稀疏网络。DTM：在这种情况下，样本来自于以下DTM模型：

DT是Dirichlet树分布。我们使用了DIABIMMUNE数据子集中的系统发育树，并根据相同的数据集设置了和的矩估计值。尽管这两种设置都在一定程度上利用了树结构，但我们注意到，在这两种设置下，LTN并没有假设相对丰度的正确分布形式。特别是，从ilr-LN生成的相对丰度在其他常用的对数比率变换下确实是高斯分布的，但从tlr变换下生成的则不是。我们故意使用与我们的模型不同的模型进行模拟，以展示即使模型有误，LTN仍然能够以鲁棒的方式描述真实的协方差结构。为了评估我们方法的性能，我们使用了一种名为COAT的先进方法作为基准，该方法最近在[42]中提出，并考虑了四种不同的损失函数——Frobenius范数、范数、逐条目范数以及估计的clr相关矩阵与真实值之间的谱范数。我们检查这些损失函数是因为COAT旨在估计clr协方差，因此在使用COAT作为基准时，这个选择是最合适的。值得注意的是，虽然我们提出的模型明确地在tlr变换下对协方差进行了建模，但我们的模型是完全生成的，因此可以基于推断的模型评估在clr变换下的诱导协方差——或者任何其他变换。这与COAT程序不同，后者不是基于生成模型的，因此只能用于估计clr协方差。这就是我们特别将我们的clr协方差估计与COAT进行比较的原因。我们注意到，在我们的方法和COAT中，关于稀疏性的假设是不同的——我们的方法假设在tlr变换下的精度矩阵是稀疏的，而COAT直接假设在clr变换下的协方差矩阵是稀疏的。在每个模拟轮次中，我们运行LTN的Gibbs采样器10,000次迭代，并丢弃前半部分作为燃烧期。结果在表2中呈现。在这里，我们在LTN中设置了。在LN-sparse设置下，LTN在所有四种损失函数下都优于COAT。在其他模拟设置下，LTN在大多数损失函数下也优于COAT。多变量正态的灵活性允许LTN即使在相对丰度分布被错误指定的情况下也能有效地描述样本间的变异性。

3.2.3 LTN的生成能力
为了展示LTN的生成建模能力，我们分析了American Gut Project（AGP）微生物组数据集[57]的一个子集，该数据集是一个通过QIIME管道处理的开放访问16S rRNA测序队列。根据[12]的分析，我们专注于2016年7月从被医学诊断为炎症性肠病（IBD）的参与者那里收集的粪便数据。在原始的OTU表（包含27,774个OTU）中，我们应用了[12]中描述的相同的OTU过滤方法：首先保留总丰度最高的75个OTU，然后选择在这些OTU中总计数至少为500的样本。这种预处理后得到最终的数据集，包含189个样本。我们使用Gibbs采样方法对这个数据集进行了两次精度矩阵规格下的LTN模型拟合——分别是稀疏型（图解套索先验）和对角型——每次迭代10,000次，并丢弃前5000次作为燃烧期。从后验样本中，我们计算出后验均值，这些均值将作为我们模拟研究的生成参数。具有对角协方差的LTN模型可以作为狄利克雷树多项式（DTM）的替代模型。为了评估LTN的生成性能并展示LTN与其他模型之间的差异，我们进行了以下对比：（1）狄利克雷多项式；（2）非生成模型COAT，其样本来自均值为且协方差等于COAT估计值的对数比率正态分布。对于每种模型，我们根据每个样本的总OTU计数条件生成合成数据集，并通过以下方法进行比较：

- 零比例分布：评估模型是否捕捉到了OTU计数中特征性的大量零值。
- NMDS图：在低维空间中比较生成样本和真实样本之间的关系。
- 堆叠条形图：提供跨分类单元相对丰度分布的视觉评估。

LTN及其他模型的生成性能在图6-8中展示。图6表明，LTN有效地捕捉到了观察到的每个分类单元的稀疏模式。实际上，生成数据的每个分类单元零比例的经验分布与观察数据之间的Kolmogorov–Smirnov距离分别为0.106、0.159、0.333和0.111，分别对应于使用图解套索先验的LTN、使用对角协方差的LTN、狄利克雷多项式以及使用COAT估计值的LN，其中使用稀疏精度矩阵（带图解套索先验的LTN生成的数据表现出与经验零比例分布最小的Kolmogorov–Smirnov距离。关于低维结构，如图7所示，LTN生成的样本在NMDS投影中占据了与观察数据相似的几何空间。最后，图8显示LTN生成的样本紧密再现了观察数据的相对丰度模式。值得注意的是，LTN准确捕捉到了拟杆菌目（最丰富的分类单元）的分布以及根瘤菌目和芽孢杆菌目的样本间的变异性。这一比较分析进一步表明，LTN和使用COAT估计值的LN模型比狄利克雷多项式方法更好地捕捉了样本间的变异性，这突显了在生成微生物组分析中灵活协方差建模的重要性。

图6展示了分类单元的零比例分布。子图(a)-(e)中的红线表示所有OTU计数中零计数的比例。(a) 观察值。(b) 使用稀疏精度的LTN。(c) 使用对角精度的LTN。(d) 狄利克雷多项式。(e) 使用COAT估计值的LN。

图7展示了使用COAT估计值的LN生成的样本（作为DTM的代理）、使用对角精度矩阵的LTN以及真实数据（“观察值”子图）的NMDS图。子图使用相同的NMDS轴。

图8展示了观察值和模拟的相对丰度。(a) 观察值。(b) 使用稀疏协方差的LTN。(c) 狄利克雷多项式。(d) 使用COAT估计值的LN。(e) 使用对角协方差的LTN。

除了生成性能之外，了解LTN学到的协方差结构提供了关于模型如何捕捉分类单元间关系的关键见解。LTN的后验均值及其对应的关联矩阵在图9中进行了可视化展示，揭示了节点之间的显著非零协方差和相关性，因此需要节点间的灵活协方差结构，而不是DTM假设的独立性。

图10显示了对角LTN、图解套索LTN和COAT的估计clar相关矩阵。LTN的不同变体表现出不同的聚类模式，反映了基于系统发育信息的协方差结构，其中对角LTN产生的相关矩阵比图解套索LTN更稀疏。图解套索先验提供的灵活协方差结构使得OTU之间的依赖网络更加丰富，捕捉到了更复杂的分类单元间关系，而更受限制的对角规范则无法做到这一点。相比之下，COAT估计没有表现出这样的聚类模式。尽管在这个真实数据集中真实的clar相关结构仍然未知，但我们之前的模拟研究（第3.2.2节）表明LTN的估计比COAT的估计更接近已知的真实相关性。这一验证表明，LTN和COAT估计之间的差异可能反映了真实的潜在相关性，而不是方法上的人为误差。

4案例研究：DIABIMMUNE数据

4.1 数据和分析概览

现在我们对DIABIMMUNE项目的1型糖尿病（T1D）队列进行了更深入的分析，以研究微生物组组成与多个变量之间的关系。这项研究的主要目的是比较那些发展成T1D或血清自身抗体的婴儿（预测T1D发作的标志物）与同一地区的健康对照组之间的微生物组差异。在研究的时间范围内，33名婴儿中有11名血清转化为自身抗体阳性，其中4名婴儿发展成了T1D。之前对该数据集的研究已经建立了微生物组组成与T1D状态以及饮食因素之间的关联。例如，研究表明T1D患者的微生物组往往具有较低的β-多样性以及拟杆菌目分类单元的较高占比[32, 58, 59]。在接下来的分析中，我们重点关注血清转化状态和饮食协变量，并研究微生物组组成与这些变量之间的关系。我们再次应用了基于LTN的混合效应模型，并尽可能将结果与MaAsLin2和DirFactor的结果进行比较。我们分析的主要发现总结如下：(i) 血清转化后采集的样本的微生物组组成在拟杆菌属、拟杆菌亚属和红球菌科的相对丰度方面与其他样本不同，这些分类单元与[32]中识别的T1D相关分类单元不同；(ii) 引入多种食物和停止母乳喂养可以显著改变婴儿的肠道微生物组。我们将在以下部分详细介绍我们的研究。

4.2 不同血清转化状态下的差异丰度分析

研究结束时个体的T1D状态被记录为三个等级之一：对照组、血清转化组和T1D发病组。这一状态用于定义原始研究中的病例组和对照组，其中病例组包括血清转化个体和临床诊断为T1D发病的个体。对于病例组中的个体，还记录了血清转化时的年龄。之前对T1D状态的分析显示，在“血清转化窗口”之后（即所有血清转化和T1D受试者血清转化年龄的第一和第三四分位数之间的时间），T1D受试者的微生物群落的α-多样性显著下降，而其他受试者的α-多样性继续增加，这种差异可以由几组细菌解释[32]。受这些发现的启发，并受到[22]中对这些数据的重新分析的鼓舞，我们对血清转化前后收集的样本进行了差异丰度分析，以识别血清转化后微生物组组成的变化。我们注意到，虽然[22]也研究了微生物组组成与血清转化状态的关联，但他们专注于估计个体分类单元的协变量效应，而没有提供正式的方法来测试分类单元的差异丰度。遵循[32]的方法，我们还包括了八个二元饮食变量，这些变量指示在样本收集时是否摄入了特定类型的饮食：母乳喂养、固体食物、鸡蛋、鱼、大豆制品、黑麦、大麦和荞麦及小米。饮食模式被认为显著影响微生物组组成，因此必须将其纳入考虑。我们还包括了收集时的年龄（对数变换）、性别和国籍作为协变量。探索性分析：图11使用Bray-Curtis差异性绘制了多维缩放图，其中收集时的年龄几乎沿着第一个轴的方向增加，芬兰和爱沙尼亚样本在第二个轴上大致分开。年龄和不同国家之间的总体趋势明显，以及同一年龄组内受试者之间的显著差异。图12展示了每个受试者随时间的变化最丰富的门类，图13展示了随时间的门类组成。它们在出生时以及不同个体之间都显示出较大的变异性，这种变异性随时间逐渐稳定。因此，我们使用对数年龄作为协变量来解释早期微生物组组成的快速变化。有趣的是，血清转化通常发生在高度变化的阶段之后。大多数一岁后收集的样本主要由拟杆菌目或厚壁菌目主导。此外，在爱沙尼亚个体收集的样本中，大约有一半的样本中放线菌目占主导地位，而在芬兰个体中则不是这样。

图11展示了样本之间Bray-Curtis差异性的多维缩放图。

图12展示了每个样本中最丰富的门类。

图13展示了随时间的门类组成。

基于模型的差异丰度分析：接下来，我们应用基于LTN的混合效应模型来识别血清转化状态下的OTU组成差异，并关注相对丰度最高的前100个OTU。具体来说，如果样本是在血清转化后收集的，则样本被标记为1；否则标记为0。我们包括了可能影响微生物组组成的环境和饮食因素作为固定效应。此外，还包括了个体随机效应来解释个体间的差异。我们使用相同的固定和随机效应集合对MaAsLin2和DirFactor进行了比较。为了测试全局零假设，我们设置了γ，对应于全局零假设成立的先验概率为0.5。相应的PJAP接近于1，表明有压倒性的强证据反对零假设，这与MaAslin2和DirFactor的结果一致，这两种方法在1000次置换下的置换p值都小于1000。对于识别差异丰度特定分类单元的更有趣任务，我们设置了θ，对应于有差异分裂概率的节点的先验预期比例为5%。这种更为宽松的多重检验调整策略更符合设置差异丰度分类单元预期数量的实践。DirFactor实际上不允许直接评估个体分类单元级别的关联测试，因此被排除在比较之外。

这是所有报告为显著的节点上的γ的平均值。在这里，我们采用了一定的决策规则来处理某些ζ。为了控制ζ，我们将PMAPs排序为?ζ，并在决策规则中设置了PMAP的阈值βζ。

图14展示了我们方法为系统发育树的所有内部节点估计的PMAPs。为了便于比较，在同一张图中，我们用MaAsLin2报告的经过BH调整的q值（和）来标记叶节点（即OTUs）。四个具有最高PMAP的节点表明，在血清转化期间，Bacteroides属、Erysipelotrichaceae科和Parabacteroides属内的OTUs的相对丰度发生了变化。特别是，我们的方法进一步表明，在血清转化后，与系统发育树上的姐妹物种相比，P. distasonis的丰度增加了，这与数据中观察到的P. distasonis的相对丰度一致（见图15a）。这些关联得到了基于其他数据集的先前研究的支持[60, 61]。虽然我们当前的分析主要关注血清转化状态的主效应，但图15a和b表明，对于某些分类单元，血清转化组和非血清转化组之间的年龄轨迹可能存在差异，这表明年龄和血清转化状态之间可能存在交互效应。LTN混合效应对模型可以通过扩展模型以包括协变量之间的交互项来自然地容纳这种交互效应，并对主效应和交互系数使用尖峰和平板先验。在这种情况下，后验推断可以类似于当前的Gibbs采样器进行，尽管存在交互项时解释会变得更加复杂。我们注意到，尽管MaAsLin2报告了与血清转化水平相关的大量显著OTUs，但鉴于我们在多OTU关联情况下的模拟研究中观察到的经验FDR膨胀（如图5所示），我们认为其中可能包含相当数量的假阳性。

图14 在图形查看器中打开

基于LTN的混合效应模型的节点的后验边际替代概率（PMAPs）和MaAsLin2报告的OTUs的BH q值，在血清转化前后样本的比较中。MaAsLin2报告的BH q值为 和 的OTUs分别用“**”和“*”标记。分支上的“+”和“”根据 的后验均值的符号进行标记，表示在案例中是左边的子节点还是右边的子节点的相对丰度增加了。

图15 在图形查看器中打开

P. distasonis和OTU 4439360的相对丰度。对于P. distasonis，省略了第一次血清转化事件之前的数据。(a) P. distasonis。(b) OTU 4439360。

4.3 对OTU内和样本内零值的拟合优度

为了检验我们提出的模型是否能够充分捕捉数据中的零模式，我们将带有所有协变量系数矩阵正态先验（即没有尖峰和平板成分）的LTN混合效应模型拟合到DIABIMMUNE数据上，并生成OTU计数的后验预测样本。图16显示了零比例的后验预测分布。总体而言，我们的模型能够在不显式建模零膨胀成分的情况下捕捉到该数据集中的OTU内和样本内的稀疏性，这与[51]和[50]中的论点一致，即在计数采样模型下适当建模交叉样本异质性通常足以描述测序数据中的大量零值。

4.4 基于饮食变量的差异丰度分析

作为展示我们方法性能的额外示例，并将其与MaAsLin2进行比较，我们对每个饮食因素的微生物组组成重复了相同的丰度分析。由于空间限制，我们将在支持信息部分D中报告所有结果。

5 结论

我们引入了一个基于树的生成模型，称为LTN，用于微生物组组成。LTN将多项式采样模型分解为系统发育树内部节点上的二项式集合，并将OTU组成转换为节点特定的对数几率，并在对数几率上采用多元正态模型，从而提供了OTUs之间的灵活协方差结构。通过基于树的分解，LTN避免了LN模型中多项式和多元正态之间缺乏共轭性所带来的计算挑战，并通过引入PG辅助变量使用Gibbs采样进行高效推断。这为采用各种多元分析模型/方法进行微生物组分析打开了大门，同时保持了计算的可行性。由于其完全生成的特性，LTN既可以用作独立模型，也可以作为更复杂模型（如涉及协变量效应和潜在结构的模型）的组成部分，这些在微生物组应用中通常是需要的。我们提供了用于组间比较的混合效应建模示例。我们展示了如何通过贝叶斯模型选择框架进行微生物组组成的组间比较。在这些示例中，我们采用了gLasso先验进行演示，但我们的模型/先验选择可以根据实际情况由实践者自行替换。最后，除了组成数据之外，基于树的对数几率分解还可以用于建模一般概率分布[49]，本文描述的建模和推断策略也可以应用于这些更广泛的情境中。

致谢

本研究部分得到了NIH拨款R01-GM135440以及NSF拨款DMS-1749789和DMS-2013930的支持。这项工作的大部分是在L.在杜克大学完成的。

资助

这项工作得到了美国国立普通医学科学研究所（拨款号R01-GM135440）和美国国家科学基金会（拨款号DMS-2013930和DMS-1749789）的支持。

利益冲突

作者声明没有利益冲突。

数据可用性声明

支持本研究发现的数据可以在DIABIMMUNE微生物组项目中公开获取，网址为https://diabimmune.broadinstitute.org/diabimmune/。本文中的分析代码可以在Github上找到，网址为https://github.com/MaStatLab/LTN_analysis1，而LTN模型的R包可以在https://github.com/MaStatLab/LTN找到。