分子系统发育分析中的贝叶斯方法

分子系统发育分析中的贝叶斯方法

引言

贝叶斯推断 (bayesian inference, BI) 于上个世纪90年代引入系统发育学。最初，贝叶斯推断引起了一波关于该方法优缺点的激烈讨论，特别是与最大似然法 (maximum likelihood, ML) 之类在系统发育学中具有稳固地位的方法之间的比较。其中，部分讨论围绕哲学或基础的问题展开，如何选择先验分布 (prior distribution) ？以及如何控制相关分析对这种选择的敏感性？后验概率 (posterior probability) 的意义是什么？如何与其它统计支持度 (如非参数自展 non-parametric bootstrap) 比较？

贝叶斯推断在系统发育基因组学和分子进化研究中的一系列问题上，很快便达到了实际应用的阶段。为此，许多工作致力于设计越来越复杂的模型和开发基于马尔可夫链蒙特卡罗 (Markov Chain Monte Carlo, MCMC) 的高效算法。基于对当前研究的实际影响，以及贝叶斯推断与进化遗传学的相关性，它已经可以被学界更好地理解和接受。事实证明，贝叶斯推断的使用极大地改变了我们对系统发育基因组学相关模型作用的观点，同时提供了传统方法无法企及的新机会。但不可否认的是，一些实际应用也指出了纯贝叶斯立场的一些弱点和局限性。

系统发育学中的贝叶斯推断

概率推断的一般原则

现在我们考虑一个简单的系统发育问题：用单个基因 (如核糖体RNA) 的多序列比对来推断一个分类群 (taxon) 的系统发育。假设该基因满足一个简单的单参数序列进化模型，如Kimura模型 (仅依赖于转换/颠换比transition-tranversion rate ratio)。

用来表示多序列比对，表示未知的树拓扑结构，表示分枝长度，表示转换/颠换比。

而后我们用似然 (likelihood) 来定义：具有转换/颠换比的序列进化过程，沿着分枝长度为的树运行，在树的顶端产生核苷酸序列的概率。此概率可以写成：

在这个系统发育问题中，假设序列中各个位点都是相互独立演化的。因此，该概率可被改写为各个位点上的概率乘积形式 (独立事件与概率乘法定理)。设表示第个 (，表示所有对齐位点的数量) 位点上的序列比对，所以有：

在最大似然法中，估计树的工作原理如下：

首先，对于给定的树，通过调整分枝长度和替代模型参数使似然最大化：

这定义了相应树的似然分数。然后，搜索下一个能提升似然分数的树，也就是说，目标是寻找一个树拓扑能够让似然分数最大化，

需要注意的是，似然分数是在多个未知参数条件下，联合最大化的，其中有我们直接感兴趣的树拓扑，也有那些我们不直接感兴趣，但多序列比对的产生概率有影响的未知参数 (称为冗余参数：分枝长度和转换/颠换比)。因此，树拓扑仅基于该树拓扑结构下进化过程的所有未知方面的最佳场景进行评分，而不考虑这些未知冗余参数的其它有细微差别的替代配置。这是区分最大似然和贝叶斯推理的一个重点。

贝叶斯方法

对于当前的例子，贝叶斯方法又是如何进行推断的呢？

首先，我们必须为树拓扑、分枝长度和替换模型参数，分别定义一个先验分布，这些先验可以表示为，和。注意，是定义在 (离散的) 树拓扑空间上的一种概率；另一方面，由于和是连续型的参数，因此和不是概率，而是概率密度。通常我们会设定简单的先验分布，如：一个具有均匀先验概率的 (树空间中的所有树等概率)，一个符合指数分布 (均值为) 的分枝长度，和在到之间均匀分布的转换/颠换比。

然后，对于给定的树拓扑，将序列比对的条件概率——在连续参数和的所有可能取值条件下，的均值，定义为似然：

似然，即给定一个树后，产生序列比对的概率。

公式 (1) 中描述的实际上是一种加权平均 (通过这两个参数分量和 的先验分布 和 加权。

最后根据贝叶斯公式，我们可以定义树拓扑的后验概率：

其中，分母是多序列比对的边缘概率 (或边缘似然)，起到对后验概率标准化的作用，根据全概率公式：

标准化作用，即让的积分等于1。

此外多数情况下，归一化因子不是必要的，因此贝叶斯定理可简化为：

这表明：树拓扑的后验概率，与其先验概率与似然 (数据贡献的证据权重) 的乘积成正比。

: 正比于符号

因此，如果代表了我们在看到数据之前关于哪些树可能是正确的信念状态，那么贝叶斯定理形式化了如何在看到数据时更新我们对正确树拓扑的信念，并生成我们的后验信念。对于非常大的数据集，后验通常会集中在一个树拓扑上。就信念而言，考虑到数据和模型，这就等于几乎完全确定了这棵树就代表了正确的系统发育。另一方面，对于较小的数据集，后验可能会分散于多个树，这代表了根据数据，我们对感兴趣分类群的系统发育，仍保留有不确定性。

上述方程可以等价地重写，通过使用贝叶斯定理定义树拓扑和连续参数上的联合后验：

这个公式 (3) 更为普遍，因为它包含了所有未知参数，包括树拓扑、分枝长度、替代模型的参数。

仍然假定我们对树拓扑感兴趣，接下来聚焦于树拓扑的边缘后验概率，它可以通过对和积分而得到：

这里所定义的树拓扑的边缘后验概率与公式 (2) 是等价的。如果我们对估计转换/颠换比感兴趣，那么我们将考虑上的边缘后验分布，它在所有可能的树拓扑和分枝长度上求和：

最后，贝叶斯推断总是简化为计算我们想知道的后验概率，给定现有的数据证据和生成模型的结构假设，整合所有其它未知参数 (求平均)。

比较公式 (4) 和公式 (5)，可加深对积分和求和之间关系的理解。

使用蒙特卡罗的实践贝叶斯推断

上述方程涉及树空间中大量的树拓扑结构，对于给定的树，需要求得分枝长度的积分和值的积分。一般来说，这些积分是难以获得解析解的，也很很难用足够的精度进行数值估计。作为一种替代方案，贝叶斯推断常使用蒙特卡罗方法来实现。蒙特卡罗的一般目标是设计针对感兴趣的概率分布的随机抽样算法，这里也就是公式 (3) 表示的联合后验分布。

目前最常用的算法是马尔可夫链蒙特卡罗 (Markov chain Monte CarloMCMC) 方法。MCMC的思想是在模型的参数空间中实现一个随机游走 (random walk)，以使得参数配置以正比于其后验概率的频率被访问。算法运行足够长时间，将会产生后验分布的样本：

其中 ，为样本数量。对于MCMC，样本之间并不是独立的 (连续样本通常相关，由马尔可夫链模型定义)。此外，它们只是目标后验分布的渐近分布。在实践中，这意味着，从任意参数配置开始，马尔可夫链只有在预烧期 (burn-in period) 后才达到稳定状态，而只有在达到这种稳定状态后，样本才能被认为是从后验分布中抽取的。

一旦得到了大量的样本后，联合后验概率上的任何边缘概率都可以简单地用样本上平均值来近似。例如，样本中给定的树拓扑的频率将成为后验概率的蒙特卡洛估计。更一般的说，在样本中发现一组给定类群呈单系的频率是该单系分支后验概率的蒙特卡洛估计。因此，对分析进行总结的一种方便方法是获得蒙特卡罗收集的所有树在多数原则 (majority-rule) 下的一致树 (consensus tree)，并用其后验概率支持的蒙特卡罗估计标记每个分枝。

另一方面，如果我们感兴趣的是一些连续参数，例如转换/颠换比，那么蒙特卡罗收集的值的直方图表示我们对这个参数的后验概率密度的估计，还可计算出95%的可信区间。

贝叶斯推断的一些重要属性

对不确定求平均

首先，对一些感兴趣的参数进行贝叶斯推断时，总是平均所有其它冗余参数的不确定性 [公式 (4)和(5)]。这一点我们可从如何处理MCMC的输出中获得一些理解。例如，一个给定分枝后验概率的蒙特卡洛估计，只是这个分枝在所有蒙特卡洛样本树中出现的频率。重要的是，所有含有该特定分枝的样本树，可能在许多其它方面有所不同——包括分枝的长度、模型参数值，还有树中其它地方存在的其它分枝。通过这一点可以看到，在我们评估一个给定类群是单系的可能性时，我们对问题所有其它方面的许多可能结果进行了平均。这与最大似然形成鲜明对比，后者在决定其点估计时，只在某个冗余参数的单一配置上下注 (最大似然法，当给定树时，使似然最大化的其它冗余参数的取值是单一的)。对不确定性的平均将导致更稳健的推断，特别是当有许多冗余参数时。

蒙特卡罗 vs. 解析积分

上面讨论的这一点也显示了蒙特卡罗与积分之间的关键关系，即：从树拓扑和其它参数的联合后验分布上抽样，然后丢弃所有参数，只保留树拓扑，相当于从树拓扑的边缘后验分布上抽样树拓扑。换句话说，蒙特卡罗自动实现了公式 (4)，而且还不用显式地计算这个积分。数学上这种有趣的观察有着重要的实际意义：每当一个似然是许多冗余变量上的复杂积分时，总是可以明确地从冗余变量中抽样，而不用明确地计算这个积分。这种参数扩展（或数据增强）的方法可以实现更广泛的模型，尤其是那些无法通过数值积分而实现的模型。

先验分布

如上所述，对冗余参数的可能性进行平均，预计将导致更稳健的推断。另一方面，这个平均值继承了一个特定的先验分布。更一般地说，基于与先验概率成正比的后验概率做决策，提出了如何选择先验分布，以及统计推断如何依赖于这一选择的问题。先验分布的选择可能是贝叶斯推理中最重要的问题，不同的选择会在概念上的和实际上的产生不同的结果。关于这个问题有大量的文献，许多问题仍然悬而未决。问题是复杂的，因为事实上从不同的哲学角度有非常不同的方法定义先验分布，不同的方法也导致不同意义的后验概率，以及不同的操作属性。

在模型设计中的灵活性

使用机制启发的和分层的先验分布，结合一般的蒙特卡罗方法，具有重要的实际意义。事实上，它使得设计分层模型，连接多个层次的过程和集成多个经验数据源成为可能。在这方面，在过去十年中出现了许多重要的进展。

贝叶斯非参数位点异质性模型

位点间的差异

在大的进化尺度上，多重序列比对最显著的特征之一是位点间的差异和不同的保守性模式。这在核酸和氨基酸水平上都能看的到。原因是明确的：在长期进化周期中保守的序列几乎肯定经历了强净化选择。然而，选择是高度上下文特异的，换句话说，选择不会对所有核酸/氨基酸位点“一视同仁”。因此，在一个基因的不同位点上，选择无论是在整体强度上，还是在偏好核苷酸或氨基酸的性质上，都可能是不同的。系统发育重建通常是使用那些可以在更大系统发育尺度上，可靠地对齐的基因/基因区域上进行的，这一事实进一步放大了这个问题。为了保证”数据中的序列变化只是由点突变引起的”这一假设的有效性，选择可靠的序列对齐结果是必要的。这一假设是目前在系统发育学中几乎所有模型都默认的假设。然而，这样做会对那些结构高度保守的基因/基因区域产生选择偏差。反过来，这意味着一个给定的对齐位点，坐落在一个非常特定的生化环境中，诱导强位点特异的纯化选择，以保持大分子的构象稳定性，而这种选择又是长期稳定的。

就由此产生的序列进化过程而言，这种跨位点的选择性约束将转换为跨位点的替代速率和替代模式的差异。这种广泛的差异将如何影响系统发育的估计？我们应该明确地在系统发育重建的模型中体现不同位点的这种差异呢？还是捕捉所有位点的平均替代过程就足够了？如果显式建模对系统发育的准确性很重要，那么我们如何设计模型来准确地捕捉跨位点的替代速率和替代模式的分布呢？这些都是最近系统发育组学中的重要问题。

跨位点的速率的变化：一个参数化的随机效应模型

1994年，Ziheng Yang最早在最大似然的框架下，实现了对跨位点替代速率差异的建模。随后，有参随机效应方法没有进行重大改动，而被移植到了贝叶斯推断中。在解决如何对跨位点替代模式异质性建模，这一具有挑战性的问题之前，详细了解这种方法的概念结构是很有用的。

Yang引入的跨位点速率模型的基本思想是将特定位点的相对速率作为随机变量，该随机变量假定服从伽马分布，均值为 (因为是相对速率)，以及依赖 (伽马分布) 形状参数的未知方差。在数学上，位点的似然也就是在所有可能的替代速率上积分：

其中是伽马分布的概率密度函数，为了符号上的简便，我们用来代表模型的所有全局参数 (除以外的，如树拓扑、分枝长度和替代模型的参数)。实践中，这个积分是也难以处理的，标准方法是通过少量离散的类速率值，来数值估计它，通常以伽马分布的分位数为中心 (一般):

最后，我们可以针对所有位点得到概率：

这给出了整个序列比对的似然。该似然可以关于和进行最大化。或者在贝叶斯方法中，可以为模型参数和定义的先验分布，然后从联合后验分布中抽样：

首先，位点特异的速率被集成在一个分布上，并对分布本身进行估计 (特别是其方差，)。在最大似然法的思想里，对所有位点的所有速率的进行似然最大化，在原则上是可能的。然而，除非分类群的数量非常大，同时树也非常长，否则这将导致过度拟合。位点特异的速率估计将存在较大的随机误差，从而导致对树拓扑的进一步估计出现误差。此外，估计的速率方差将会大于真实的速率方差，因为速率方差的估计还包含速率估计误差的额外贡献。

相比之下，将速率作为随机效应来处理，会自动折扣这种额外的采样方差，并通过，公平地估计不同位点的真实速率方差。该随机效应模型在树的拓扑结构和全局参数的估计方面也取得了更高的精度。这是关于随机效应模型的一个重要观点：在将模型与数据拟合后，随机效应的值仍然有很大的不确定性，然而，在许多情况下，我们仍然会对它们的分布和模型的全局参数 (特别是树拓扑) 实现渐近一致的估计。

其次，上面考虑的伽马速率模型是一个有参数的随机效应模型，因为我们假设位点之间的速率分布属于一个预先指定的参数族 (这里是一个伽马分布)。没有什么能保证这一假设在实践中不会被违反。跨位点速率的真实分布可能是任意的，而这种真实分布和模型所假设的分布之间的不匹配，原则上可能对系统发育估计的准确性产生不可忽略的影响。一般来说，伽马分布（或一定比例不变位点和伽马分布的混合）将提供一个足够好的描述跨位点速率分布的描述。值得注意的是，现在已经有替代方法来放宽有关速率的假设，其中一些与下面考虑替代模式异质性的方法相似。

跨位点的氨基酸偏好：非参数化模型

不仅是速率会呈现跨位点差异，核酸/氨基酸的替代模式也会呈现跨位点差异 (下面重点讨论氨基酸序列)。考虑到不同位点表现差异的主要因素是选择，需要对这种跨位点差异建模的最重要特征或许就是氨基酸偏好了 (作为氨基酸适应度的代理)。在数学上，位点特异的氨基酸偏好，可以通过演化过程中氨基酸的平衡频率所构成的20维向量来表示。该向量通常被称为氨基酸频率侧写 (amino acid frequency profile)。

与位点特异性速率的类比表明，我们应该将氨基酸侧写建模为位点特异的随机效应；而且，我们应该有一种方法，允许足够准确地估计氨基酸侧写在各位点间的真实分布。然而，位点特异性速率和位点特异性氨基酸侧写之间存在着重要的技术差异，因此用于速率的方法不能直接推广到目前的情况。首先，基于分位数的离散化方法不能很好地扩展到高维随机效应。另一个问题是，氨基酸频率侧写的真实分布是非常复杂的，也可能是多模态 (multimodal)，因此可能不能用任何已知的简单参数分布来很好地描述。

混合模型

一种替代方法是使用有限的混合模型mixture model。混合模型背后的基本原理是，在序列对齐位置上表现出来的氨基酸偏好模式的多样性，可能有望被少数几个的典型氨基酸频率侧写所表示 (例如疏水性，极性，带负电荷，芳香族等)。考虑个分量，每个分量都有自己的20维频率侧写和权重，那么位点的似然是所有混合分量的加权平均：

进而对于所有的位点有：

该似然依赖于频率侧写的集合和权重向量，以及统一用表示的模型中的其它参数。在最大似然法中，这种似然通常会通过、和使其最大化。或者，可以预先基于一个数据库来估计一组 (个) 频率侧写，并固定下来，然后在系统发育推理中使用。这就是所谓的经验混合模型empirical mixture models。

研究者花了很多时间来推导出可以在系统发育学中常规使用的经验混合模型。然而，到目前为止，这种方法已经产生了混杂的结果。一个主要问题是，为了获得良好的经验拟合和足够的系统发育精度，似乎需要的不同频率侧写的数量很高 (分量较大)，意味着氨基酸偏好的真实分布可能过于复杂或分散，无法用少量典型氨基酸频率侧写来描述。而且，如果允许大量的分量在模型中存在，将会快速地增加计算和统计上的挑战。从计算上讲，在每个位点上针对混合模型中所有的频率侧写求平均，对于较大的是难以实现的。从统计学上看，混合模型中的较多频率侧写会很快变得冗余，因为许多频率侧写只有小细节不同，通常会给出未知经验分布的基本等效近似，从而导致难以识别的模型。

然而，这些问题在贝叶斯框架中并不是那么关键，有两个不同的原因，与贝叶斯推理处理复杂模型的方式有关。首先，在贝叶斯的MCMC中，参数扩展可以用来避免在每个位点对所有频率侧写显式地求和。相反，我们可以显式地在MCMC过程中对混合模型中不同的频率侧写进行抽样。MCMC策略的复杂性对混合模型中的频率侧写数量相对不敏感。其次，贝叶斯方法可以通过平均所有可能模型配置的后验分布，来自动处理混合模型的冗余性。

非参数随机效应模型

实际上，我们可以在一个完全不同的指导思想下使用混合模型，即不需要保持分量的数量尽可能低 (代价是不正确地捕获生化特征侧写的真正经验多样性)，而可以使用非常冗余的混合模型。这样做会带来更大的灵活性。足够丰富的模型可以近似任意分布与任意精度，他们的冗余性并不重要，只要一个有效的MCMC能够通过平均代表性样本替代混合模型配置平滑掉这种冗余，所有这些可基本上实现对真实分布的等效近似。这就是贝叶斯非参数随机效应模型背后的基本思想。

非参数推断的原始目标是放宽关于真实分布应该先验地属于一个预先指定的参数族的假设。原则上，非参数方法应该给出任意随机效应分布的渐近一致的结果。在贝叶斯方法中，这是通过设计先验混合模型来实现的。狄利克雷过程Dirichlet process就是这样一个非参数化的先验。从技术上讲，狄利克雷过程通过实现无限混合，将冗余混合模型的想法推到极限。无限混合在所有可能分布的空间中都是密集的，因此，先验狄利克雷过程将考虑任何概率分布，当然包括跨位点的随机效应的真实分布。然后，将模型套在一个足够大的数据集上，将导致一个集中在真实分布附近的后验分布。最后，使用基于参数扩展的MCMC方法实现这一强大的非参数推断想法，原则上在任意（可能是多维）随机效应分布下可实现渐近一致性。

软件与平台

目前有大量的软件程序可用于进行系统发育相关的贝叶斯推断。这些程序通常有非常不同的具体目标或专业化：系统发育重建、分子年代测定、系统地理学和系统动力学、系统发育密码子模型、比较研究、基因树/物种树调和，或物种划分。

尽管目前向整合建模方向发展，系统发育基因组学中的应用研究，仍在使用更经典的超级矩阵范式，其中包含大量的单基因，同时假设它们是沿着同一物种树进行演化的。

有四个主要的软件程序被用于系统发育基因组学分析，特别是使用超级矩阵重建物种树：MrBayes，PhyloBayes，ExaBayes和P4。理想情况下，在模型设计中结合这些不同层次的单一实现将是非常有用的 (即允许同时在基因、位点和分支上，考虑速率和模式异质性)，所有这些似乎都是对实现准确的系统发育重建至关重要的。然而，到目前为止，还没有这种综合的实现方法/软件。造成这一空缺的一个主要原因是这些多种变异源所固有的计算复杂度，这将在联合实现中变的更加复杂，并因系统发育组学中的数据集大小而进一步加剧。

结论与展望

在多个方面，贝叶斯推断已经彻底改变了我们在系统发育学中的实践。理论上，贝叶斯推断为引入先验信息提供了一个灵活的框架。然而，在实践中，这并不是贝叶斯推断在进化遗传学中取得成功和流行的主要原因。相反，是分层模型和通用蒙特卡罗方法处理复杂随机效应和多层次进化过程的结合发挥了最重要的作用。

为跨位点的替代模式异质性建模代表了一个特定的例子，其中复杂的随机效应模型对实际的系统发育组学有重要的影响。这个问题在两个方面具有挑战性：第一，因为要建模的随机效应 (氨基酸偏好) 是高维的；其次，因为这些随机效应在各个位点上的分布本身是未知的，而且显然是复杂的。这种组合使得使用蒙特卡洛的贝叶斯推理特别适合，而更简单的方法，如有参随机效应模型或有限混合模型，现在已经显示出在计算成本上的高消耗，而且也不那么准确。

然而，贝叶斯推断也受到几个限制。首先，目前的蒙特卡罗算法不能随着数据大小的变化而很好地扩展，因此，贝叶斯推断可能不能适应大型系统发育数据集，对于基于狄利克雷过程先验的非参数模型尤其如此。其次，贝叶斯推断为建立新的、复杂的多层次模型所提供的灵活性在理论上是很好的，但在实践中，它需要人们为可能想要考虑的每个新模型进行大量的编程工作。这也提出了软件实现的可靠性问题，因为通常很难保证蒙特卡罗算法的给定实现确实是从预期的目标分布中进行采样的。

在这个方向上，当前的趋势是泛型编程平台 (generic programming platforms) 的开发。模型组件可以像积木一样组成复杂的层次结构，使用模块化的蒙特卡罗方法来探索所得到的后验分布，这使得贝叶斯推理的泛型编程相对简单，至少在概念上是这样。这种泛型编程平台在一般应用统计学和最近的系统发育学中都被提及，代表了一个很有前途的方向，通过为应用进化科学界提供面向用户的工具，以便可靠地设计问题特异的集成模型；并将其应用于经验数据的任意组合 (遗传序列、形态学数据、时间序列等)。然而，计算上的挑战是巨大的，在蒙特卡罗算法和软件开发中还有许多工作可做。

全文翻译自：Nicolas Lartillot. The Bayesian Approach to Molecular Phylogeny. Scornavacca, Celine; Delsuc, Frédéric; Galtier, Nicolas. Phylogenetics in the Genomic Era, No commercial publisher | Authors open access book, pp.1.4:1–1.4:17, 2020.