摘要：

        PAML软件包中的CODEML程序已被广泛用于分析蛋白质编码基因序列，以估计同义和非同义率(dS和dN)，并检测驱动蛋白质进化的达尔文正选择。对于不熟悉分子进化分析的用户，该程序有一个艰难的学习过程。在这里，我们提供了一个循序渐进的方案来说明程序中可用的常用测试，包括分支模型，位点模型和分支-位点模型，它们可以分别用于检测驱动适应性蛋白质进化的正选择，影响物种系统发育的特定谱系，影响蛋白质中的氨基酸残基子集（趋同进化），以及影响预先指定谱系的子集。本文以10种哺乳动物和2种鸟类黏液病毒(Mx)基因的数据集为例。我们讨论了CODEML中的一个新特性，它允许用户对同一组分类群的多个基因执行正选择测试，这在现代基因组测序项目中很常见。在GNU许可下，PAML包在https://github. com/abacus-gene/paml上发布。https://groups. google.com/g/pamlsoftware讨论站点提供支持。此协议中使用的数据文件可从https://github.com/abacus-gene/ paml-tutorial获得。

前言

        分子进化的中性理论认为，大多数观察到的物种内部和物种之间的变异不是由于自然选择，而是由于具有小适应度意义的突变的随机固定(Kimura 1968;King and Jukes 1969)。换句话说，有利的突变在分子水平上是罕见的。然而，有利的基因和基因组突变最终决定了物种的形态、行为和生理，或者导致了物种的分化和进化创新。因此，检测分子适应使我们能够更好地理解进化过程。鉴于大量基因组数据和计算资源的可用性，研究分子适应目前比以往任何时候都更令人兴奋，因为现在有可能系统地整合基因组数据，以在广泛的生物体中寻找正选择的信号。

        在蛋白质编码基因中，我们可以区分同义或沉默取代(不改变编码氨基酸的核苷酸取代)与非同义或替代取代(改变氨基酸的核苷酸取代)。由于自然选择主要在蛋白质水平上起作用，同义突变和非同义突变在非常不同的选择压力下以非常不同的速率固定。因此，以同义率作为参考点，人们可以判断作用于蛋白质的自然选择是否加速或减缓了群体中非同义突变的固定。因此，同义和非同义取代率的比较可以揭示作用于蛋白质的自然选择的方向和强度(Kimura 1977;Miyata and Yasunaga 1980)。非同义取代率加快的基因，其非同义/同义的比率dN/dS >1,表示是在经历正选择。这种测试在检测多样化选择或平衡选择方面特别有效，因为它使用过多的非同义替换作为自然选择帮助固定非同义突变的证据。基于dN/dS的测试在应用于来自同一物种的数据时可能效果较差，因为缺乏序列差异，并且由于dN/dS比率解释的复杂性(Kryazhimskiy和Plotkin 2008)。

        有了来自不同物种的蛋白质编码基因序列，人们可能会提出以下问题。是否有密码子在dN/dS>1的正选择下进化?如果有，我们如何识别它们?在某些谱系中，是否存在正选择推动着密码子的快速替换，例如，在基因复制之后或当一个物种适应新环境时?什么样的统计检验适合回答这些问题?如何估计选择的强度?

        许多马尔可夫链模型已经被开发出来，用于检测影响基因、氨基酸残基或进化谱系的正选择(Kosakovsky Pond等人，2005;杨2007)。在本使用指南中，我们讨论了密码子模型在序列比对的最大似然系统发育分析中的使用，以检测驱动优势非同义突变固定的正选择。

密码子取代模型下的适应性分子进化检测
       几个密码子替换模型已经在软件CODEML中实现，它是PAML包的一部分(Yang 2007)，作为Goldman和Yang(1994)以及Muse和Gaut(1994)模型的扩展。与核苷酸或氨基酸取代模型相比，密码子取代模型的主要特点是，密码子三联体被认为是进化的单位(Goldman and Yang 1994)。常用的模型版本(如Yang 1998;Yang和Nielsen 1998)忽略了氨基酸之间的化学差异，并使用相同的非同义/同义比率(即ω = dN/dS)，该比率不依赖于密码子编码的源和目标氨基酸。这一假设大大简化了对模型的解释。具体来说，ω比值测量了氨基酸选择变化的方向和幅度:ω 的值:<1， = 1，>1，分别表示为负选择(净化选择)、中性进化和正选择。然而，在一个基因的所有位点和系统发育上的所有谱系(分支)上的平均ω比率很少大于1，因此，它用于检测正选择实际上是没有力量的。相反，检测影响特定分支或单个位点的正选择被证明更有用。

        分支模型对不同分支的系统发育假定不同的ω比值参数(Yang 1998;Yang和Nielsen, 1998)。它们可以用来检测作用于特定谱系的正选择，而不需要在整个系统发育树中平均ω比率。例如，它们对于检测基因复制后的阳性选择是有用的，在这种情况下，复制的一个拷贝可能获得了一种新功能，因此可能以更快的速度进化。

        位点模型将基因中任何位点(密码子)的ω比率视为来自统计分布的随机变量，因此允许ω在密码子之间变化(Nielsen and Yang 1998;Yang et al. 2000)。正选择定义为存在某些密码子的ω >1. 执行似然比检验(LRT)与null模型的一般密码子比较，不允许任何密码子ω>1。有两对位点模型特别有效，它们构成了使用CODEML进行正选择的两种广泛使用的测试(表1):(1)M1a(接近中性)vs. M2a(正面选择)(Nielsen and Yang 1998;Wong et al. 2004;Yang et al. 2005)和(2)M7 (beta) vs. M8 (beta&ω) (Yang et al. 2000)。LRT统计量，或两倍于两个比较模型之间的对数似然差(2Δ)，可用于卡方检验，其自由度为两个模型之间自由参数数量的差异。例如，M1a有2个自由参数，M2a有4个，因此自由度为2，卡方检验的临界值为χ²_2,5%显著性水平下= 5.99,1%显著性水平下χ²_2,1% = 9.21。M7-M8的比较也有2个自由度(表1)。

        如果LRT倾向于M2a或M8，我们可能会问:基因中的哪些位点处于正选择之下?利用贝叶斯定理计算每个位置来自ω>1的正选择类的后验概率来回答这个问题;1. 在CODEML中实现了两种方法。Naïve Empirical Bayes (NEB)方法(Nielsen and Yang 1998)通过使用模型中参数的最大似然估计(MLEs)计算每个位点来自不同位点类别的后验概率，而不考虑其抽样误差或不确定性。贝叶斯经验贝叶斯(BEB)方法是对NEB方法的改进，并适应了MLEs中的不确定性(Yang et al. 2005)。虽然CODEML报告两种方法的结果，但应该使用BEB而不是NEB。

        分支位点模型旨在检测只影响预先指定谱系中少数位点的正向选择(Yang和Nielsen 2002)。接受正选择测试的分支称为前景分支，而树上的所有其他分支都是背景分支。对于背景分支，有两类站点:0 &lt的保守站点;ω0 & lt;和ω1 = 1的中性点。For the back ground branches, there are two classes of sites, the con served sites with 0 < ω0 < 1 and the neutral sites with ω1 = 1. For the foreground branches, some of those sites become under positive selection with ω2 > 1. Positive selection or the presence of sites with ω2 > 1 is tested by comparing this model with a null model in which ω2 = 1 is fixed, using a 50:50 mixture of 0 and χ²₁ as the null distribution (Yang et al. 2000 ; Zhang et al. 2005 )。与位点模型一样，BEB方法可以用于识别前景支上可能处于正选择的密码子位点。

        这种常用的分支站点模型是Yang和Nielsen(2002)的老分支站点模型A和B的改良版。旧模型A将保守位点的ω0 = 0和中性位点的ω1 = 1固定下来，这是不现实的，因为它不能适应0 &lt约束下的位点;ω& lt;1.   旧模型B同时使用ω0和ω1作为自由参数，当估计的ω1略大于1时，当弱约束下的站点被分配给这样一类正选择时，可能会有高假阳性。修改后的分支站点模型A通过使用0 <ω0 <1作为保守位点的自由参数，并且通过固定ω1 = 1来解释接近中性或弱约束的位点(Yang et al. 2000,2005;Zhang et al. 2005)。

        重要的是要注意前景支应该预先指定。如果在没有先验生物学假设的情况下使用相同的数据集对树的多个分支进行正选择测试，则可能需要对多个测试进行校正(Anisimova和Yang 2007)。Bonferroni修正可能过于保守，而Rom程序(Rom 1990)的功效略高，是首选(Anisimova and Yang 2007)。人们还可以使用控制假阳性(FDR)的程序，它是所有被拒绝的零假设中真零的预期比例，或者是所有阳性测试结果中假阳性结果的比例(Benjamini和Hochberg 1995,2000)。需要注意的是，如果这些序列非常分散或存在严重的模型违规，则多重测试校正可能是不可靠的(Anisimova and Yang 2007)。

        在指南中，我们提供了在上述模型下执行分析的详细步骤。我们特别关注使用CODEML进行四种类型的分析:(1)在所有位点和分支均为1个 ω的同质模型下，计算ω作为对基因的平均选择压力的度量，(2)检测影响编码序列中在正选择下进化的位点子集的正选择(位点模型)，(3)检测在正选择下进化的系统发育的一个或多个分支(分支模型)，以及(4)检测特定分支的位点子集(分支-位点模型)。我们使用哺乳动物和鸟类黏液病毒序列比对来说明分析结果(Hou et al. 2007)。最后，我们讨论了CODEML的最新实现，它可以使用由数千个基因比对组成的基因组数据集进行正选择测试。

操作指南

1、获取和编译PAML套件程序
        指南的所有分析都使用最新版本的CODEML进行，这是PAML包的一部分(撰写本文时为PAML v4.10.6)，可以从PAML GitHub存储库下载（https://github. com/abacus-gene/paml)。在不同的操作系统上安装PAML的说明可以在http:// abacus.gene.ucl.ac.uk/software/找到。

        在整个指南中，我们展示了如何在UNIX(例如Mac OSX)和Linux环境(例如Ubuntu, Centos等)中从命令行执行CODEML。Windows用户可以从Windows命令提示符运行CODEML，也可以从Windows 10或更高版本的Windows Subsystem for Linux (WSL)中启用的Ubuntu终端运行CODEML。以下连结提供使用终端机的简短教程:

http://abacus.gene.ucl.ac.uk/software/CommandLine. Windows.pdf

http://abacus.gene.ucl.ac.uk/software/CommandLine. MACosx.pdf

        我们假设可执行文件的路径已经正确地导出到用户的系统中(即，添加到path全局变量中，该变量通常可以通过~ /.Bashrc或~ /.bash_profile隐藏文件)，因此只要在终端上输入可执行文件的名称(即CODEML)就可以执行CODEML。如果不是这种情况，则应该提供可执行文件的完整路径。

2、运行CODEML

        对于本指南中提出的所有分析，我们使用了一个由10种哺乳动物和2种鸟类黏液病毒(Mx)基因的12个序列组成的数据集，该数据集由Hou等人(2007)生成并分析(见图1)。Mx基因参与宿主物种的抗病毒反应。侯等人(2007)假设该基因在鸡和鸭谱系中是在正选择下进化的，而在哺乳动物谱系中是在强净化选择下进化的。

        我们已经创建了一个GitHub存储库(https://github. com/abacus-gene/paml-tutorial/tree/main/positive-selection)作为补充材料，我们提供了一个循序渐进的教程，详细介绍了如何下载和过滤黏液病毒序列数据，然后生成我们在整个协议中使用的基因比对。

3、控制文件
        控制文件包含CODEML运行任何分析所需的信息。PAML文档中提供了控制文件中变量的详细描述。在这里，我们关注与本协议中详细分析相关的变量(见图2)。我们注意到一些关于控制文件的规则:

            1)控制文件逐行读取。如果同一变量被指定两次，第二个变量将覆盖第一个变量。

            2)控制文件中的一些规范对所有替代模型都是通用的。

            3)空白行被程序忽略。同一行中星号(*)之后的任何内容都被视为注释。、

4、输入/输出文件
        示例控制文件(.ctl)的第一个块(参见图2)指定了输入文件、序列比对文件(seqfile.phy)和系统发育树文件(treefile.newick)的路径，以及保存主要分析结果的输出文件(outfile)的路径。对齐文件应该是phy格式。标题(对齐之前的行)详细说明了分类群的数量和对齐的长度。例如，本协议中使用的Mx基因比对文件有12个物种和1989个碱基对，因此标题为“12 1989”。如果输入文件中包含多个基因对齐，则假定每个基因对齐从标题行之后开始。树形文件包含Newick格式的树形拓扑(如果存在，应该删除引导值)，并且可能包含一个头文件。例如，对于包含12个分类群的Newick树，应该在Newick树之前添加标题“12.1”。分支长度不是必需的，并且考虑到它们可能会干扰某些模型中的附加信息(例如，节点标签)，我们建议不将它们包含在Newick树中。如果在树文件中列出了多个树，则需要一个头文件。有关数据格式的讨论，请参阅补充材料中的“对齐，序列数据文件和树”一节和我们的GitHub存储库(https://github.com/abacus-gene/paml-tutorial/tree/main/ positive-selection/00_data)。

5、屏幕输出
        要在屏幕上或输出文件中打印的信息量由变量noisy和verbose控制。

6、输入序列数据
        要在推理中使用的数据类型(例如，核苷酸、氨基酸或密码子，它们将被CODEML翻译成氨基酸)使用变量seqtype指定。基因或排列的数量是用变量ndata定义的，而遗传密码是用变量icode指定的。此外，变量cleandata可用于决定歧义数据和比对间隙的位点是否保留(cleandata = 0)或删除(cleandata = 1)。在使用CODEML之前，序列必须正确对齐，内含子，非编码区和终止密码子被删除。为了保存阅读框，一种有用的策略是首先对蛋白质序列进行比对，然后相应地构建密码子比对。我们建议删除主要是间隙或难以对齐的对齐列或区域，然后使用cleandata = 0来保留数据中的信息。另请参阅补充资料中的“比对、序列数据文件和树文件”和“基因树与物种树”章节中的建议。

7、替换模型
        使用几个变量来指定分析中要使用的模型。首先，ω可能在谱系(模型)、跨位点(NSsites)或跨谱系和位点(模型和NSsites)之间变化。其次，密码子频率(CodonFreq)可以相等(每个1/61)或不同，以解释密码子使用偏差。FmutSel模型(CodonFreq = 7)为通用遗传密码的每个密码子分配了60(= 61-1)个密码子适应度参数(Yang and Nielsen 2008)。FmutSel0模型(CodonFreq = 6)是FmutSel的特例，它对同义密码子赋予相同的适应度值，因此只使用19(= 20-1)个氨基酸适应度参数。该模型假设氨基酸频率由蛋白质的功能需求决定，但同义密码子的相对频率仅由突变偏倚参数决定。在这些突变选择模型中，变量estFreq指定使用从数据中观察到的密码子频率或通过最大似然估计。根据用这些参数定义的模型，ω的值可以固定或通过变量fix_omega和omega估计。可变时钟用于加强分子时钟(即，谱系之间的速率常数)。由于我们在这里讨论的模型都是时间可逆的，并且我们没有假设分子钟(clock = 0)，因此应该使用无根树。请注意，几乎所有的系统发育程序，包括RAxML (Stamatakis 2014;Kozlov et al. 2019)， IQ-TREE (Minh et al. 2020)和MrBayes (Ronquist et al. 2012)生产了无根树。严格时钟和本地时钟模型在CODEML中通过使用clock = 1和clock = 2来实现，但是在本指南中没有使用这些模型。关于树是否生根的讨论，请参见补充材料中的“有根树与无根树”一节。

8、跨谱系和位点具有均匀ω的一个比值
        在CODEML中实现的最简单的密码子模型是M0 (one-ratio)，它假设所有位点和所有谱系的ω比为1(图3A)。在大多数情况下，这个假设是不现实的，因为大多数位点都受到ω <1的约束. 因此，试图在M0模型下识别正选择的证据将会失败。事实上，当我们感兴趣的是检测正选择时，最佳实践是在协议后面描述的位点或分支位点模型下估计ω。然而，M0模型为我们提供了一个零假设，可以作为比较的参考，以确定更复杂的模型是否能更好地拟合数据。M0下的ω估计也给出了在位点和谱系上平均的选择性约束的总体度量。

        使用控制文件CODEML - m0运行CODEML。图2的CTL，我们打开一个终端，导航到包含输入文件(序列比对文件、树文件和控制文件)的文件夹，然后输入: codeml codeml-M0.ctl

        屏幕输出应该是不言自明的，其中包括有关分析、对齐和ML迭代进度的信息。如果输入文件中有错误，或者程序由于其他原因而失败，屏幕上也可能会打印错误消息。分析完成后，还将在输出文本文件(out_m .txt)中打印有关数据和结果的一些信息。该输出文件分为五个部分，如下所述。

out_m .txt_9、输入对齐中的位点模式概要
        输入对齐及其压缩版本打印在输出文件的顶部，其中每个位点模式只表示一次，相应的频率(模式计数)在下面的块中显示。例如，在这里使用的基因比对文件(Mx_aln.phy)中，大多数位点模式是唯一的，频率为1，而两个位点模式重复了2次和4次，如下所示。如果使用cleandata = 1，则将打印删除间隙和模糊位置之前和之后的对齐，然后是修剪对齐的位点模式计数。请注意，如果删除有间隔的列，则位点将重新编号，这可能会影响位点或分支位点模型下的输出。

out_m .txt_10、输入对齐和所选模型的摘要
        接下来，用于分析的CODEML版本连同对齐文件的名称一起打印出来。然后是对照文件中指定的模型的详细信息，比对文件中物种数(ns = 12)和密码子数(ls = 663，碱基对数[1989]分为3):

out_m .txt_11、核苷酸和密码子频率综述
接下来，每个序列中观察到的核苷酸和密码子频率以及它们在所有序列中的平均值被打印在输出文件中，然后是模型下的密码子频率，可以使用evolver等软件用于模拟目的(Yang 2007)。

out_m .txt_12、树得分和估计模型参数值的总结
        对数似然值、自由参数个数、估计分支长度以及模型下的其他参数(本例中为M0)如下面的前四行所示。自由参数包括分支长度、平衡频率、替换/颠换速率比(κ)和ω分布中的参数。注意,无根树与n个分支：2×n−3分支长度(即2×12−3 = 21我们的树的分支长度n = 12个分类群)和3 mutation-bias参数(四个频率之和为1)。它可能很难匹配的ml 4号线到相应的分支系统学的长度或模型参数,但Newick树中的分支长度再次打印。最后一块显示了κ(替换/颠换率比)和ω的估计，以及突变偏差(核苷酸频率)参数。突变选择模型通过模拟突变偏差和选择下突变的固定来适应不同的密码子频率:较高的突变偏差参数，例如π * A，意味着突变过程偏向于A (Yang和Nielsen 2008，等)。1和4)。下面是每个分支的t(分支长度)、dN(非同义速率)和dS(同义速率)的估计值。输出如下所示:

        由于指定了同质模型(M0) (model = 0和NSsites = 0)，因此在dN/dS列下的每个分支报告相同的ω值(0.3357)(见上面的粗体值)。非同义位点的树长dN = 1.1386，同义位点的树长dS = 3.3921， ω = dN/dS = 0.3357，表明黏液病毒基因总体上处于纯化选择（ ω<1）阶段，非同义突变进入固定的几率(ω = dN/dS = 0.3357)平均为同义突变的三分之一。

        在接下来的部分中，我们将重点介绍用于检测正选择的三种类型的模型:位点模型、分支模型和分支位点模型。

13、跨位点的异质位点模型
        在本节中，CODEML在各种位点模型下运行，这允许ω在密码子之间变化(图3B)。我们可以在一次运行中分析几个位点模型下的相同数据。在这个例子中，我们在同质M0模型下运行CODEML(和以前一样)和以下位点异构模型:M1a, M2a, M7和M8(见表1)（NSsites=0 1 2 7 8 *每个数字间隔一个空格）。我们在图2的控制文件中编辑以下控制变量，并将其重命名为CODEML -sites.ctl:

        现在我们可以从命令行执行CODEML: CODEML -sites.ctl
        输出文件包含控制文件中定义的五个模型中的每个模型的一个部分，其中提供了对数似然值和总参数的数量，后者用于确定每个模型比较的自由度。

        为了比较嵌套的位点模型，我们使用LRT统计量，定义为零假设和备选假设之间的对数似然差的两倍，2Δ ρ = 2(ρ 1−ρ 0)，其中ρ 0是零模型的对数似然评分，ρ 1是备选模型下的对数似然评分。LRT统计量与χ2分布进行比较，其自由度等于两个模型之间自由参数(在每个模型的输出文件中给出)数量的差异。结果总结在表3中。

        模型M0(one-ration)和M1a(近中性)是嵌套的，可以使用LRT进行比较。这是对氨基酸位点之间的选择压力变异性的测试，而不是对正选择的测试。M1a拟合数据比M0要好得多，其2Δ = 559.26 >χ²_1,5%，说明ω反映的选择压力在不同位点间差异很大。与M1a相比，M2a增加了一类ω2 >1的正选择位点(按p2的比例)。这并没有显著改善模型的拟合，因为2Δ r = 0(表3)。我们通过比较M7 (β，零模型)和M8 (β&ω，替代模型)对正选择进行了额外的检验。在5%的显著性水平上，M8比M7更适合数据，2Δ χ²_2,5% = 5.99，表明存在ω >1的正选择位点. 因此，对正选择位点的测试是模棱两可的，M1a-M2a和M7-M8的比较给出了相互矛盾的结果。当正选择的证据存在但不是很强时，注意到M1a-M2a检验更为严格，因为弱正选择的位点往往被集中到ω1 = 1的位点类别中(Zhang et al. 2005)。

        各位点模型下各参数的最大似然值见表4。可以在输出文件的相同部分中找到参数估计，这些部分已经用同质模型进行了说明。您可以使用bash脚本从输出文件中提取这些值。或者，我们在GitHub存储库(https://github.com/ abacus-gene/paml-tutorial/tree/main/positive-selection/00_ data)的分步教程中包含代码片段。

      由于指定了同质模型(M0) (model = 0和NSsites = 0)，因此在dN/dS列下的每个分支报告相同的ω值(0.3357)(见上面的粗体值)。尽管M0假设相同的ω为所有基因的密码子(即模型= 0和NSsites = 0),该位点模型假设几个位点类(图3 a和B)。例如,在M8(β&ω),94.2%的位点(p0 = 0.942的估计价值,见下表4 M8)的ω值服从β(0.404,0.750)分布(贝塔分布的参数的估计价值,p = 0.404和q = 0.750,见下表4 M8),而5.8%的位点(p1的估计价值= 0.058,(见M8下表4)的ω = 1.841（即ω>1），表明存在少量正选择的氨基酸残基。这些信息可以在输出文件中以术语MLEs开头的行之后找到。当在一次运行中使用多个位点模型时，每个模型都有自己的输出信息块，其中报告了参数估计。

        在M2a和M8下，使用BEB方法计算来自不同位点类别的每个位点的后验概率。正选择类的后验概率高的位置可能处于正选择之下。M8下的输出如下:

        在每一行中，第一列显示位点位置(例如，10、25、108和123)，后面是第一个序列中该位点上的氨基酸(这是为了识别序列中的位点)。第三列(Pr (w >1))表示该位点来自正选择类的后验概率(即ω>1).最后一列表示该站点ω分布的后验均值和标准差。请注意，BEB计算仅在正选择模型(即M2a和M8)下进行，而不是在零模型(即M1a或M7)下进行。

        在我们的例子中，14个位点中具有50%的后验概率的ω >1，上面列出了其中的4个。例如，位点10在第一个序列中含有丝氨酸，来自正选择类的后验概率为0.564，该位点ω的后验分布均值为1.468,SD为0.634。这里，由于ω >1的后验概率较小（0.564），BEB计算仅为正选择位点提供了微弱的证据较低，ω的后验分布在每个部位呈弥漫性分布。此外，LRT在M1a-M2a比较中不显著，在M7-M8比较中仅在5%水平上显著。综上所述，这些结果为Mx基因中的正选择位点提供了一些证据，尽管这一证据不是很有力。请参阅补充材料中的“BEB分析”部分，以获得正选择位点的示例。

14、具有跨分支异构ω的分支模型
        在本节中，我们将展示如何在分支模型下运行CODEML(图3C)，假设ω在树的各个分支上变化。分支模型通过在树文件中使用标记来指定:#0(默认)、#1等。该模型目前允许最多8种不同ω比的支路类型。这里我们只考虑两种类型。前景分支被标记为标签#1(并分配比例ω1)，而所有其他分支都是背景分支，默认标签#0(并分配比例ω0)。注意，用户不需要将默认标记包含在树中。用户可以手动添加标签来定位Newick树中的前景分支，使用内部脚本或代码片段，例如我们在GitHub存储库(https://github.com/abacus-gene/paml-tutorial/tree/ main/ posiposi_selection /00_data)中提供的代码片段，或者使用图形工具，如PhyloTree (Shank等人2018)或EasyCodeML (Gao等人2019)轻松定位和标记前景分支。有些工具可能会输出带有不同标记的Newick树，与CODEML中使用的标记不同，因此用户应该在运行CODEML之前检查标记格式。

        在这里，我们进行了四项测试，指定不同的分支作为前景:(1)鸡分支，(2)鸭分支，(3)鸡和鸭分支，(4)鸡和鸭分支以及有根树中鸡和鸭的分支祖先(即后续通常称为鸟类分支);表5)。零模型在所有测试中都是M0，它假设所有分支都是相同的ω。除检验4中的备择假设外，所有模型都使用无根树，在备择假设中，树的根周围的两个分支被赋予不同的ω比，因此需要有根树。详细信息请参见补充材料中的“有根树与无根树”一节。

    在第一个测试中，我们评估了鸡系是否处于正选择状态。树文件如下所示，用标记#1将鸡分枝标记为前景分支：

        然后，我们在其他三个假设下进行同样的分析。此分析的控制文件与之前类似，除了以下变量:

        通过设置model = 2，我们明确了前景分支可能在不同于背景分支的ω值下进化。在运行其他三个分析时，请注意树和输出文件名应该相应地更改。我们可以将更新后的控制文件保存为codeml-M0_branch.ctl并从命令行执行CODEML:

        输出具有与以前相同的格式，唯一的区别是输出文件中有一个额外的块，其中包含输入树的每个分支的估计dN/dS比率。例如，假设鸡谱系作为前景分支的分析输出文件包括以下内容:

        在此树中，每个分支的估计ω比显示为分支标签，在“#”符号之后。这里，所有背景谱系ω0 = 0.293，而前景谱系ω1 = 1.171（即鸡）。

        在将鸭子谱系标记为前景分支的分析中(表5，测试2)，前景分支的估计为ω1 = 999，背景分支的估计为ω0 = 0.300。数值999是程序中设置的上限，表示无穷大。如果相关分支缺乏同义词替换，则可能发生这种极端估计(Hou et al. 2007)。注意，在这种情况下，即使很难估计ω1的精确值，LRT仍然有效。当鸡和鸭谱系都被标记为前景分支(测试3)时，估计值为ω0 = 0.287和ω1 = 0.711，当鸟类分支的三个分支都被标记为前景分支时，我们得到的估计值与测试4相同。结果表明，鸡和鸭两个世系的ω比值高于其他世系，表明可能存在正选择。

        为了检验我们假设的重要性，我们可以使用LRT(表5)。根据LRT，我们得出结论，分支模型比M0模型更适合所有测试假设的数据。换句话说，在每个假设下测试的谱系(鸡、鸭或鸟类进化枝)的ω比率与背景分支的ω比率显著不同。

        在这里，我们展示了(1)如何在分支模型下执行CODEML，当对树上不同的进化谱系或分支假设不同的选择压力(即不同的ω比)时，以及(2)如何使用M0作为零模型进行影响预先指定分支的正向选择的LRT。我们进行了四个测试，部分是为了说明，如果模型为根周围的两个分支指定了不同的进化过程，则可以使用有根树。注意，测试中前景分支的识别取决于被测试的生物学假设，这应该是先验的。如果在没有任何先验假设的情况下，将树的每个分支依次指定为前景进行分支测试，则需要对多重测试进行校正。

############################################################################

先验假设应该怎么理解？做出的先验假设必须正确，还是不一定正确？既然是假设，就存在错误的可能性，推断出一个概率判断错误的可能性，错误的可能性越小，就能拒绝假设。又为何对每个分支进行前景测试？概率论知识亟需补充，遇到这方面的不好理解。（胡言乱语）

############################################################################

15、跨分支和位点的异构分支-位点模型（branch-site model，大部分文献选择最多）
        在本节中，我们在分支位模型下运行CODEML，其中ω假设在谱系和位点之间都是变化的(图3D)。这种模型可用于检测沿预先指定的前景分支影响特定氨基酸位点的正选择。

       备择假设： 树文件将具有与分支模型相同的格式，前景分支标记为#1。我们使用前面为分支模型创建的树文件，并按如下方式编辑分支模型的控制文件:

        输出包括一节，列出了分支-位点模型Ma中假设的k = 4个位点类别的比例估计值以及背景和前景分支的ω值(表6)。位点类别0：从0到1变化(0 <ω0 <1)，该类位点正在对背景和前景分支进行净化选择。在第1类位点中，ω总是固定为ω1 = 1。这些地点处于中性进化状态。在位点类别2a或2b中，前景分支进行ω2≥1的正选择，背景分支进行0<ω<1的净化选择(2a类站点)或ω1 = 1的中性进化(2b类站点)。

        ω的估计值(表6)表明，鸡和鸭谱系中存在正选择的位点。事实上，正选择似乎影响了鸟类进化中的所有三个分支(图1)。

       null 零假设：在一个新目录中，我们复制上面使用的控制文件，但更改以下变量，使位点类别2a (ω2)的前景分支的ω值不被估计，而是固定(fix_ω = 1)到1 (ω = 1):

        该程序将像以前一样估计ω0和站点类的比例。我们可以将新的控制文件保存为codeml branchsite_null.ctl并执行CODEML:

        我们从上面的分析中收集零假设和备择假设下的对数似然得分进行LRT。结果如表7所示，鸭和鸡黏液病毒基因均处于正选择状态。

        在相应的输出文件中报告了所有假设的BEB分析结果。表8显示了谱系中被选为前景的位点，其被正选择的概率大于95%或大于99%。

        在这里，我们展示了(1)如何在分支-位点模型Ma下执行CODEML (Yang et al. 2000;Zhang et al. 2005)和(2)当使用ω2 = 1的分支-位点模型Ma作为空模型(即少一个自由参数)时，如何运行LRT进行正选择。四项正选择试验的结果是明确的:鸭和鸡谱系(以及鸟类进化的所有三个分支)的黏液病毒序列都处于正选择状态。

16、基因组尺度数据集分析
        如今，对来自同一组物种的数千个蛋白质编码基因进行相同的CODEML分析是很常见的。我们认为，在大规模运行这种分析之前，用户理解模型假设和测试的基本原理是至关重要的。例如，首先详细分析一个基因排列可能会有所帮助，正如我们上面的协议所解释的那样。

        我们修改了CODEML，使用来自同一组物种的多个基因比对进行相同的正选择测试，允许某些物种可能缺少某些基因。所有的基因序列将在相同的输入序列文件中，一个接一个，基因序列的数量使用控制文件中的ndata变量指定。树文件中提供了所有物种的主树，CODEML通过修剪基因的缺失物种来生成每个基因比对的树。

        具体来说，变量ndata是用四个选项实现的，这里简要介绍一下。我们在PAML发行版的examples/ndata/文件夹中提供了说明性示例，在我们的GitHub存储库(https://github.com/abacusgene/paml-tutorial/tree/main/positive-selection/02_ extra_analyses)中包含了一个带有其他示例的教程。选项如下，所有选项都假设一个示例序列文件有三个对齐:

Case a: ndata= 3                           #3个序列比对，共用1个物种树，3个序列比对的物种都在物种树中存在

Case b: ndata= 3 separate_trees  #3个序列比对，分别对应3个物种树

Case c: ndata= 3 maintree 1        #3个序列比对，共用一个主树；其中3个序列比对中，存在丢失基因，需要从主树中修剪去除该物种，计算ML重新建树；对于多物种基因组比较分析，会大量存在基因丢失现象；基因丢失，修剪树可能会改变树枝分布，或合并树枝，容易报错，谨慎使用

Case d: ndata= 3 maintree 0        #与c相似，没有ML过程，但最终结果有一个树文件

        在情况a中，序列文件有三个对齐，所有这些都使用相同的树进行分析(即，一个树块，由标题和Newick格式的树组成)。此选项仅在所有基因比对具有树中存在的所有物种或序列(即没有丢失数据)并且具有相同名称的情况下才有效。

        在情形b中，每个对齐都有自己的树块。CODEML将读取序列文件中的第一个基因比对和树文件中的第一个树块以运行分析。然后，它将移动到下一个对齐和下一个树块，并将继续，直到所有对齐被分析。在这种情况下，某些物种可能缺少某些基因，并且物种或序列名称可能在比对中不同。不同基因序列和不同树的输出信息分别由“Dataset 1”和“TREE #1”这样的标题分开。

        在情况c中，树文件只包含一个主树，它包含序列文件中存在的所有物种。该程序将逐一读取基因序列，并通过从主树上修剪掉缺失的物种来生成对应于每个基因序列的子树。然后，它将运行每个基因序列的ML分析，估计模型参数并计算对数似然值。

        情况d与情况c相同，只是CODEML不执行ML分析。相反，CODEML将把这些基因树写在一个名为genetrees.trees的文本文件中。与案例c相比，案例d是一次演练，用于检查序列文件中的错误。

        如果所有物种都出现在序列文件中的每个基因序列中，那么所有四种情况应该是相等的，情况a将是最简单的。

        当ndata用于分析分支和分支-位点模型下的许多基因比对，但某些物种缺少特定基因时，用户应谨慎使用。在这些模型下，使用标签将主树的分支分为前景分支和背景分支。然而，为了一个特定的基因而从主树上修剪掉缺失的物种可能会导致主树上的两个或多个分支合并成一个分支。如果这些合并的分支具有相同的标记，CODEML将保留结果分支的标记。但是，当这些分支的标记不同时，就会产生歧义，在这种情况下，CODEML将终止并显示错误消息。关于这些问题的示例在我们的GitHub存储库(https:// github.com/abacus-gene/paml-tutorial/tree/main/positiveselection/02_extra_analyses)中有进一步的讨论。

讨论
检测正选择的系统发育方法
        在本协议中，我们重点研究了在CODEML中实现的用于检测正选择的四种类型的密码替代模型:(1)所有位点和分支都有一个ω的one-ration模型，(2)编码序列中不同位点的ω异质的位点模型，(3)系统发育分支之间不同ω比的分支模型，以及(4)分支和位点之间ω异质的分支-位点模型。

        在这里，我们简要概述了用于检测蛋白质编码基因正选择的其他系统发育软件工具(见表9这些工具的选择)。请注意，我们的方案不包括用于检测选择性清除(slection sweep)的群体遗传学方法，这种方法适用于从同一物种的多个个体中取样的基因组数据(Nielsen 2005)。

        许多作者已经开发了模型来计算蛋白质中氨基酸位点之间的可变选择压力，类似于CODEML中的位点模型。ML方法为基因中的每个密码子位点分配和估计自由ω比率，并使用一个位点的数据应用零假设ω = 1的LRT。该方法之后是SLR项目中的Massingham和Goldman(2005)的现场似然比方法或HyPhy (Kosakovsky Pond et al. 2005)的固定效应似然(FEL)方法;参见铃木(2004)。然而，这个过程的一个困难是它在模型中使用了许多参数(每个位点有一个ω参数)，但是如果在比对中有大量的序列可用，它可能是有效的。

        HyPhy程序实现了几个类似于CODEML中的位点和分支-位点模型的模型。HyPhy中实现的位点模型(Kosakovsky Pond et al. 2005;另见Mayrose et al. 2007)使用三个位点类别进行同义替换速率，以及三个位点类别进行非同义替换速率，这允许同义替换速率在位点之间变化，一些人认为这一功能很重要(Wisotsky et al. 2020)。使用Nielsen和Yang(1998)的NEB方法计算每个位点属于这些位点类别的后验概率，允许ω >1。HyPhy中的进化混合效应模型(MEME) (Murrell et al. 2012)和自适应分支-位点随机效应(aBSREL)模型(Smith et al. 2015)允许ω比值在谱系和位点之间变化，类似于本协议中讨论的分支-位点模型。然而，这些模型不需要预先指定前景分支，而是寻找一些位点似乎处于ω >1. 这种方法可能对探索性数据分析有效，但在解释LRT时应谨慎(Murrell et al. 2012)。请注意，在标准假设检验中，零假设应该是先验地制定的，而不应该从数据中推导出来，而数据反过来又用于检验假设。

        FastCodeML中提供了分支-位点模型A的更有效实现(Valle et al. 2014)。并行版本可能比CODEML效率高100倍，使得分析大型系统基因组数据集成为可能。一些具有友好图形用户界面(gui)的在线工具和资源已经将CODEML集成为检测正选择的工具。例如，EasyCodeML (Gao et al. 2019)是CODEML的包装器，因此用户不必在命令行运行CODEML。类似地，e- evolution在CODEML和SLR中使用预先配置的进化模型以并行方式自动化分析，以及内置的LRT和绘制结果的工具(Huerta-Cepas et al. 2016)。

 假设、限制和展望
        本协议中讨论的分支模型、位点模型和分支-位点模型是最简单的密码子替换模型。它们通常用于基因组序列数据的分析，以检测在分子水平上可能处于正选择驱动适应的基因或位点。我们的协议是针对初学者的，我们不试图提供一个全面的主题覆盖。感兴趣的读者可以参考许多讨论这些模型的理论，局限性和应用的评论，例如Yang和Bielawski (2000)， Nielsen (2001)， Anisimova和Kosiol (2009)， Yang(2014)的第12章，以及Cannarozzi和Schneider(2012)编辑的书，以及CODEML程序手册和我们的GitHub存储库，其中包含再现本协议中描述的结果所需的代码和数据。

        像任何假设的统计检验一样，这里讨论的正选择检验可能存在两种类型的错误:假阴性错误(或缺乏权力)和假阳性错误。已经进行了一些模拟，以检查基于本协议中讨论的站点模型和分支站点模型的两种类型的测试误差(例如，Anisimova和Yang, 2007年)。请注意，所有正选择测试都考虑到过度的非同义替代——超出了基于遗传密码表和DNA水平突变过程特征的预期——作为正选择作用于给定蛋白质的证据。考虑到对任何功能性蛋白质的限制，偶发性正选择可能无法将非同义替换速率提高到足以检测到信号的程度，在这种情况下，测试将具有低效。当位点或分支-位点模型下的测试返回一个不重要的结果时，测试能力的缺乏总是可以解释这样的结果。同样，当基础模型的假设严重违反时，测试可能会产生过多的假阳性(Anisimova and Yang 2007)。

        首先，所有模型都假定对齐是正确的。请注意，已经注意到过度的对准误差会导致模拟中基于位点模型的测试误报(Fletcher and Yang, 2010)。特别是，一些比对方法倾向于将非同源位点合并到一列中，在位点上产生明显的非同义替换，因此在推断基因比对之前删除难以比对的区域可能是谨慎的。或者，可以使用插入和缺失模型以及密码子替换模型来对齐序列并联合拟合密码子模型，如BAli-Phy (Redelings 2021)程序所示。不幸的是，对于大型数据集，这在计算上可能不可行。有趣的是，经历适应性进化的具有许多非同义替换的蛋白质区域(例如，HIV env基因的高变区)也被注意到具有许多插入和缺失，这表明插入和缺失突变也可能处于正选择下。

        其次，常用的密码子取代模型假设一次只改变一个核苷酸的简单突变或取代。同时改变两个或多个核苷酸的复杂取代已被注意到会导致过多的假阳性(Kosiol等人，2007;Jones et al. 2018;Venkat et al. 2018)。

        第三，基于位点的检测被证明对同一基因位点之间的低水平重组是稳健的(Anisimova等人，2003)，尽管在高重组率下假阳性率可能非常高。

        最后，有偏差的基因转换也可能导致假阳性(Galtier and Duret 2007;Ratnakumar et al. 2010)。在进行正选择测试之前，已经做了一些努力，试图过滤掉受比对错误、重组或基因转换影响的数据。除了模拟分析，正选择测试也被应用于众所周知的加速进化的基因，如HLA (Hughes和Nei 1990)。在某些情况下，从统计检验中产生的生物学假设促进了实验验证，提供了令人兴奋的分子适应案例研究(例如，Sawyer et al. 2005)。

        自从二十年前密码子替换模型被引入以来(Goldman and Yang 1994;缪斯和Gaut 1994)，它们一直在积极发展，并在许多方面得到扩展。这些模型为研究基因序列进化中自然选择的方向和强度以及检测分子适应提供了重要的框架。它们还提供了更现实和更具解释性的氨基酸替换模型，可用于推断系统发育树或重建祖先蛋白质序列(Yang et al. 1998)。最近的研究利用突变选择(mutt - sel)公式来提高生物现实性和模型的解释力(Halpern和Bruno 1998;Yang and Nielsen 2008)。特别是，mt - sel密码子模型已被用于通过使用大型系统发育数据集来估计单个氨基酸位点的选择系数分布(Rodrigue et al. 2010;Tamuri et al. 2014;Tamuri和dos Reis 2022)，以估计突变偏差(Latrille和Lartillot 2022)，并测试基因序列进化速度与物种生活史性状之间的关系(Latrille等人，2021)。通过明确考虑突变偏差、同义和非同义突变的固定概率，mutt - sel模型具有更大的解释力，并具有整合群体遗传学和系统遗传学的潜力，以进一步了解蛋白质编码基因的进化过程。

###########################################################################

补充材料

##########################################################################

基因树与物种树
        有时，基因树(例如，使用基因比对推断的ML树)和已建立的物种树可能不同。正选择分析应该基于基因树还是物种树?这个问题没有一个简单的答案。请注意，测试中使用的模型假设系统发育树代表了序列的真实进化关系。因此，应该使用最可能正确的树。在同源基因和相似基因的重复分析中，物种树可能不适用，因此推断的基因树是唯一的选择。同样，在病毒序列分析中，物种树不存在，因此基因树是唯一的选择。如果基因序列较短或不包含太多的系统发育信息，则推断出的基因树可能无法解决或不正确，而物种树将更可取。如果趋同进化可能会误导基因树的重建，那么物种树将是更好的选择。

        当系统发育树有疑问时，建议使用几种可能的树(例如，包括基因的ML树和物种树)来评估树拓扑的影响。在模拟分析中，发现基于位点的正选择测试对于系统发育的微小变化是稳健的(例如，Yang等人，2000年)。在分支或分支-位点测试中，如果前景分支是很好的谱系，并且系统发育的不确定性涉及指定为背景分支的分支内部的细节，那么树的拓扑结构可能不会对测试产生重大影响。有关确保正选择推论质量的其他检查，请参阅Álvarez-Carretero和dos Reis 2020)。

有根树与无根树

        在PAML中，有根树在根处使用三分岔表示，而无根树在根处使用二叉树表示。例如，在图S1A中，左边的树“(A, B)， C))”是一棵有根的树，有四个分支长度，其中根周围有两个分支长度，而右边的树“(A, B,C);”是一棵无根的树，有三个分支(根周围的分支长度b3a和b3b合并为一个分支长度b3)。可以使用文本编辑器或各种脚本删除Newick表示法中的一对括号，从而将有根树转换为无根树。我们在GitHub教程(https://github.com/abacus-gene/pamltutorial/tree/main/positive-selection/00_data#readme)中包含了一个代码片段。

        在分析中是否使用有根树或无根树取决于替代模型能否识别出树的根。特别是，如果替代模型是时间可逆的，则替代过程是时间均匀的:核苷酸、密码子或氨基酸的频率是固定的，不同的谱系有自己的速率(即，没有分子钟的假设)。在这种情况下，根的位置无法识别，应该使用无根树。实际上，所有的系统发育程序，如RAxML或IQ-TREE (Minh et al. 2020)都假设了时间可逆的替代模型，并且没有时钟。因此产生了无根的树。即使树形绘制软件(例如FigTree)可能出于视觉目的而显示有根的树，如果在无法识别根的模型下推断出该树，则应认为该树无根。

        几乎所有在文献中发展的密码子模型，包括这里讨论的，都是时间可逆模型。此外，我们不假设分子钟在整个方案中描述的分析。因此，在使用分支或分支-位点模型时，通常应该使用无根树，但有一个例外(参见协议中的示例)。在这种情况下，如果我们假设根周围的两个分支正在经历不同的进化过程(例如，具有不同的ω)，则根的位置是可识别的，并且应该使用有根树。如果假设根周围的两个分支按照相同的过程进化(例如，两个分支都是前景分支或两个分支都是背景分支)，则根是不可识别的，并且应该使用无根树。在模型M0(one-ration)和位点模型(如M1a、M2a、M7、M8)下，总是假设根周围的两个分支按照相同的过程进化，因此应该使用无根树。图S1展示了几种场景，其中黑色分支代表前景分支，灰色分支代表背景分支。

        在运行任何分析之前，用户可能会问:当应该使用无根树时，错误地使用有根树会产生什么影响?为了回答这个问题，让我们考虑将模型M0 (1 -ratio)拟合到图S1A中左侧的有根树。所有模型参数，如替换/颠换率比(κ)、非同义/同义率比(ω)、分支长度b1和b2将被识别和正确估计，并将正确计算对数似然值。分支长度b3a和b3b是不可估计的，尽管它们的和b3 = b3a + b3b是可估计的。如果多次运行CODEML, b3a和b3b的估计值可能会随运行而变化，但b3的估计值将是稳定的。如果我们对零假设(ω= 1)进行LRT，并在零假设和备选假设中都使用有根树，我们将把参数的数量多计算1(即，用于对树进行有根的额外分支长度)，但自由度将被正确计算，并且LRT仍然是正确的。例如，如果在M1a和M2a下都使用了有根树，即使参数的数量多了一个，LRT也会是正确的。请注意，此场景也适用于位点测试。在这些情况下，理想情况下应该使用无根树，尽管使用有根树不会像前面解释的那样造成任何严重损害。

        如果我们现在考虑图S1D中指定的分支或分支-位点模型，假设根周围的两个分支具有不同的进化过程。该模型只能通过使用有根树来表示，因为使用无根树将指定一个完全不同的模型。

BEB分析
        在协议中使用的示例中，BEB分析没有列出任何概率大于95%或大于99%的正选择位点。下面，我们展示了一个例子，说明在上述限制下，如果一个站点被正选择，输出将是什么样子:

        在这个虚构的场景下，序列中的第10个位点的后验概率为98.4%，来自正选择类𝜔>1 ,该位点𝜔的近似后验分布均值为1.468,SD为0.634。同样，位点25的后验概率为99.9%，来自正选择类，𝜔的近似后验均值为1.464,SD为0.638。在CODEML输出中，后验概率P >0.95以*和P >0.99用**表示。

参考文献：

Beginner's Guide on the Use of PAML to Detect Positive Selection。

############################################################################

疑问：

1、位点存在正选择，是否可以认为基因处于正选择？应该是这样。

2、指定前景支，其余皆为背景支；如果不想要这么多背景支，是否必须修剪去除一些物种才可以做到？