1、PAML与 Branchmodel以灵长目动物溶菌酶编码基因适应性进化分析为例解读BranchModel1.什么是 Branch model？Branch model 是 PAML 软件 CODEML 程序中通过 likelihood ratio test （LRT）进行不同支系间（lineages）适应性进化检测的一种模型。该模型通过限制（constraint）系统发育树中不同分支上的omega(dN/dS)值的异同，并对不同的限制进行显著性分析（PAML软件中的 Chi2 程序） ，进而得到较为可靠地分析结果。在该法提出之前，不少学者通过简约法（parsimony method ）或者似然
2、法（likelihood method）先重建祖先序列（ancestral sequence ） ，然后通过对构建的祖先序列的 omega 值估算进而预测不同支系的适应性进化特征。诸如 Prof. Messier 等对于灵长目动物溶菌酶的分析便是如此。Prof.Yang 认为，从统计学的角度而言，这种将预测的数据当做真实观测数据的分析理念存在一定的随机误差（random errors）和系统误差 (systematic errors), 本身并不是一种严谨的统计学方法。Prof.Yang 所提出的 Branch model 巧妙地避开了直接利用ancestral sequence 进行支系间适
3、应性进化检测的流程，而是通过平均统计每一个节点（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 程序的 contr
4、ol file 文本中如何选择？Branch model 主要是对系统发育树中的不同支系的 omega 值的异质性进行界定，主要的 model 有：one-ratio model，即系统发育树中所有支系的 omega 值是相等的；free-ratio model，该模型指的是系统发育树中所有支系的 omega 值是不相等的。这两个假设是不同支系 omega 取值的两个极限。此外，还可以设定前景枝（ foreground clade） ,假定其与其余支系（又称背景枝 background clade）的 omega值不同。前景枝可以根据需要设置多个。在 control file 中，Model=
5、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
6、软件中，#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 值是
7、否显著大于 1（greater than one）这主要通过比较 two-ratio Model 中存在与不存在 11 这一约束的两种情况下所对应的 likelihood values 进行说明。4. 如何对不同 Model 的差异性进行比较？该比较主要在 PAML 软件 Chi2 程序中进行，首先在 mlc 文件中查找不同 Model 所对应的 likelihood values 并计算不同 Model likelihood 差值绝对值的两倍，即 2l=|l 1-l2|。打开 Chi2 程序，界面如图 1 所示。图 1 Chi2 程序主界面Fig1 The main of Chi2 prog
8、ram通过上图查询 df=10 时，2l 值所对应的显著性水平，小于等于 0.05 时，被认为是存在显著性差异的, 如图 1 中绿色框所标注。注：该法与 linxiao.name（http:/linxiao.name/archives/172 ） 网站中所述方法有异（将 df 值与 2l 值输入程序中，回车查看显著性水平） ，但网站中是在 Linux 平台下操作，而本法是在 windows 平台下操作。另外，在 Chi 程序的 windows 版本中并未发现任何输入的光标。5.Prof.Yang 对灵长目动物溶菌酶不同支系的适应性进化分析5.1 数据和方法（Data and methods ）
9、5.1.1 数据（Data）本文所涉及的数据主要分为两部分内容，首先是大数据集（large dataset） ，这主要包括 Prof.Messier 分析的 24 条灵长目动物溶菌酶编码基因中有显著不同的 19 条序列（Distinct sequence） ，其系统发育树如图 2 所示。6.langur Sen&Sve7.langur Tob&Tfr8.Douc langur Pne9.probiscis Nla5.colobus Cgu&Can10.baboon Pcy11.mangabey Cat12.rhesus Mmu13.Allen Ani14.talapoin Mta15.pata
10、s Epa16.vervet Cae1.human2.chimp bonobo gorilla3.orangutan Ppy4.gibbon Ggo17.squirrel m18.tamarin Soe19.Marmoset Cja0.02hc颊囊猴人科动物图 2 Prof.Yang 所分析的灵长目溶菌酶系统发育树Fig2.Phylogeny of the primate lysozyme analyzed by Prof. Yang叶猴新大陆猴在图 2 中，Branch h 和 Branch c 是本文分析是所选取的前景枝（foreground clade） 。文中这两个前景枝的选取是根据
11、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 认为，对于这种大数据集和代表性的小数据集的分析比较能够之时取样方法对于所得结果的一个
12、敏感度。 （原文：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 fileFig 4 The control file of lysozyme small dataset在图 4 中，seqfile，treefile 和 outfile 分别指
13、代的是分析所用的序列文件，树文件以及主要结果输出文件。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 s
14、equence，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 co
15、nstancy 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
16、值，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。
17、icode 用来指定遗传密码子的类型, 可取值有 010 共 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 a
18、nd 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 c
19、lock 模型下计算枝长的迭代算法。 Method=0 表示该迭代算法为老式算法，Method=1 表示该迭代算法为新式算法， 注意当clock=1,2,3 时，Method 只能取 0.PAML与 sitesmodel以人类HIV病毒适应性进化分析为例解读SitesModel1.什么是 Sites Model？Sites Model 是 PAML 软件 CODEML 程序的一个正选择作用分析模型，其主要观点是同一序列不同位点的 omega 值不同。在进行 sites Model 分析时，需要设置 control file 中的 Model=0，Nssites 命令在此是一个变量，根据不同Mo
20、del 的选择设置不同的值。值得注意的是，以此可以选择多个 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 分布趋势； M
21、7 假定所有位点的 omega 属于矩阵（0,1）且呈 beta 分布； M8 是在 M7 的基础上增加另一类 omega 值，该值可通过计算得到，可以大于 1. M8a 与 M8 类似，只不过新增加的 omega 值等于 1.3.不同 Model 的比较可以得到什么样的结果？首先是 M0 与 M3 的比较，该比较与 Branch Model 中的 Model 之间的比较是一样的，首先计算 2l 值，然后在一定 df 值下进行显著性水平计算。 这里需要注意的是，在参阅了 Prof.Yang 关于 PAML 的一些参考材料之后，我们发现 sites Model 比较的 df 值一般取 2.在 s
22、ites 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
23、M7- M8 comparsion）此外，还有一类比较是 M8 to M8a，其中 M8 和 M8a 是两个极为类似的Model。在涉及到 positive selection 的文章中，对于这两个 model 的比较并不常见，而且说明书中也并未给出明确的比对结果意义。4. 如何检测 positive sites？ 在 CODEML 中，positive sites 的检测流程主要如图 1 所示.图 1 positive sites 的检测流程Fig1 The process of positive sites其中 CODEML computation 主要是对 control file 中的
24、命令值进行设定之后，运行 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
25、model 下运行。5.以 example 中 control file 文件为参考，解读 site model 下的 control file。图 2 site model 下的 control file 截图Fig2 The fig of control file under site modelCODEML computationLikelihood ratio testPP value computationSite model 计算的 control file 与 Branch model 中大致相似，但在 site model中应当注意，model=0 是一个恒定值，Nssites
26、 命令可以设置不同的模型参数。PAML与Branchsitesmodel1. 什么是 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=
27、 （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 Mod
28、el 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 之间的显著性水平检测可以直接用来检测
29、 foreground Branch 中的 positive selection 位点）在 Branch- site Model 中，Model A 还可以与 site Model 中的 M1a 进行显著性水平比较。若这两者之间存在显著性差异，则可以说明要么 foreground Branch 受到的选择约束较为宽松，要么 foreground Branch 受到明显的正选择作用。Model B 通常与 site Model 中的 M3 进行显著性比较，对应的 df 值为 2Model C 与 site Model 中的 M1a 进行比较，对应的 df 值为 3.Model D 通常与 sit
30、e 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 种。
31、另外，在 Model C 中，其 1 是一个定值，而在 Model D 中， 1 则是一个自由变化的参数。Maximum likelihood methods for detecting adaptive evolution after gene duplication导读：基因组数据的不断报道使得基于大数据分析的基因家族进化成为可能，在此基础上，本文提出了一种基于 Maximum likelihood methods 的方法，对基因家族形成过程中的适应性进化进行检测，并以灵长目的 ECP-EDN 基因家族为例对此方法进行了说明。1. 选择压力支系特异性检验（Detecting lineage
32、- specific changes in selective pressure）若基因家族的分化是通过正选择作用推动的，则在复制事件之后，会立刻出现非同义替换速率大于同义替换速率的现象。但是，若这一基因家族在其余的时间都是受到 purifying selection 的作用，那么仅仅是两个序列之间的比对是很难发现 omega 值大于 1 的位点。鉴于这一问题，本文提出了用于检测基因复制后适应性进化的最大似然法，旨在完成以下目标：（1）同一个基因进化历史中不时间点的选择压力；（2）这些选择压力是否不同。1.1 支系间选择压力可变模型（Models of variable selective p
33、ressures among branches）1.1.1.模型介绍(1)one- ratio model ,即 Phylogenetic tree 中所有 sites 的 omega 值是一个恒定值。0=1=2=3；(2)Two- ratio model, 复制事件之前的 omega 值和复制事件之后的 omega 值不相等。01=2=3；90图 1 假设基因家族的进化树Fig1 Phylogeny of a hypothetical gene family(3)Three- ratio model, 该模型共假设三个 omega 值，即复制事件之前的 omega 值，复制事件之后的 ome
34、ga 值，以及后期分支中的 omega 值。01 2=3。(4)Four- ratio model, 该模型共假设四个 omega，不仅复制事件前后的 omega 值不同以外，随后的分支上的两个 omega 值也不相同。即 0123。1.1.2 模型比较可以揭示的问题（likelihood ratio test）(1) one ratio model with two ratio model : 揭示基因家族的复制前后所受的选择压力是否相同；(2) two ratio model with three ratio model: 揭示基因家族复制后到分化前以及分化后这段时间的选择压力是否发生变化
35、。(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
36、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 subcl
37、ade 与 EDN subclade 所处的进化选择压力不同。1.2 鉴定适应性进化的氨基酸位点（Identification of amino acid sites under positive selection）从蛋白质进化的角度而言，为了维持蛋白质功能的稳定性，其平均非同义替换的数量往往较少。只有少数一部分氨基酸位点在进化过程中受到适应性进化的选择。因此，计算整个支系的平均 omega 值往往很难检测到 positive selection。鉴于此，我们需要对小部分的 condons 所受的正选择压力进行检测并确定存在这些正选择压力的位点的位置。为完成上述目标，我们可以根据已知的结构域
38、和功能域的信息，将氨基酸位点分成若干具有独立 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 sitesModel A: 指定 0=0 且 1=1，因此正选择作用只约束在了前景枝中；通常选择 M1 与 Model A 进行 LRT 检测，自由度 df=2；Model B：0 和 1 为自由参数（free parameters） ，所分析的整个支系中均有可能出现正选择作用位点。通常选择 M3 与 Model B 进行 LRT 检测，自由度df=2；
