多物种蛋白编码基因的适应性进化检测方法

引言

20世纪下半叶，当来自不同物种的蛋白质编码DNA序列可以实现比对时，进化生物学家就开始试图比较，不改变编码氨基酸的替代 (同义替代) 和那些导致氨基酸变化的替代 (非同义替代) 的替代率了。早期的方法是基于简单的计数方案 (counting schemes)。他们的目标之一是解释并非所有密码子状态在同义和非同义的替代方面都具有相同的潜力。例如，编码色氨酸的密码子没有发生同义替代的机会，因为它只编码该氨基酸，而亮氨酸可由6个编码子编码，因此有很高的同义替代潜力。这些早期的方法依赖于多序列比对结果中所有或大多数两两序列的比对。然而，到了20世纪90年代，在完整的系统发育框架内开始出现一些统计模型。这些模型没有采用成对比较的计数方案，而是基于拟合参数 (例如，通过最大似然) 的想法，以表示运行在系统发育树分枝上相关序列的密码子替代过程的特征。

开发检测蛋白质编码基因适应性的密码子替代模型主要有两个思路。第一类包括基于参数的模型，对应于非同义 () 和同义 () 替代的速率比 ()。对的传统解释是，模型拟合的序列实例得到，即 (接近)中性演化，而表示纯化或负选择，而则对应于适应性演化或正选择。下面将讨论这种模型的详细形式，以及常用的扩展，允许不同密码子位置取不同的值。第二个建模思想集中于更好地捕捉净化选择的微妙之处，在后来被称为的突变-选择框架 (mutation-selection framework)。这些方法试图解释不同位点的氨基酸适应度分布异质性的问题，通过形成一个零模型 (null model)，以测试在突变-选择平衡下大于预期的非同义替代率。我们也将详细描述了突变-选择框架的基本原理，并用一个模拟研究来进一步描述该方法，还将其在数千个真实数据集上的结果与经典密码子替换模型进行对比。最后，我们讨论了一些违反密码子替代模型假设的事实，并概述了未来值得广泛关注的研究方向。

密码子替代的经典模型

在1994年9月的MBE杂志上，发表了Muse and Gaut (MG) 和 Goldman and Yang (GY) 的背靠背研究。他们将基于似然的系统发育分析中有关核苷酸状态空间，扩展到了三联密码子的状态空间，用的密码子替代率矩阵代替的核苷酸替代率矩阵。终止密码子在状态空间中是不允许的 (所以矩阵是维，而不是维的)。另一个内建的假设是密码子间的替代遵循点突变的过程，即相差两个或三个核苷酸的密码子之间的替代率被设置为零。然而，后一种假设是单核苷酸替代模型中所固有的，密码子的替代模型只是简单继承了这个假设。因为在给定的时间间隔内，随着区间时间趋近于零，两个核苷酸位点经历核苷酸替换的概率也趋近于零。

MG型-模型

得益于点突变假设，我们可以将核苷酸水平的模型框架，如GTR模型，直接应用于三联密码子的状态空间。令 表示一组对称的核苷酸替代速率参数 (nucleotide relative exchangeability parameters)，并有；令，表示一组核苷酸平衡频率参数 (nucleotide equilibrium frequency parameters)，并有。在GTR模型中，替代速率矩阵中的数据项由构成，只是对于核苷酸替代来说该速率矩阵是维的。同样的模型可以重写为一个的替代速率矩阵，用于描述从密码子到密码子的替代速率：

其中表示密码子的第位 ()，所以分别指代4种核苷酸。等式1和GTR模型之间的区别只是索引方式不同，但等式1中给出的公式表明，我们可以进一步识别不同类型的密码子替代。例如，在不考虑终止密码子的情况下，将其简化为的速率矩阵，就可以识别同义和非同义替代速率之间的区别，并将矩阵中的数据项指定为：

Muse and Gaut 1994 原文中，作者用和分别表示同义替代速率和非同义替代速率。和就是这里所谓的“乘法参数”。

等式2中给出的模型与Muse and Gaut最初提出的模型非常相似。最显著的区别是，在他们的开创性论文中，Muse and Gaut为非同义替代速率和同义替代速率调用了一个新的乘法参数，而等式2中，我们有效地将同义替代速率的乘法参数设置为，并在非同义替代速率上调用单个参数，因此就有了。与原始模型的另一个区别是，等式2包含了控制核苷酸可交换性的参数 (如GTR模型)，而Muse and Gaut最初只使用目标核苷酸的频率参数 (如F81模型)。密码子的平稳概率，可以看作是长时间运行替代过程后，出现密码子状态的比例，即：

分母是所有61个密码子状态的频率之和。换句话说，等式2给出的模型的平稳性几乎与在三个核苷酸位置上的GTR模型几乎相同，但在状态空间中没有停止码子的情况下，有轻微的重新归一化。

GTR模型的替代速率矩阵：

F81模型的替代速率矩阵：

对这个公式广泛使用的扩展之一，通常称为F3x4，是调用三个不同的核苷酸频率向量, 和，用于表示密码子内三个不同的位置 (Yang, 2006)。这样可得速率矩阵

平稳概率：

因为它涉及三组三维频率向量，F3x4背后的想法是解释由基因编码的三个位置的结构；例如，如果数据中有高度倾斜的组分偏差 (compositional bias)，它最能反映在第三密码子位置中，其状态通常对编码的氨基酸没有影响，而第一和第二密码子位置，其状态通常决定氨基酸，会有不同的核苷酸频率。换句话说，这三个位置的核苷酸频率的差异是在氨基酸水平上施加的不同选择压力的特征。虽然在核苷酸水平上用一个有丰富参数的模型来解释在氨基酸水平上选择的特征似乎有些不那么自然，理由是现象的：无论在氨基酸水平上的机制细节是怎样的，这种建模方法都在试图捕捉这些机制的净效应，至少是部分的。

Pond et al. (2010) 提出了一种修正方案，将这些参数的值设置为在从真实数据的三个密码子位置观察到的核苷酸频率。这些参数也可以用最大似然值来估计，或成为贝叶斯推断的一部分。

GY型-模型

由Goldman and Yang发展的模型与上述MG型-模型在几方面有所不同。其中最显著的区别是，GY型-模型的替代率矩阵中的项与目标密码子的频率 (或平稳概率) 成正比：

这与等式2中给出的MG型-模型形成了对比，MG型-模型的项与目标核苷酸的频率成正比。

大多数研究者根据核苷酸水平的规格来设置61维向量的值。也就是通常使用一个核苷酸频率的单一向量，表示为F1x4，并设置。然而，这种做法导致了密码子替代率规范过程中的特殊影响，且有时可能与建模意图相矛盾。Rodrigue et al. (2008) 指出了一个例子：假定一个更容易突变成或的演化语境，换句话说在这个演化语境中和要大于和；在所有其他条件不变的情况下，从到的替代速率 (即朝向高频状态T)，将低于从到的速率 (即朝向低频状态G)。同样，与核苷酸水平上的组分偏差 (到最终状态G) 的比率会更高，仅仅因为演化语境不同。目前还不清楚研究者是否完全意识到GY型F1x4模型的这些特性，或者他们是否认为它们对感兴趣的推断可以忽略不计，通常是 (或其分布)。

当然，F3x4的设置，在GY型-模型中也被广泛使用了。然而， Huelsenbeck and Dyer (2004)，以及Rodrigue et al. (2008)，已经表明这种的近似通常远远超出了61维向量全贝叶斯推断的95%置信区间 (通常表示为F61)。另一方面，GY-F61模型的主要缺点是混淆了核苷酸、氨基酸和密码子偏好。正如下面所要讨论的，MG型模型提供了在突变-选择框架内单独解释这些偏好的机会。

的概率分布

除了与MG和GY，或F1x4、F3x4和F61相关的问题之外，密码子替代模型的主要目标是描述非同义替代率与同义替代率的比率。如上所述，人们对的情况特别感兴趣。在前一节中，我们介绍了经典密码子替代模型的“均质”版本；对齐的每个密码子列都被认为是一个完全相同的马尔可夫过程的实现，在系统发育分枝或各对齐位置上保持不变。当将这些全局模型匹配到真实数据中时，我们几乎无法遇到的情况，因为适应性进化不太可能在所有密码子状态，位点和分枝中都起作用。

跨位点可变的

第一个解释整个密码子对齐中值变化的模型的灵感来自于Yang (1993; 1994) 的随机效应方法，这种方法在核苷酸水平模型中很好的模拟了不同位点的速率。他们认为每个密码子的对齐列都是由一个符合参数统计规则的模型 (以为参数) 产生的。然而，同样的问题仍然存在，即没有内在的理由来指导选择哪一个统计规则 (即概率分布)。在Yang、Nielsen和合作者的开创性工作中，他们所采取的方法是经验式的：探索许多不同的统计规则，并执行基于似然的模型比较，以确定最合适的规则。

捕获未知分布的最简单策略是将其离散为一个有限的混合模型。假设在密码子对齐位点上，我们允许个不同的参数，。第个密码子对齐列，在参数为的给定模型下的概率，可以表示为混合模型个分量上的加权平均值：

其中，且，是一组与每个分量有关的权重。这些权重可以看作是由每个混合模型的分量生成一个特定的对齐密码子的先验概率。

在他们所谓的中性模型 (neutral model) 中，Nielsen and Yang (1998) 设定 ，即存在两个不同的，分别，并通过最大似然法来推断出剩余的参数。 换句话说，他们的中性模型假设有两类密码子替代模型混合，一个是非同义替换率与同义替换率相匹配或相等 ，另一个是不允许发生非同义替代事件的 。他们的正选择模型增加了第三个分量 ，所以。我们可以在这两个模型之间进行似然比检验，以确定是否正选择的证据。这样的检验是在最大似然框架下检验适应性进化的一般方法的一个例子，通常接着是经验贝叶斯方法来识别高概率有的位点。

捕获夸位点值的未知分布也可以使用连续的参数分布。这与利用两个分量的中性模型分别对不同的位点进行建模的思想不同，我们也可以用在之间的连续分布，如beta分布。在此模型下，位点的似然会有以下形式：

其中为beta分布 (以和为参数，也是需要进行最大似然估计的参数) 下的概率密度，相当于等式7中的。在实践中，等8中的积分常通过离散化技术进行近似，将其简化为一个很像等式7中的加权和，但总体方法仍然允许一个紧凑的参数化来解释各个位点的异质性。

当然，我们也可以设想由连续分布和离散类别的混合而构建的模型，以允许具有正选择的位点，从而导致具有以下形式的位点似然：

与前面一样，后一种模型，连同前一个单独基于beta分布的模型，可以形成具有适应性进化特征的位点存在的似然比检验的基础。

基于离散和连续分布的混合模型，或几个连续分布的模型，可以以类似的方式建立，这导致了Yang et al. (2000) 提出了一套类模型，并对适应性替代体系进行了相关的似然比检验。实际上，后一种检验基于附加了组分的beta分布模型与只有beta分布的模型进行比较，被称为M8对M7的检验 (the M8 versus M7 test)。

对的变化进行越来越微妙的建模

从历史上看，密码子模型的发展主要以上述方式发展：重点是提出模型，允许适应性进化越来越微妙的表现，如在特定位点和/或特定分枝上运行的适应性进化机制，通常通过使用各种统计方法控制位点和分枝的值。这一焦点可能已经将注意力从更丰富的突变特征建模中转移开，同时也阻碍了对使用作为检测适应性进化的适当方法的更普遍的质疑。

事实上，对值的解释一直是争论的焦点。Nielsen and Yang (2003) 在这个参数与基本的种群遗传学理论之间建立了间接联系。它们将的值与尺度选择系数 (scaled selection coefficient, ) 联系起来，通过表达式来量化与特定氨基酸替代相关的适合度 (fitness) 的变化。这种关系导致了一些奇怪的情况。例如，意味着所有非同义替代（在给定的位点和/或分枝）都有；某氨基酸位点从“L”变为“I”的有正选择系数，反过来从“I”到“L”的替代也有。相反，暗示每一个非同义替代都会降低适应度 ()，包括“D”到“E”和“E”到“D”。

与此同时，另一种建模思想成了一个完全不同的焦点：设计一个在蛋白质编码DNA序列上普遍的潜在纯化选择的更好的表示方法。

突变-选择框架

1998年，Halpern and Bruno ("HB") 利用植根于种群遗传学的概念，发展了一个可以更直接解释适应性演化的新模型。他们的基本思想是明确地识别密码子替代的初始状态和最终状态，而不是简单地区分同义和非同义替代事件，以及核苷酸或密码子的最终状态。Yang and Nielsen (2008) 明确地介绍了该模型的概念，并引入了密码子的适合度 (fitness) 参数。从野生型状态到突变状态的变化就意味着一个选择系数。突变的固定概率为 (近似)：，其中是有效的染色体群体大小 (effective chromosomal population size)。在涉及单倍体的背景下，直接对应于有效种群大小，而对于二倍体，我们的符号意味着是有效种群大小的两倍；换句话说，我们在这里提到的包括一个依赖于倍性的乘法因子。突变过程可以被指定为，例如核苷酸水平的GTR模型，其中密码子到的突变率由给出，因此染色体群体水平的突变率为。这两个概念用乘法组合起来，可以定义替代率：

通过将最左边的乘以固定概率 (见以上近似公式)，用替换，可定义为缩放的选择系数 (scaled selection coefficient) (以有效染色体群体大小的两倍缩放)，等式10给出的模型可以重写为：

其中，而为密码子缩放的适合度。我们可以通过将其中一个参数固定在上，并用这个约束来估计剩余的适合度参数；毕竟我们关心的主要是相对的缩放适合度。然而，另一种选择是设置，其中，且，是一个密码子频谱 (codon profile)。在突变-选择框架，、和的解释分别对应负选择、中性选择和正选择系数。

profile, 直译为侧面，半面；外形，轮廓；[航]翼型；人物简介。多翻译为“侧写”。其意义为通过一种方式描述事物的特征。

这里根据的数学定义，它所描述的是某个位点上密码子的比例/频率。所以这里我们翻译为“频谱”。

该模型下密码子的平衡概率为

或

在等式12和13中，我们定义了控制核苷酸偏好的参数和控制密码子适合度的参数。平衡概率是基于和计算的，但模型的呈现表明，需要被估计的是和 (或)。然而，在Halpern and Bruno (1998) 提出的模型中，被估计的是和，其中 (或) 是隐含的。具体地说，他们构造的替代矩阵为：

将等式13引入等式14会产生等式11，可以很好的说明这些方法的等价性。

HB型-模型

突变-选择框架的形式是以全局式模型的方式呈现的。然而，Halpern and Bruno (HB) 模型的一个核心思想是为每个位点设定一组独特的密码子频谱。此外，在实践中，他们将模型简化为只包含氨基酸频谱 (amino acid profiles)；这将为每个位点生成一个模型：

其中为在位点上的氨基酸频谱，函数返回的是由密码子编码的氨基酸的指数。注意，我们可以将等式15重写成类似等式11的形式：

因此，序列比对中有多少个密码子列，HB型-模型就会指定多少个密码子替代矩阵，虽然每列都共享一个核苷酸水平的参数化，但每个列都会有一组独特的氨基酸适合度参数。

该方法的应用受到其高维数的阻碍，直到十年后才再次有人使用这些位点特异的参数 (Holder et al., 2008; Tamuri et al., 2012, 2014)。这里还有另一个担忧，即将位点特异的氨基酸频谱作为用于ML估计的真正 (bona fide) 参数，这并不符合似然法推断的渐近理论的条件。位点特异的方法没有渐近条件，因为将它们应用于越来越多位点的数据集意味着引入新的氨基酸频谱，从而会改变模型；将它们应用于具有更多序列的数据集，也会因为需要在不同树上添加额外的分枝长度，而改变模型。另一种建模策略，通常应用于大多数系统发育分析，是随机变量的方法。

突变-选择模型的随机变量方法

正如前面聚焦于的经典密码子模型中所描述的，随机变量方法也已经应用于氨基酸适合度谱 (amino acide fitness profiles)，包括参数的 (Rodrigue, 2013) 和非参数的 (Rodrigue et al., 2010b; Rodrigue and Lartillot, 2014) 方法。与以前一样，氨基酸适合度谱被视为随机变量，可以用统计学规则描述。以前的研究已经使用过的统计规则，包括了20个氨基酸状态的纯平狄利克雷分布 (plain flat Dirichlet)；一个本身拥有控制它的中心和集中度的参数的自由狄利克雷分布 (free Dirichlet)；具有经验值 (empirically derived values) 的有限混合模型。

非参数的版本，利用了氨基酸适合度谱的狄利克雷过程，并通过“中餐厅 (Chinese restaurant)” (Rodrigue et al., 2010) 和“折棍子 (stick-breaking)” (Rodrigue and Lartillot, 2014) 两种构造过程实现。这两种方式都使用一个辅助变量，表示为密码子位点所分配的“活性”氨基酸适合度谱 (个集合之一)，位点的替代模型通常表示为：

其中为，基于分配给位点的氨基酸频谱 () 的，缩放的选择系数。等式17中给出的模型有时表示为MutSelDP。

尽管被称为是非参数的，Dirichlet过程确实涉及到参数 (通常被称为超参数)，指定一个类似于均值的基分布和一个粒度参数，控制着未知分布估计的粗糙度。似然函数可以被认为是无穷组混合模型上的积分，基于这些超参数；该积分可用蒙特卡罗方法有效地逼近。到目前为止，在研究Dirichlet过程不同超参数方面的工作相对较少；之前的研究赋予了他们自己的简单统计规则 (超先验)，并将它们视为推断的自由元素。更丰富的超先验，如由混合的Dirichlet过程作为基分布本身，应该在未来的工作中进行研究。

突变-选择框架作为检验适应性的零模型

Halpern and Bruno提出的突变-选择框架背后的基本动机是定义一个更好的零模型作为理解蛋白质编码基因进化特征的起点。具体地说，该框架专注于以一种位点异构的方式捕获净化选择。Spielman and Wilke (2015) 清楚地阐述了这些模型如何在平衡状态下引起小于1具有影响。我们将由突变-选择模型产生的称为，通过检查两种极端情况，可以直观地看到为什么它被限制在0到1范围内。首先，假设由于某种原因，所有的氨基酸都几乎具有相同的适合度；然后，所有的值都将接近0，进而，这样，当考虑到突变机会时，非同义率和同义率非常匹配 (即接近1)。另一个极端情况是，假设氨基酸适合度谱由单一氨基酸强烈支配；然后，在平衡状态下，处于这种主要氨基酸之外的状态很少，所有可能的突变导致非常负，因此，从而导致非同义替代速率接近0 (即接近0)。在这两种极端情况下，氨基酸谱的其他配置会导致处在0到1的范围内。

这些想法提出了一种检验适应性演化的替代方法：与检测哪里发生了的情况不同，我们可以瞄准检测整个大于在纯突变-选择模型下预期值的情况。一种方法是对非同义替代事件引入一个乘法参数，我们表示为 (星号用来区分经典密码子模型中的)，导致以下Rodrigue and Lartillot (2017) 的模型：

等式18中给出的模型称为MutSelDP-。有了这个公式，平衡状态下模型的整体可以被认为是。或者，写成，这个新参数可以看作是整体对纯突变-选择模型预期的偏差的度量。表明，在突变-选择平衡下，整体非同义率大于预期，即存在一个潜在的适应性演化。这种方法比经典的密码子模型要求低得多，即要求非同义速率实际上超过同义速率，因此有潜力检测到其他可能被忽视的适应表现。更根本的是，该方法展示了一个不同的建模视角：制定一个更好的零框架，以发现与这个新零模型的更微妙的偏差。

模拟研究

（略过）

基于真实数据的突变-选择框架与经典模型的比较

（略过）

结论

虽然MutSelDP-是迄今为止提出的最丰富的密码子替代模型之一，但它为检测适应性演化所做的参数化是在非同义替代率上施加了一个简单单变量乘子 ()。然而，它已经能够在模拟中检测到适应性演化的微妙实例，并有潜力检测到真实数据中类似的、与具有非同义替代率乘子分布的经典模型相似的、适应性演化信号。

当然，扩展MutSelDP-模型，以允许在位点和/或分枝之间进行值的分布，将是非常有趣的。这与过去几十年用经典密码子模型探索的精神一样。这种具有多种独立类型异质性的模型 (如氨基酸谱的分布和的分布) 提出了重大的挑战，尤其是它们的计算负担。一些捷径，如利用经验的氨基酸谱的混合，以及预设的值网格 (preset grid of )，对于快速的对大数据集进行第一通分析，是值得考虑的。

探索更多的模拟，将更多进化过程的已知特征纳入源自数据的模型也是有吸引力的。真实数据时常会违反等式10的模型假设，即时间同质的有效群体大小 () 的假设。了解MutSelDP-模型，或其最终扩展，如何对树上有效群体大小变化的模拟数据做出反应，将是重要的第一步。最近的其他模拟表明，CpG超突变性会误导一些密码子替代模型在同义密码子使用方面检测选择。更普遍地，应用模拟来评估这些违反模型的影响，对于基于突变-选择模型是为了更仔细地校准它们所导致推论的可靠性。

最后，长期建模目标可能应该寻求将最近的创新整合到一个单一模型中，该模型可以适应密码子使用不均匀、系统发育的可变有效群体大小、环境依赖的突变率和组织效应 (基因内和基因之间)。一些初步的工作阐述了开展该项目的技术手段，包括近似贝叶斯计算和nested-MCMC系统。将这些想法结合在一起可以完成一个构建模型的框架，这些模型越来越忠实于现代生物学对分子进化的理解。

全文翻译自：Christine Lowe, Nicolas Rodrigue. Detecting Adaptation from Multi-species Protein-coding DNA Sequence Alignments. Scornavacca, Celine; Delsuc, Frédéric; Galtier, Nicolas. Phylogenetics in the Genomic Era, No commercial publisher | Authors open access book, pp.4.5:1–4.5:18, 2020.