微生物组数据系统发育分析的方法

Methods for phylogenetic analysis of microbiome data

Nature Microbiology, [14.1], Review, 2018-05-24

原文链接: https://www.nature.com/articles/s41564-018-0156-0

第一作者：Alex D. Washburne 1,6*, James T. Morton2,3,6

通讯作者：Alex D. Washburne 1,6*

其它作者：Jon Sanders 3, Daniel McDonald3, Qiyun Zhu3, Angela M. Oliverio4,5 and Rob Knight 2,3

主要单位：1. 蒙大拿州立大学，微生物与免疫系；2. 加州大学圣地亚哥分校，计算科学系；3. 加州大学圣地亚哥分校，儿科；4. 科罗拉多大学，生态与进化生物系；5. 科罗拉多大学，环境科学联合研究所；

热心肠导读

Nature子刊：菌群系统发育分析的方法

Adam刘永鑫-中科院-宏基因组 09-14 热心肠日报

原标题：微生物组数据系统发育分析的方法

① 四种研究方法依次是：比较方法用来分析进化上不同起源的微生物之间的多种表型差异；

② 祖传表型推断基于物种进化史和亲缘关系推断共祖群体的共有表型；

③ 系统发育参数通过构建具有生物学意义的变量简化或概括群落特征；

④ 还可结合亲缘关系描述样本间距，用于样本聚类、拟合等；

⑤ 此类分析基于物种间的系统发生关系，受所选基因、基因水平转移等因素影响很大；

⑥ 必须小心选择研究问题和研究假设，具体可参考作者在线提供的实战教程。

主编评语：复杂菌群中的个体并非是孤立存在的，而是具有一定亲缘关系的。因此，在分析菌群表型时理应将物种之间的系统发育信息考虑进来。Nature Microbiology上的这篇综述总结出4种不同的研究问题和分析方法，更在网站上在线提供了具体分析流程，是进行菌群数据系统发育分析时非常有价值的参考材料。

摘要

了解微生物的进化历史如何影响我们对微生物数据集的分析？根据研究问题，微生物的共同祖先可能是混杂变异的来源，或用于推断的支架。例如，当对性状进行回归时，共同祖先是观测值（观测变量）之间依赖性的来源，而当搜索具有相关丰度的进化枝时，共同祖先是推断的支架。微生物的共同祖先及其基因以树状结构呈现-所谓的系统发育树 - 而这个系统发育的信息可以并且应当将其纳入微生物数据集的分析中。虽然最近关于系统发育的分析工具愈来愈多，但对于哪种方法能够最好地回答哪些生物学问题却鲜有指导。在这里，我们回顾了微生物组数据集的系统发育知识的分析方法，同时也回顾了选择适当方法的注意事项以及这些方法中固有的挑战和局限性。我们介绍了这些工具概念，将其分类为系统发育比较方法、祖先状态重建、系统发育变量和距离分析 ，并在补充在线教程中提供示例。对研究问题的仔细考虑以及生态学和进化假设将有助于研究人员选择系统发育和适当的方法来产生准确的、更有生物学意义的和新的见解。

背景

高通量测序可以提供有关微生物群落的信息，其数量超出了我们理解它的能力。大多数微生物类群从未被培养或通过实验表征。对于多数实验而言，我们只有序列片段和一些远亲的全基因组序列，以及捕获微生物进化历史的生命之树。我们如何组织和分析有关未被描述的微生物及其序列片段的海量信息（The deluge of information）呢？

组织生物多样性的两个基本工具是分类和系统发育（taxonomy and phylogeny）。分类学基于从三个界域（细菌，古细菌和真核）到几百万的种的层级分类名称来对微生物进行分类。系统发育是对微生物进化历史的估计，并通过对应于估计的进化事件的一系列分裂节点对每个生物进行分类，其中最近的共同祖先（a most recent common ancestor）被指定形成两个子物种（daughter species）。

微生物分类学和系统发育最终可能是等同的，系统发育中的每个进化枝都具有分类学名称。然而，目前的分类学分类很粗糙；现代分类只是对系统发育中的一小部分的分枝进行了分类命名。目前，系统发育是用于微生物分类的更详细的支架（骨架）。

系统发育是组织和理解微生物世界的工具。由于相关生物往往倾向于具有相似的性状，即使我们无法直接度量它们，系统发育也可以将这些性状纳入我们的分析中。系统发育是对谱系进行分类和推断功能性生态性状的支架，即使对于未经分类学或生理学分类的谱系也是如此。通过高通量分类可以加速微生物生态学，并且可以通过系统发育进行推断。资源消耗，栖息地关联和物种相互作用是性状的原因和结果，使用系统发育推断或潜在的与性状一起作用可以增强我们对影响人类健康、生物地球化学和气候变化的微生物群落控制的能力。

系统发育如何帮助分析微生物组数据？不同的研究问题需要考虑如何使用系统发育来修正统计分析呢？例如，性状之间关联的统计检验应该将系统发育视为观测值之间依赖性的来源，而寻找更简单的合并物种方式的研究应该将系统发育视为可能的物种合并的骨架。有越来越多的文献报道了关于分析带有系统发育结构化数据的方法，这些方法在他们寻求回答的问题中存在细微却本质的差异。因而需要将多领域的方法简化为几类概念上不同的方法集，从而为用于分析系统发育结构化数据的方法提供指导、比较和开发框架。

在本综述文章中，我们通过讨论主要的几类方法来组织系统发育结构化数据分析领域。我们首先强调该领域的一个基本问题：估计的系统发育并不完善。然后我们定义了进行估计的四类方法：（1）比较法; （2）祖先状态重建（ancestral state reconstruction）和后代性状推测插补（trait imputation）; （3）系统发育变异; （4）系统发育感知距离（phylogeny-aware distances）（表1）;并提供补充在线教程的示例（https://knightlab-analyses.github.io/phylogenetic-tutorials/ ）。大多数统计工具可以进行系统发育感知的重分析，但我们涵盖的类别囊括了最常用的和积极开发维护的几大类。

我们讨论了微生物组数据的系统发育感知分析的挑战，包括水平基因转移（HGT）和选择在构建系统发育时使用哪些基因。通过将文献划分为几类理念上不同的方法，我们为微生物组数据分析中这些重要方法的开发和实施提供了一个通用框架。有关相关术语的词汇表，请参见方框1。

正文

表1 使用系统发育的方法和对应回答的科学问题

方法	描述	示例	应用
比较方法	找出性状之间的联系，在系统发育上控制进化	16S rRNA基因拷贝数与体内生长速率相关吗?	PGLS18、Paired t-test
祖先状态重建	在系统发育过程中对历史谱系的性状值进行估算，并利用祖先性状对当代物种的性状值进行估算	根据OTU亲属的16S rRNA拷贝数，对OTU的16S rRNA基因拷贝数的最佳估计?	PICRUSt2
系统发育变量	利用系统发育来构造生物学上可解释的变量(例如，一个分支的丰度)，并简化/总结群落中的特征算	生境组合P个不同独立类群的哪些内部边缘? Faith的系统发育多样性是如何随pH值变化的?	Diversity analyses、Taxonomic analyses、Phylofactorization、EdgePCA、PhILR
系统发育感知距离	利用系统发育来构造样本之间的距离，然后用这些距离来修改用于分类、正则化回归等的统计工具	两种微生物群落有何不同?	UniFrac、Inner product methods

方框1 |专业术语

Ancestral state 祖先态。祖先物种的性状，通常是祖先生物的表型和基因型的估计。

Ancestral state reconstruction 祖先态重建。在系统发育的不同点上推算插补祖先状态。（Imputing the ancestral states）

Bayesian inference 贝叶斯推断。给定祖先状态的先验置信集以及观测到的现有物种的表型/基因型，贝叶斯方法将试图获得对祖先状态的更准确的估计，以及预测的置信度。

Blomberg’s K.与Pagel’s λ（见下文）相比，更为常见的当今系统发育信号测度，从0到无穷大，表示随时间演变的加速或减速程度。

Bootstrapping 系统发育自展支持率。 重采样，统计学家Efron 提出，Felsenstein提出的，将其应用到系统发育领域，重复随机重建系统发育树，通常用于评估每个进化枝的重建自展支持百分比。

Brownian motion 布朗运动。连续随机游走，其中跳跃是正态分布的随机变量。通常在PCM中用作从祖先节点到树的叶子节点的连续性状进化的零模型，其中随机游走在系统发育中的那些节点分支。在布朗运动进化模型下，物种观察到的性状之间的协方差与其共有祖先的分支长度成正比。

Classification 分类。回归或其他预测类别因变量的方法，有监督的学习方法。

Clustering 聚类，创建识别变量组的分类器（分类变量），例如将具有比较高的组内相似性和较低的组间相似性的组进行判别归类，无监督的学习方法。

DNA amplicons DNA扩增子。人工扩增事件的DNA产物，例如16S rRNA基因的聚合酶链反应扩增的所得产物，其随后被测序并计数以组装微生物组数据集。有时可以使用扩增子构建微生物的准确系统发育。

Edge 边。系统发育中的结构，表示在某个时间点内假设的不同的，不间断的谱系。

Edge lengths 边长。边长可以表示历史沿持续的时间或将其祖先与子节点分开的突变事件的数量。

EdgePCA 一种对一组变量vi进行主成分分析的方法，对应于沿每条边的丰度差异，即i。

Epistasis 上位。当两个或更多个遗传基因座相互作用以确定表型性状时。

Evenness均匀度。alpha指标的一般术语，表明群落与所有物种具有相同丰度的接近程度。

ILR Isometric log-ratio等距对数比。对数变换数据的算术平均值的标准化差异。通常用于在对数尺度上更适合分析微生物数据集，但非对数转换数据的类似差异是双样本t检验的t统计量。

Maximum likelihood 最大似然。最大似然法将祖先状态视为未知参数。给定概率进化模型，最大似然方法将尝试优化这些参数以试图找到产生我们在当前已知物种中观察到的性状的最可能的祖先状态。

Maximum parsimony 最大简约。最大简约性试图通过最小化祖先和当前后代之间的性状变化的数量来重建祖先状态。

Monophyletic 单系。一组物种被称为“单系物种”，相对于更大的物种，如果它们最近的共同祖先除了物种组之外没有其他后代。

Node 节点。当一个谱系分裂成两个或更多个不同的谱系时，系统发育中的结构表示物种形成的假设时间。

Pagel λ。系数发生信号的测量值，范围在0和1之间，表示近亲之间性状相关性与特性进化的布朗运动模型相匹配的相对程度。

PhILR。系统发育等长对数比。需要完全解析的系统发育（即没有多重结构）的数据转换。 PhILR变换不是用每个物种的一个变量表示数据，而是表示系统发育中每个节点有一个变量的数据。使用ILR变换构造变量以对比从每个节点下降的姐妹分支。

Phylofactorization。一种通过广义图分区算法选择变量的方法。通过首先考虑沿边的对比来构建变量，例如对比鸟类和非鸟类的差异或ILR，然后找出哪个变量使研究者的目标函数最大化。系统发育沿着该边分开，并且重复该过程，仅限于发现边的亚系统发育的对比（例如，在分割鸟类/非鸟类之后，分隔鸽子/非鸽子的边代替鸽子/非鸽鸟）。

Phylogenetic comparative methods 系统发育比较方法。校正近亲之间性状观测的相关性的统计方法。 PCMs通常使用进化模型来计算随机进化下预期的近亲观察值之间的相关性 。

Phylogenetic distance系统发育距离。沿着连接系统发育中两个物种的路径的边长的总和。

Phylogenetic inference系统发育推断。估计一组基因的进化历史。

Phylogenetic variables系统发育变量。借助于系统发育构建的变量（包括所有物种来源于同一多歧节点（Polytomy）的星状系统发育star phylogeny）。与系统发育距离相反，变量指示具有生物学意义的变化方向和轨迹。

Phylogeny 系统发育。一组基因进化史的图解假设。系统发育树可以是有根树，暗示基因组的最原始的共同祖先，也可以是无根的。

Polytomy 多歧（支）节点。具有两个以上子系的节点。通常，polytomies代表了历史物种形成事件的精确时间的不确定性。

Regression 回归，估计变量之间关系的一个数学预测模型。

Shannon diversity 香农多样性。均匀度H的一种特殊度量，定义是通过计算一组总和为1的相对丰度pi，通过一个xlogx 函数计算而得：H = -sum（pi log（pi））。

系统发育推断

生命之树是未知的 - 它是被估计出来的，并且系统发育的准确性提高了系统发育结构化数据分析的准确性。通常通过收集基因序列，基于同源性比对序列，并使用突变模型来推断最可能的进化历史来估计微生物系统发育。而估计的系统发育可能会根据测序的基因，序列位点的比对方式，以及使用哪种进化模型和推断历史的方法的不同而大相径庭。系统发育推断中的错误可以蔓延到系统发育结构化数据分析中的错误。在这里，我们讨论系统发育推断和系统发育感知分析之间的相互影响;对于系统发育推断方法的综述，读者可以查阅该主题的重点综述文章Nei, M. & Kumar, S. Molecular Evolution and Phylogenetics，Oxford Univ.
Press, Oxford, 2000 和 Yang, Z. & Rannala, B. Molecular phylogenetics: principles and practice.
Nat. Rev. Genet. 13, 303 (2012)

人们可以构建任何基因的系统发育；不同的物种在同一基因的拷贝数，以及系统发育的分辨率和各种性状的系统发育信号方面各不相同。 16S rRNA基因通常用于细菌和古细菌的系统发育推断，但也可以用其他基因，如β-内酰胺酶及其亲属构建系统发育，产生一棵以出现抗生素抗性性状逐渐增加的趋势为边的系统发育树。微生物真核生物同样具有许多可用于系统发育推断的基因，18S rRNA基因是最常用的。

选择用于系统发育推断的基因决定了与系统发育相关的一组性状。细菌基因组树通常与16S rRNA基因（16S）衍生的系统发育相关，但16S树和基因含量之间的相关性随谱系和系统发育深度而变化 。基因的水平转移（HGT）破坏了16S树和基因含量之间的相关性，因为它允许具有远距离16S基因的细菌共享同一个重要的性状，例如致病岛（pathogenicity islands）和抗生素抗性基因。此外，16S序列具有多个可变区，并且可以在单个基因组内有多个不同的拷贝，使得系统发育推断举步维艰。更复杂的情况，例如当上位性作为功能性生态性状的基础并且其中一个上位基因可以水平传播时，就会阻止我们去明确使用哪个基因树。

用于分析系统发育结构化数据的不同方法使用系统发育的不同性状。在系统发育中通过将许多分支进行向上聚合的距离法和系统发育比较法（PCM）对系统发育推断中的错误更加稳健和鲁棒。依赖于少数分支的方法对系统发育推断中的错误更敏感 。对于依赖于少数内部节点或分支的方法，系统发育的推断的不确定性 - 特别是临界分支的单系（monophyly of critical branches）的自展支持（bootstrap support） - 对于并入下游数据分析可能是重要的。在单系是至关重要且允许多歧分支（polytomy）的情况下，研究人员可以将已解析的细粒度节点折叠（collapse）成多歧节点，以改善整个树的引导自助支持 。更确定但粗粒度的系统发育可能优于不太确定但完全解析的系统发育 。

图1|系统发育比较方法(PCM)控制由沿着系统发育树的性状进化产生的性状之间的统计依赖性

a，一个极端的系统发育与两个远缘相关的进化枝。如果将性状进化模拟为对系统发育的随机游走，则两个远缘相关的进化枝将驱动性状之间的协方差。未能纠正随机性状进化的影响可能导致高的假阳性率。诸如PGLS的方法校正了随机性状进化下预期的残差协方差，并产生更准确的关联统计检验。

b，我们使用16S拷贝数（＃16S）和回归系数与环境元数据（如干扰频率（β））之间的假设检验来说明PGLS。在检验性状之间的关联时应该使用系统发育广义最小二乘法(Phylogenetic generalized least squares, PGLS)，包括性状量，例如来自丰度的回归系数：元数据关联。为了实施PGLS，需要假设一个性状进化模型。我们用估计Blomberg’K 来系统发育信号; K> 1的值表明亲属在布朗运动模型下比预期更相似。然后我们使用来自R包phytools的corBlomberg函数来控制残余协方差结构。无论所使用的性状是否已知或通过祖先状态重建估算，都应使用PGLS。来自F检验的P值表明，关联的重要性取决于是否考虑了在随机性状进化下预期的协方差结构。

系统发育比较方法

PCMs在比较微生物的多种性状时使用。由于来自共同祖先的遗传，密切相关的生物通常具有相似的性状; 在不同微生物中的性状之间的依赖性可以影响（性状：性状）以及（栖息地：性状）的关联统计检验。

例如，16S拷贝数（性状）和pH偏好（栖息地）之间的关联可以通过16S拷贝数与1000种微生物的pH偏好测量值之间的相关性找到（图1a）。然而，如果被调查的分类群由一组具有低16S拷贝数和低pH偏好的密切相关的酸杆菌和一组具有高16S拷贝数和高pH偏好的密切相关的梭杆菌组成，则该组合的重要性可能被夸大。直观地，这些性状的系统发育信号减少了我们的样本量，因为观察到的性状代表来自两个谱系的样品，而不是1,000个独立样品。使用PCM进行性状关联的鲁棒检验（图1b）

广义最小二乘法（GLS）可以在进行回归分析时控制观测变量之间的依赖性。在GLS中，残差（预测值和观测值之间的差异）预期是协变的，并且协方差矩阵用于修正最小二乘计算。随机进化产生近亲，由于在共同祖先期间获得的共享变异，其观察到的性状将是协变的。系统发育感知的广义最小二乘（PGLS;图1），一个实现了具有由进化模型定义的残差协方差的GLS的用于 性状：性状 和 性状：栖息地 关联分析的工具。

针对没有系统发育信号的零模型（null model），PCM中常见的第一步是去估计和检验系统发育信号（alternative model）。 Pagel的λ参数或Blomberg的K参数通常用于系统发育信号假设检验的统计量（test statistics）。对于PGLS，必须假想一个进化模型; 布朗运动 （Brownian motion），即从树的根节点的祖先值到叶子节点的性状值的分支默认是以随机游走（random walk）进行。进化模型定义了残差的协方差矩阵（图1b）。在布朗运动进化模型下，两个物种的性状值残差之间的协方差与共享进化历史的数量成正比; 即使在随机进化的零模型下，更密切相关的物种也具有更密切相关的性状。对于更复杂的进化模型，可以联合估计进化模型和回归系数的参数。

PCMs涉及扩展到许多统计检验。在没有系统发育信号的前提下，检测细菌孢子的体积是否小于子细胞的体积，这将是一个很简单的配对t检验。系统发育感知的配对t检验（phylogenetic paired t-test）便是考虑了系统发育信号后的，配对t检验的衍生物。在比较性状时，存在许多性状进化模型，系统发育信号的度量以及控制系统发育信号的方法。最近学术版的现代PCM提供了对当前研究领域和方向的综述。

PCMs不常用于微生物组研究，尽管最近的一项研究采用PCMs来鉴定与人类肠道定殖有关的16S基因（性状：栖息地）。在性状检验中未能校正系统发育依赖性，性状：性状，性状：栖息地之间的关联可产生高假阳性率（图1b）。为了修改这一点，我们建议研究人员熟悉并使用PCMs。许多方法可以通过R包，ape , phangorn，phytools，picante，caper，Geiger和phylolm实现。在补充在线教程中，我们使用这些包来模拟性状进化，检验性状之间的关联，并说明这些方法对水平基因转移(HGT)的敏感性。

祖先状态重建

估计或重建祖先性状值有助于对未被描述的物种中的性状值进行插补（推测填补imputation），并确定出现主要性状差异的谱系。在微生物学中，祖先状态重建通常用于使用一组参考基因组估计现存群落的遗传和代谢谱。在微生物组研究中，这通常使用PICRUSt进行，PICRUSt使用祖先状态重建来估算性状值，例如编码糖苷水解酶活性的基因，用于性状未知的分类群。

PICRUSt 操纵 16S基因系统发育，连接测序的基因组和环境序列。测序基因组中的性状信息用于推断祖先性状谱。然后使用祖先谱来推测（impute）环境样本中微生物的概况。通过添加样品中操作分类单位（OTU）丰度的乘积及其相应的预测谱来估计输入样品的预测宏基因组谱。由于PICRUSt在很大程度上依赖于参考数据库和可用的测序基因组，因此在已知很少或没有基因组数据的环境中性能较差。相反，PICRUSt预测人类粪便样品中全基因组鸟枪测序谱，Spearman相关系数 r2> 0.9（参考文献38）。

祖先状态重建背后的方法与PCM有关，因为它们都需要进化模型。而在分析中主要有三类算法用于进化推断：最大简约性（maximum parsimony），极大似然（maximum likelihood）和贝叶斯推断（Bayesian inference）。最大简约性通过最小化祖先和当前后代之间的性状改变的次数来重建祖先状态。这种方法假设性状变化（比如16S序列变化）缓慢，并且不考虑快速进化的情景。此外，最大简约性将所有分支视为相同，并最小化每个分支上的改变次数; 这种计算可能是有问题的，特别是如果不是所有物种都被观察到的话。极大似然和贝叶斯推断通过将显式的进化模型（例如沿着树的性状进化的布朗运动模型）整合到祖先状态的估计中来改进最大简约性。这些方法并没有假设改变次数很少，而是可以考虑变化更频繁发生的一些情况 - 例如，假设同义替换比非同义替换更频繁 - 并且在给定估计的系统发育的情况下将参数拟合到这些模型里。然而，极大似然法通常低估了单个分支内的变化数量，并且可能产生次优解，特别是如果进化速率在整个系统发育中发生变化时。贝叶斯方法可以通过对可能的进化树的深度采样来计算进化参数，并评估更复杂的进化模型，这些模型可以考虑到非均匀进化速率的情况。虽然贝叶斯方法可以产生比最大简约性或极大似然性更准确的结果，但它们在计算上开销很高。因此，PICRUSt使用最大简约性或极大似然估计微生物祖先状态。至于PCMs，HGT可以快速混淆祖先状态的估计（参见补充在线教程），这些方法在微生物数据集的应用应该考虑所感兴趣的基因家族的观测到的转移速率。

图2|系统发育定义了群落生态数据的几何形状

就像球体定义GPS数据的几何形状一样

a，更改变量可以允许更复杂拓扑的自然描述。球形地球表示球面坐标。系统发育变量使用树作为构建与系统发育性状相对应的坐标的支架。 Phylofactorization构建对比组G1和G2的坐标，由边分开，其中出现诸如飞行之类的性状。

b，两点之间的默认路径是直线，但球体上更有意义的路径是测地线 - 即沿球体表面的最短路径。同样，诸如UniFrac之类的系统发育感知距离定义了进化路径以及它们在一个群落与另一个群落之间的距离。

c，PhILR构建对比姐妹分支之间的坐标。

d，可能的系统发育变量和距离的空间是无限大的。如此处所示，远距离进化枝之间的比率是可行的，但目前尚未使用的系统发育变量。研究人员应该考虑新变量和距离的生物可解释性，以及它们为未来研究提供信息的能力。图标（袋鼠，鸭嘴兽，蜥蜴，伸出的蛇，鹰，猫头鹰和鹤）由Freepik从www.flaticon.com制作。

系统发育变量分析

地球表面上的位置可以用三个笛卡尔坐标（xyz）来描述，但是更自然地是使用两个球坐标（纬度和经度）描述。类似于球体，系统发育让人联想到自然坐标系（natural coordinate）。系统发育变量用于减少群落生态数据的维度，简化距离计算，并描述群落中有意义的性状和变化方向（图2）。

我们用“系统发育变量”这个术语来描述使用系统发育中的性状构建的变量，对比和总结系统发育树中物种的数据（图2）。变量和距离是相关的，但包含不同的信息：比如说一个城市在东方并不表示它有多远，并说它是80公里以外却不知道它是哪个方向。通过系统发育变量描述方向（图2a），通过距离测量变化幅度（图2b）。系统发育变量包括taxa丰度，多样性测度，沿所有边的丰度差异，进化枝之间丰度的差异（图2a-d）等。

系统发育变量通过将数据的维度减少到携带生物信息的一些变量来简化微生物组数据集。如果一些单系分支在环境梯度上解释了微生物组数据集的大部分变异，那么每个分支的成员之间可能存在一些共同性状，这些性状决定了丰度并且是沿着环境梯度观察到的群落组成变化的基础。

可能的系统发育变量集是无限大的。因此，研究人员在选择系统发育变量时必须慎重考虑 - 对于进一步研究有哪些重要的变化方向？沿着系统发育变量（例如α多样性）的方向进行的群落变化不一定传达有用的生物信息或容易建议未来的研究方向。系统发育变量分析中的两个常见挑战可以指导它们的选择和发展：统计依赖性和生物可解释性。

统计独立性或充分表征的依赖性促进了鲁棒的多变量统计分析。例如，当检验种丰度和环境元数据之间的关联，并重复属，科，目，纲和门的过程时，所分析的变量具有嵌套依赖性：如果一个分类级丰度增加，那么则牵一发而动全身，其他分类级也会跟着变化，所有它的更高分类群的丰度也都增加。再例如，如果每个序列都是新物种，则n个序列的Shannon多样性将是H = log（n），并且物种丰富度和均匀度（richness and evenness）将是相关的，会相应的改变。在进行多重假设检验（multiple hypothesis tests）时，未能考虑系统发育变量之间的依赖性会增加错误率。

具有明确生物学解释的系统发育变量可以促进未来的研究设计和理论发展。单系分支的丰度变化表明可遗传的性状可驱使丰度的变化; 未来的实验可以集中在分支上，以寻找可能的功能性生态性状。在宏观生态学中，理论证明了各种多样性测度（diversity metrics）作为灭绝率，岛屿生物地理过程，生态系统稳定性和保护目标的代名词（proxies）的合理性。系统发育变量的理论依据将系统发育变量（例如，多样性和元数据之间的关联）的分析与实验设计和生物学理论联系起来。

两个最近开发的方法 - 系统发育等距对数比（PhILR，phylogenetic isometric logratio ）和phylofactorization- 说明了系统发育变量分析的挑战。鉴于测序序列计数数据（count data）的组成型特征（compositional nature），两种方法通过系统发育中两个进化枝之间的丰度的平均对数比来构建变量。 PhILR变量计算的是姐妹进化枝之间的差异（图2c），而phylofactorization是迭代构建变量，测量由树中的边分开的进化枝之间的差异（例如图2a，d）。 PhILR坐标的变化表明一个性状是区分姐妹分支，而phylofactorization的坐标变化表明沿着识别的边出现了一个性状。在这两种方法中，系统发育变量和元数据之间的显著关联促进未来的工作去比较两个进化枝的基因组以寻找功能性状。 PhILR可用于比较姐妹进化枝（例如，胎盘哺乳动物与有袋动物，或鸟类与鳄鱼），而phylofactorization将比较由临界边分开的进化枝（例如，鸟类与非鸟类）。在补充的在线教程中，我们用系统发育变量分析说明了这两种方法，展示了如何构建这些变量，并将它们与EdgePCA进行比较，EdgePCA是一种对与边缘丰度差异相对应的变量进行主成分分析的方法。为了说明这些方法，我们分析了一个模拟数据集，其中rRNA基因拷贝数驱动与土壤中干扰频率的关联，并解释结果。

分析系统发育变量的目的是确定微生物组数据有意义的变化方向。主成分分析确定数据变异的主要方向/轴，而系统发育变量识别微生物组数据的变化方向，解释了群落组成的变异，并对灭绝风险，需要培养哪些微生物，抑或要比较的基因组等有影响。

图3 |系统发育感知距离

a，物种丰度的热图，红色表示高丰度，黄色表示不同环境下的丰度低。进化历史用系统发育树表示，环境A和环境B之间的主要差异由进化枝A和进化枝B的丰度驱动.

b，虽然变量包含每个样品的信息，但距离涉及两个样本。绘图是样本之间成对的UniFrac距离; 与来自环境A的样本之间的距离或来自环境B的样本之间的距离相比，来自环境A的样本与来自环境B的样本之间的距离更大.

c，来自环境A的样本与所有其他样本之间的UniFrac距离说明距离用于样本地点分类如何有用。 UniFrac的距离高表明两个群落之间的系统发育新颖性或系统发育代表性的变化大。图标（袋鼠，鸭嘴兽，蜥蜴，伸出的蛇，鹰，猫头鹰和鹤）由Freepik从www.flaticon.com制作

使用系统发育感知距离

量化不同物种之间以及包含这些物种的不同群落之间的差异可以促进元数据的准确分类（例如患者是否患有疾病），样本聚类以及群落功能的推断。森林中的树木将木材中的碳固定在一起，而草则不会。因此，测量含有树木的群落和含有草的群落之间的距离可能表明森林和草原的生态系统功能存在差异。对于微生物世界而言，驱动生态系统功能的性状通常是未知的，但疾病状态的准确分类可能对人类健康产生重大影响，并且当类似于木质生物量的遗传性状构成栖息地关联时，将系统发育纳入距离测量可以有助于分类（图3）。系统发育感知距离将数据集（图3a）转换为样本之间的距离矩阵（图3b），可用于对样本进行分类（图3c）。

用于微生物组数据的系统发育感知分析的最广泛使用的方法之一是分析样品之间的UniFrac距离 。与标准欧式（Euclidean） 和Bray–Curtis距离相比，UniFrac距离被用作群落之间更具生物学意义的距离。 UniFrac背后的直觉以及大多数系统发育感知的距离是，包含更多系统发育不同物种的群落与物种更密切相关的群落更为不同。将系统发育距离纳入发生群落变化中可以更好地量化群落之间的功能差异。

已经有许多UniFrac的衍生或者扩展被探索，目的是控制计数数据中的统计误差（statistical artifacts）并调整物种丰度值在UniFrac距离中的重要性。如果计数在物种之间随机分布，则具有更多物种的进化枝将具有更高的总计数方差，因此对UniFrac距离的影响大于物种较少的进化枝。为了弥补这种影响，方差校正加权（VAW）的-UniFrac方法可以稳定UniFrac距离的方差。VAW-UniFrac对广义的UniFrac距离进行了扩展，该距离包含一个用于增加或减少丰度在群落之间距离的重要性的可调参数。

还有许多其他系统发育感知的距离指标，例如索伦森指数（Sorensen’s index），Rao’s D和Rao’s H，它们的差异在于如何结合进化信息。此外，诸如线性回归这样的标准统计方法可以纳入其中，以惩罚近亲之间的差异。系统发育是许多变量的支架，可以作为许多有用的距离度量的基础。那么，在所有可能的距离测度中 ，哪一种距离是微生物学家感兴趣的呢？

构建和使用系统发育感知距离有两个非排他性的研究目标类别：改进样本点分类或可视化，并提供具有生物学意义的群落差异。对这些研究目标的明确认识可以增加对系统发育距离分析的影响。

如果样本点分类或可视化是实验的目标，研究人员可能倾向于搜索一个可能的距离空间，直到找到看起来最好的一个距离，而不管这种距离的生物学意义。然而，搜索太多距离有可能是冒着在过度挖掘数据（risks dredging the data）并只呈现统计学显著的模式（pattern）的假阳性风险，因为这些模式你可能是通过检验多个候选者而并未进行多重假设检验校正而得到的。而要纠正这种多重假设检验将面临同样的挑战，以为我们需要解决在分析多个系统发育变量时出现的检验之间的依赖性问题（传统的多重假设检验指不考虑多个假设之间的依赖性，假设彼此独立）。

虽然许多现有距离可以成功地分类，以及在一系列的地点类别和临床变量中也可以可视地分离样本，但发现的差异的生物学意义通常不清楚。距离越远表明一个群落的生物修复难度越大吗？更远的距离是否意味着生态系统功能或患者发病率的差异更大？在UniFrac距离较大的情况下，应该进行哪些后续实验以更好地了解群落差异的生化和微生物原因？

构建新的系统发育感知距离及其在修正后的统计方法中的应用应考虑相对于现有方法的性能增益（performance gains），以及它们是否对发现的差异提供了新的解释。仔细证明新的距离可以改善结果的生物学解释。例如，宏观生态学家讨论了β多样性如何用于保护。此类讨论可以改进对现有和新开发的系统发育感知距离的解释，并帮助研究人员了解群落之间距离大或者小的影响。此外，高质量的树对于揭示生态相关模式至关重要。与PCMs和系统发育变量一样，系统发育感知距离受益于对生态和进化模型的明确考虑，以帮助对其结果进行生物学解释。

方框2 |对微生物组数据进行系统发育分析的挑战

Horizontal gene transfer 水平基因转移。 HGT破坏了基因进化历史之间的相关性，并提出了关于哪些基因树用于系统发育分析的重要问题。虽然16S基因树与微生物中的大部分基因组含量相关，但重要的水平传播基因如β-内酰胺酶具有与16S不同的系统发育。分析β-内酰胺酶基因树将允许分析β-内酰胺酶性状 - 这样的树可能适合于研究环境中抗生素抗性基因的组成。我们将在本文的“系统发育分析的挑战”一节中进一步讨论这一点。

Phylogenetic inference 系统发育推断。 16S rRNA基因树是最常用的，但其他基因，如β-内酰胺酶基因，可用于构建系统发育。无论什么基因，系统发育推断都是对进化历史的估计，物种采样规模越大越均匀，则分类估计最准确。不均匀的分类单元抽样会导致错误的系统发育，导致类似的性状被误解为同源性。使用许多种蜥蜴的皮肤和骨骼结构进行系统发育重建，一种鸟类和一种蝙蝠可能错误地估计鸟类和蝙蝠是姊妹类群，而包括哺乳动物在内的更完整的分类群采样可能正确地将蝙蝠与哺乳动物分组。在每个研究地点从头构建进化树，每个地区的分类单元抽样有限，可能会产生许多错误的树，这些树难以在不同的研究中进行比较。从常用基因构建的全局共识的进化树，均匀的分类单元采样和用于向现有树添加新序列的标准化方法，可以确保研究人员对微生物进化历史的进行可比较的推断是建立在同一个合理，准确的进化骨架上。

Ancestral state reconstruction 祖先状态重建。与系统发育推断一样，稀疏的分类物种抽样可以增加祖先状态重建的错误率。诸如PICRUSt之类的方法利用来自相对良好采样的环境（例如人类微生物组）的生物体的基因组和性状，可能在不太良好采样的环境中对微生物具有高错误率。

Vast number of species 海量的物种。最近的一项研究估计有超过一万亿微生物物种。虽然存在宏观生物的大型数据集，但微生物生态学中物种丰富数据集的规律性以及收集许多样本的容易性需要特别考虑，通常为较小数据集开发的方法的计算成本，可视化和解释。并行化，强调共同知识或重要性的谱系，微生物学家之间对系统发育的共识，以及通过折叠进化枝简化的系统发育表示，可以使研究人员对微生物大数据进行彻底的系统发育分析。

Evolutionary model for microorganisms 微生物的进化模型。没有微生物化石记录。在宏观生态学中，化石记录用于校准系统发育推理和祖先状态重建所必需的进化速率。虽然我们知道不同物种中的不同基因具有不同的突变率，但微生物进化模型的校准和验证仍然是一个开放的研究领域。正确的进化模型可以产生准确的效应大小和不确定性的测量（显著性，置信区间等），确保系统发育结构化数据分析推论的准确性和可重复性。

系统发育分析的挑战

系统发育结构化的数据分析充满了挑战，包括HGT，用哪一种基因树，在系统发育推断误差的敏感性和鲁棒性，以及准确考虑该用哪种生态和进化模型。在这里，我们讨论系统发育分析的更广泛的挑战; 对于与微生物和微生物组数据集特别相关的挑战，见框2. 微生物基因组之间的HGT使系统发育树中捕获的垂直传播的进化故事复杂化。 HGT提出了一个问题，即使用哪种系统发育以及系统发育对于研究问题的信息量准确性如何。对于PCMs，HGT可能导致不正确的校正和校准不良的统计检验（见补充材料在线教程）。驱动数据变异的主要性状的HGT可降低系统发育变量或所用距离的适当性。有利的是选择对HGT不敏感的基因家族来推断系统发育。研究已经根据功能和生态性状评估了HGT的可能性，为此任务提供了指导。也许在整个生命树中没有完全HGT-free的基因，包括16S。在系统发育推断中使用多个基因可以最小化HGT的负面影响，并在选定的分类群中揭示受HGT影响的基因。计算工具可用于评估基于物种或基因树调节的假定HGT事件的概率。基因组的状况，性状序列和非典型的同源性搜索结果的探索也有利于跟踪HGTs。

HGT不会使微生物组数据的系统发育感知分析无效。对于功能性状的HGT我们可以通过系统发育感知分析，假设HGT只和很少的系统发育信号有强烈的影响。如果phylofactorization识别出与抗生素暴露相关的树的异常大量的叶子节点，则HGT可能驱动数据的变化并且可以通过比较所鉴定的系统发育因子中的基因组来进一步检验。尽管如此，在分析系统发育结构数据时需要考虑HGT。许多方法对功能性状水平转移的敏感性目前尚未得到充分研究。最后，所有方法都面临着生物可解释性和提高我们对微生物系统知识的挑战。

为此，新方法应明确考虑生态和进化模型，以了解性状如何发展并驱动数据中的模式。一项研究模拟了树上的性状进化，并将PGLS与系统发育特征向量回归方法进行了比较，这些方法使用来自系统发育距离矩阵的特征向量作为自变量，而并没有对应一种明确的进化模型。该研究发现PGLS产生了更可靠和更好校准的统计结果。考虑方法开发中的进化和群体遗传模型可以促进对系统发育感知分析表现良好的假设的准确理解，并在潜在生物过程的背景下通知对发现的解释。系统发育分析的可解释性很重要。新颖的和复杂的方法开发，研究人员应该了解机器学习和人类理解之间的权衡：前者可能会产生短期更精确的预测，而后者产生的理论，可以产生在长期更准确，更普及的预测术语。

讨论

微生物的共同祖先可能是我们数据中混杂变异的来源，或者是我们做出推论的脚手架。有许多现有的和新兴的方法用于根据进化分析微生物组数据集，选择正确的方法需要精确的研究问题表述（表1）。

首先，决定使用哪棵树。通常情况下，微生物研究使用细菌和古细菌的16S树和真核微生物18S树，但有一个系统发生的每个基因和一些问题，更好地与其他基因树分析。获得的系统发育将是一种估计，系统发育推断的不确定性可以转化为下游系统发育结构数据分析的不确定性。

如果研究问题使用性状作为因变量，则系统发育可能是混杂变异的来源。诸如PGLS之类的PCMs校正了在树的进化的空模型下人们期望的性状之间的依赖性。为了研究历史性状值或主要性状差异产生的边，需要祖先状态重建。如果检验估算性状之间的关联，研究人员需要将祖先状态重建与缺失性状的归因与PCM相结合，从而纠正混淆变异。为了简化群落组成的模式，系统发育是一种支架，可用于产生有生物学意义的变量和变化方向。变量的选择应根据其捕获数据性状，统计依赖性及其生物可解释的能力来进行。

为了区分微生物群落的样本，系统发育可以定义样本之间的距离。通过重新定义距离，系统发育几乎可用于修改任何统计方法，但选择使用的距离应基于样本点分类的研究目标或差异的生物解释。微生物组数据的系统发育分析可以使研究人员对未分类的微生物进行分类，检验关于性状关联的进化假设或驱动栖息地关联的性状，并更好地了解微生物群落的差异以及它们如何随时间，空间和治疗方式而变化。根据进化论 ，分析微生物组数据有几种可能的方法。仔细考虑研究问题以及生态和进化假设，研究人员可以识别现有方法或生成解决其研究问题的新方法，并产生空前的、准确的和有生物学意义的见解。微生物序列的海量信息正在产生系统发育结构化数据，如果有正确的工具，可以加速我们对微生物群落结构和功能的理解。

原文链接：https://doi.org/10.1038/s41564-018-0156-0

https://www.nature.com/articles/s41564-0

猜你喜欢

系列教程：微生物组入门 Biostar 微生物组  宏基因组

专业技能：学术图表 高分文章 生信宝典 不可或缺的人

一文读懂：宏基因组 寄生虫益处 进化树

必备技能：提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读 分析流程 统计绘图

在线工具：16S预测培养基 生信绘图

科研经验：云笔记  云协作 公众号

编程模板: Shell  R Perl

生物科普:  肠道细菌 人体上的生命 生命大跃进  细胞暗战 人体奥秘  

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外3000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。PI请明示身份，另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。

学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”