摘PAMl学习笔记 - 蜗牛

PAML Branch model

以灵长目动物溶菌酶编码基因适应性进化分析为例解读——Branch Model 与1.什么是Branch model？

Branch model是PAML软件CODEML程序中通过likelihood ratio test （LRT）进行不同支系间（lineages）适应性进化检测的一种模型。该模型通过限制（constraint）系统发育树中不同分支上的omega(dN/dS)值的异同，并对不同的限制进行显著性分析（PAML软件中的Chi2程序），进而得到较为可靠地分析结果。 在该法提出之前，不少学者通过简约法（parsimony method）或者似然法（likelihood method）先重建祖先序列（ancestral sequence），然后通过对构建的祖先序列的omega值估算进而预测不同支系的适应性进化特征。诸如Prof. Messier等对于灵长目动物溶菌酶的分析便是如此。Prof.Yang认为，从统计学的角度而言，这种将预测的数据当做真实观测数据的分析理念存在一定的随机误差（random errors）和系统误差(systematic errors), 本身并不是一种严谨的统计学方法。

Prof.Yang 所提出的Branch model巧妙地避开了直接利用

ancestral sequence进行支系间适应性进化检测的流程，而是通过平均统计每一个节点（each node）中可能的ancestral sequence，根据其相对发生似然率（relative likelihoods of occurance）进行加权分析。此外，Branch model还考虑到了（take into account）密码子转移/颠换速率偏差（transition/transversion rate bias）和非均匀密码子（nonuniform condon usage）这些与omega值计算有着显著关系的影响因素。

2.Branch model中存在哪些假设模型，在CODEML程序的control file文本中如何选择？

Branch model主要是对系统发育树中的不同支系的omega值的异质性进行界定，主要的model有：one-ratio model，即系统发育树中所有支系的omega 值是相等的；free-ratio model，该模型指的是系统发育树中所有支系的omega值是不相等的。这两个假设是不同支系omega取值的两个极限。此外，还可以设定前景枝（foreground clade）,假定其与其余支系（又称背景枝 background clade）的omega值不同。前景枝可以根据需要设置多个。

在control file中，Model=0 表示one-ratio model， Model=1表示free-ratio model. Model=2表示系统发育树中不同omega值得个数，其中所选择前景枝的个数为(n-1)。值得注意的是，当设置

Model=2,3，……，n时，需要在tree file中标记所要设置的前景枝，可以标记一个，也可以标记多个。树标记格式如下所示：

((1, 2), 3) #1, 4, 5); 该tree file表示Clade 1，2 and 3为前景枝，其对应的omega值为ω1（用#1表示），其余Clade为背景枝，对应的omega值为ω0（用#0指定，但在PAMl软件中，#0为默认值，故不需要在树中注明）。在result file mlc 文件中，我们可以得到两个不同的omega值。

3.通过Branch Model 可以得到什么样的结论？ 3.1不同支系间的omega值是否显著不同

这主要通过比较one-ratio Model 和free-ratio Model对应的likelihood values的差异进行说明。 3.2 前景枝和背景枝的omega是否显著不同

这主要通过比较one-ratio Model和two-ratio Model的likelihood values 进行说明。

3.3 前景枝的omega值是否显著大于1（greater than one）

这主要通过比较two-ratio Model中存在与不存在ω1<1这一约束的两种情况下所对应的likelihood values进行说明。 4. 如何对不同Model 的差异性进行比较？

该比较主要在PAML软件Chi2程序中进行，首先在mlc文件中查找不同Model所对应的likelihood values并计算不同Model

likelihood 差值绝对值的两倍，即2△l=|l1-l2|。打开Chi2程序，界面如图1所示。

图1 Chi2程序主界面 Fig1 The main of Chi2 program

通过上图查询df=10时，2△l值所对应的显著性水平，小于等于0.05时，被认为是存在显著性差异的, 如图1中绿色框所标注。 注：该法与linxiao.name（http://linxiao.name/archives/172） 网站中所述方法有异（将df值与2△l值输入程序中，回车查看显著性水平），但网站中是在Linux平台下操作，而本法是在windows平台下操作。另外，在Chi程序的windows版本中并未发现任何输入的光标。

5.Prof.Yang对灵长目动物溶菌酶不同支系的适应性进化分析 5.1数据和方法（Data and methods） 5.1.1数据（Data）

本文所涉及的数据主要分为两部分内容，首先是大数据集（large dataset），这主要包括Prof.Messier分析的24条灵长目动物溶菌酶编码基因中有显著不同的19条序列（Distinct sequence），其系统发育树如图2所示。

6.langur Sen&Sve 7.langur Tob&Tfr

叶猴

h 10.baboon Pcy 11.mangabey Cat 8.Douc langur Pne 9.probiscis Nla 5.colobus Cgu&Can 颊囊猴 12.rhesus Mmu 13.Allen Ani 14.talapoin Mta 15.patas Epa 16.vervet Cae 1.humanc 2.chimp bonobo gorilla 3.orangutan Ppy 人科动物 17.squirrel m 4.gibbon Ggo

18.tamarin Soe

19.Marmoset Cja

新大陆猴

0.02

图2 Prof.Yang 所分析的灵长目溶菌酶系统发育树 Fig2.Phylogeny of the primate lysozyme analyzed by Prof. Yang

在图2中，Branch h和Branch c是本文分析是所选取的前景枝（foreground clade）。文中这两个前景枝的选取是根据Prof.Messier于1997年的研究结果。

小数据集包括7条序列，其来源是从图2中四个分支中各挑选出几条具有代表性的序列，重新进行分析，其系统发育树如图3所示。

图3 从图2挑选出来的四个分支代表序列的系统发育分析

Fig3 Phylogeny of a subset of seven primate lysozyme selected from those of Fig 1. to represent the four

major groups of species.

Prof.Yang认为，对于这种大数据集和代表性的小数据集的分析比较能够之时取样方法对于所得结果的一个敏感度。（原文：Differences between the

analyses of the two data sets will give us an indication of the sensitively of the results to species sampling.） 5.1.2方法（Methods ）

第一，对小数据集control file文件的解读

图4 溶菌酶小数据集的control file Fig 4 The control file of lysozyme small dataset

在图4中，seqfile，treefile 和outfile 分别指代的是分析所用的序列文件，树文件以及主要结果输出文件。

Noisy命令控制屏幕中输出文件的数量，一般取值为9； Verbose命令控制文件中输出结果的数量；

Runmode命令控制是否使用树文件，常见的值有0， 1， 2， 3， -2等，常见取值为0，表示从树文件中获取树的拓扑结构，-2是pairwise alignment的命令，若选择runmode=-2，则不需要输入树文件。

Seqtype命令是指代输入系列文件的类型，可取值是1 ，2， 3，seqtype=1表示序列文件为condons sequence，seqtype=2表示序列文件为amino acid sequence，seqtype=3表示序列文件为通过condons sequence 翻译得到的amino acid sequence；

CodonFreq命令是指定密码子使用频率模型的类型，其所取数值0,1,2,3分别代表不同的密码子使用频率模型。一般取值为CondonFreq=2, 有人对不同

CondonFreq的取值进行了比较，发现CondonFreq=1和CondonFreq=2并未有显著性的差别。

Clock命令用来指代不同支系间的速率（指的是什么速率？）是恒定的还是变化的（原文：clock specifies models concerning rate constancy or variation among lineages）。可取的值有0，1，2 和3，其中0表示支系与支系之间的速率是不断变化的；1表示所有支系的速率相同；2表示系统发育树中大多树支系的速率是相同的，只有少部分特意指定的支系速率不同；3用于对多个基因或者多分区数据进行联合分析时所选取的参数。Runmode =1 or 2时，tree file必须是rooted tree。在常见的分析中，clock=0的命令居多。

Model命令主要是对系统发育树中不同支系的omega值进行限定。可取值有0,1,2。0表示系统发育树所有支系均有相同的omega值，1表示omega值是自由的，不同支系间各不相同。2表示指定的前景枝和背景枝之间的omega值不同，一般而言，有n个前景枝，便有n+1个omega值。对于Branch Model，一般选择Model=2，但在该模式下，需要对树文件进行前景支的定义。另外，对于小数据集的Branch Model分析还可以选择Model=1，但由于Model=1所涉及的参数数量较多，不建议对大数据集进行Model=1的分析。

NSsites是在Site Model和Branch-site Model中进行设置的参数，在Branch Model中一般设置为0。

icode用来指定遗传密码子的类型, 可取值有0——10 共11个数字，0表示通用密码子，也是最为常见的一个设置。该命令主要根据所分析的序列来源进行具体的设置。

fix_kappa命令用来指代K80，F84和HKY85中的kappa值是一个固定值还是根据数据进行迭代分析得到的。Fix_kappa=1表示kappa值是一个固定值。Fix_kappa=0表示kappa值是通过数据迭代运算得到的。而kappa是定义一个kappa的初始值，这个需要使用者自行设置。

Fix_omega and omega，fix-alpha and alpha 以及fix-rho and rho等的设置参照fix_kappa进行。

注：在本文中，对于图5所涉及的一些模型，需要改变fix_omega and omega值进行运算。

图5 本文中的一些模型 Fig5 some models of this paper

若要限定omega值等于1 ，则需要在control file中设置fix_omega=1 and omega=1. getSE命令是来设置是否需要计算估算参数的标准误（strand errors）,若需要则getSE=1，不需要则getSE=0. Rateancestor=0 or 1, 该值常设为0。

Method命令主要控制在no clock 模型下计算枝长的迭代算法。Method=0表示该迭代算法为老式算法，Method=1表示该迭代算法为新式算法，注意当clock=1,2,3时，Method只能取0.

以人类与PAML sites model 病毒适应性进化分析为例解读——HIVSites Model 1.什么是Sites Model？

Sites Model是PAML软件CODEML程序的一个正选择作用分析模型，其主要观点是同一序列不同位点的omega值不同。在进行sites Model分析时，需要设置control file中的Model=0，Nssites命令在此是一个变量，根据不同Model的选择设置不同的值。值得注意的是，以此可以选择多个sites Model。如Nssites=0 1 3 7 8.

2.不同的sites Model 表示什么意思？

? M0即one-ratio Model，值得是所有位点的omega值是恒定的；

? M1表示加假定有一部分位点的omega值为0，其他位点的omega值为1； ? M2是在M1的基础上增加了第三类omega值，该类omega是通过数据计算

得到的，有可能大于1； ? M3假定所有位点的omega值呈简单的离散分布趋势； ? M5假定所有位点的omega值呈简单的gamma分布趋势； ? M7假定所有位点的omega属于矩阵（0,1）且呈beta分布；

? M8是在M7的基础上增加另一类omega值，该值可通过计算得到，可以大

于1. ? M8a与M8类似，只不过新增加的omega值等于1. 3.不同Model的比较可以得到什么样的结果？

首先是M0与M3的比较，该比较与Branch Model中的Model之间的比较是一样的，首先计算2△l值，然后在一定df值下进行显著性水平计算。这里需要注意的是，在参阅了Prof.Yang关于PAML的一些参考材料之后，我们发现sites Model比较的df值一般取2.

在sites Model 中，M0表示one ratio for all sites， M3表示所有位点的omega值呈简单的离散分布。对于这两个模型的比较并非用于正选择作用的检测，而是用于位点间omega值是否一致的检测。

M1 and M2 以及M7 and M8是用于正选择作用的检测，但Prof.Yang认为，The M1-M2 comparison 与 the M7- M8 comparison相比，更加的稳定。（原文：The M1-M2 comparison appears to be more robust (or less powerful) than the M7- M8 comparsion）

此外，还有一类比较是M8 to M8a，其中M8和M8a是两个极为类似的Model。在涉及到positive selection 的文章中，对于这两个model的比较并不常见，而且说明书中也并未给出明确的比对结果意义。 4. 如何检测positive sites？

在CODEML中，positive sites 的检测流程主要如图1所示.

PP value computation Likelihood ratio test CODEML computation

图1 positive sites的检测流程 Fig1 The process of positive sites

其中CODEML computation主要是对control file中的命令值进行设定之后，运行CODEML程序，并在result file中查看运算结果。Likelihood ratio test如question3所示，即对两个模型进行显著性水平比较。PP value computation主要是指位点后验概率的计算，该结果是显示在main result file- mlc文件中。CODEML程序中常见的计算后验概率的方法有BEB和NEB。与BEB相比，NEB在计算的过程中往往会忽略抽样误差。因此，Prof.Yang建议在读取运算结果时，可以直接将NEB result忽略，但值得注意的是，BEB只能在M2a和M8 model下运行。

5.以example中control file文件为参考，解读site model下的control file。

图2 site model下的control file 截图

Fig2 The fig of control file under site model

Site model 计算的control file与Branch model中大致相似，但在site model中应当注意，model=0是一个恒定值，Nssites命令可以设置不同的模型参数。

PAML

与Branch sites model 1. 什么是Branch-site model？

Branch- site model其实是Branch model and site model的合集，在该model下，不仅假定位点间的omega值是变化的，同时也假定支系间的omega值是变化的。该model主要用于检测前景枝中正选择作用对部分位点的影响。 2. Branch-site model中都有哪些模型？

Model A=（model=2 NSsites=2）；Model B= （Model=2 Nssites=3）；Model C= （Model=3 Nssites=2）；Model D= （Model=3 Nssites=3）。注意Model D 需要ncatG参数设置位点的类型（以omega值进行分类），可以使用的ncatG值有2或者3，在ModelA，ModelB and ModelC中，ncatG的取值是被自动忽略的。 3. 在Branch- site model中如何进行模型比较？

Model A与ω2=1的null Model进行比较，该null Model可以通过设定

fix_omega=1 and omega=1进行确定。The comparison of Model A and null Model 对应的df值为1；（注：在null Model中omega值的设定是ω2=1，这便限定了在null Model中，background中的sites均处于negative selection，foreground Branch 中的sites均处于neutral selection。而在alternative Model中，foreground Branch中的site对应的omega值是大于1的。因此，alternative Model and null Model之间的显著性水平检测可以直接用来检测foreground Branch中的positive selection 位点）

在Branch- site Model中，Model A还可以与site Model中的M1a进行显著性水平比较。若这两者之间存在显著性差异，则可以说明要么foreground

Branch受到的选择约束较为宽松，要么foreground Branch受到明显的正选择作用。

Model B通常与site Model中的M3进行显著性比较，对应的df值为2 Model C与site Model中的M1a进行比较，对应的df值为3.

Model D通常与site Model中的M3进行比较，若树文件为无根树，则df=1，若树文件为有根树，则df=2.

就Model的提出理念而言，Model A and Model B侧重于寻找在进化过程中受正选择作用的点，而Model D则与之不同，其不再局限于正选择作用，而是涉及到多种选择作用(divergent selective pressures)。 4. 在Branch-site Model中的其他注意事项

在Model C and Model D中，不同omega值的Branch类型不局限于两种，使用者可以自行设置多种Branch types，最多可以设置10-12种。

另外，在Model C中，其ω1是一个定值，而在Model D中，ω1则是一个自由变化的参数。

Maximum likelihood methods for detecting adaptive

evolution after gene duplication

导读：基因组数据的不断报道使得基于大数据分析的基因家族进化成为可能，在此基础上，本文提出了一种基于Maximum likelihood methods的方法，对基因家族形成过程中的适应性进化进行检测，并以灵长目的ECP-EDN基因家族为例对此方法进行了说明。

1. 选择压力支系特异性检验（Detecting lineage- specific changes in selective pressure）

若基因家族的分化是通过正选择作用推动的，则在复制事件之后，会立刻出现非同义替换速率大于同义替换速率的现象。但是，若这一基因家族在其余的时间都是受到purifying selection 的作用，那么仅仅是两个序列之间的比对是很难发现omega值大于1的位点。

鉴于这一问题，本文提出了用于检测基因复制后适应性进化的最大似然法，旨在完成以下目标：（1）同一个基因进化历史中不时间点的选择压力；（2）这些选择压力是否不同。

1.1 支系间选择压力可变模型（Models of variable selective pressures among branches） 1.1.1.模型介绍

(1)one- ratio model ,即Phylogenetic tree中所有sites的omega值是一个恒定值。ω0=ω1=ω2=ω3；

(2)Two- ratio model, 复制事件之前的omega值和复制事件之后的omega值不相等。ω0≠ω1=ω2=ω3；

90

图1 假设基因家族的进化树

Fig1 Phylogeny of a hypothetical gene family

(3)Three- ratio model, 该模型共假设三个omega值，即复制事件之前的omega值，复制事件之后的omega值，以及后期分支中的omega值。ω0≠ω1≠ω2=ω3。

(4)Four- ratio model, 该模型共假设四个omega，不仅复制事件前后的omega值不同以外，随后的分支上的两个omega值也不相同。即ω0≠ω1≠ω2≠ω3。 1.1.2 模型比较可以揭示的问题（likelihood ratio test）

(1) one ratio model with two ratio model : 揭示基因家族的复制前后所受的选择压力是否相同；

(2) two ratio model with three ratio model: 揭示基因家族复制后到分化前以及分化后这段时间的选择压力是否发生变化。

(3) three ratio model with four ratio model: 揭示基因家族中两个旁系同源分支之间的选择压力是否发生变化。

1.1.3 应用举例——ECP-EDN基因家族在基因复制之后非同义替换速率是否发生改变。

ECP，即嗜酸性粒细胞衍生神经毒素，EDN，即嗜酸性粒细胞阳离子蛋白。二者均是核糖核酸酶类，但其特异性功能发生了分化，ECP是阳离子毒素，对于寄生虫和细菌都有非特异性毒性。但EDN是一种有效地抗病毒药物通过有效地核酸降解。本文中对ECP-EDN基因家族在基因复制之后的非同义替换率进行了检验。

图2 ECP-EDN基因家族系统发育树

Fig 2 The phylogenetic tree of ECP-EDN gene family

本分析共设置了四个ratio Model，分别为one ratio Model， two ratio

Model， three ratio Model and four ratio Model。其中R1 和R2的LRT检测结果是显著的（P=0.0001），这说明在基因复制事件之后，非同义替换速率发生了显著的增加；R2和R3的LRT检测结果也是显著的，这说明在基因复制事件之后，ECP-EDN基因家族经历了适应性检测；R3和R4的LRT检测结果也是显著的，这表明ECP subclade 与EDN subclade所处的进化选择压力不同。 1.2 鉴定适应性进化的氨基酸位点（Identification of amino acid sites under positive selection）

从蛋白质进化的角度而言，为了维持蛋白质功能的稳定性，其平均非同义替换的数量往往较少。只有少数一部分氨基酸位点在进化过程中受到适应性进化的选择。因此，计算整个支系的平均omega值往往很难检测到positive

selection。鉴于此，我们需要对小部分的condons所受的正选择压力进行检测并确定存在这些正选择压力的位点的位置。

为完成上述目标，我们可以根据已知的结构域和功能域的信息，将氨基酸位点分成若干具有独立omega值的小位点合集。但是在不了解蛋白质的结构域和功能域的情况下，可以对 all amino acid sites 设计一个omega值分布（比如gamma 分布，beta 分布等）。通过对null Model 以及其衍生模型的LRT检验，验证是否存在positive sites的假说，若positive sites存在，则通过NEB或者BEB法检验这些位点的后验概率。

1.3 Lineage- specific changes in selective pressure at specific amino acid sites Model A: 指定ω0=0 且ω1=1，因此正选择作用只约束在了前景枝中；通常选择M1 与Model A进行LRT检测，自由度df=2；

Model B：ω0和ω1为自由参数（free parameters），所分析的整个支系中均有可能出现正选择作用位点。通常选择M3与Model B进行LRT检测，自由度df=2；

本文来源：https://www.bwwdw.com/article/c9ra.html