使用PAML软件利用密码子序列进行正选择分析

PAML软件中的codeml命令可以使用Maxmum Likelihood方法对密码子或氨基酸的多序列比对结果进行分析和进化参数计算 。这两种数据共用一个命令,程序运行时的参数绝大部分是公用的,但仍需要分清两种数据各自的参数意义。PAML软件中进行正选择分析的示例数据位于examples/lysozyme/目录下。

1.输入文件(1/3):含有拓扑结构的树文件

输入 一个树文件input.trees。该文件可以是有根树,也可以是无根树,也可以带有枝长信息。codeml命令用于分析各分枝上的omega值,若是有根树的话,对于分析root节点所在的分枝,其结果可能是不对的。因此,作者推荐使用无根树进行正选择分析。无根树是3分叉的树,最少需要3个物种;而3个物种的无根数有3个分枝(2n – 3),对其进行分析,仅能得到某个物种所对应分枝的omega值,这种情况下没太大意义;而使用4个物种的无根树进行分析,可以得到两个物种共同祖先分枝的omega值。故作者推荐使用至少4~5个物种进行分析。一个树文件内容示例:

7  1
((Hsa_Human, Hla_gibbon) #1,((Cgu/Can_colobus, Pne_langur) #1, Mmu_rhesus), (Ssc_squirrelM, Cja_marmoset));

文件内容分两行:第一行表述树中有5个物种,共计1个树,两个数值之间用任意长度的空分割。此外,Newick格式的树尾部一定要有分号,没有的话程序可能不能正常运行。

2. 输入文件(2/3):密码子序列比对结果文件 input.txt

7  390
Hsa_Human AAGGTCTTTGAAAGGTGTGAGTTGGCCAGAACT...
Hla_gibbon AAGGTCTTTGAAAGGTGTGAGTTGGCCAGAACT...
Cgu/Can_colobus AAGATCTTTGAAAGGTGTGAGTTGGCCAGAACT...
Pne_langur AAGATCTTTGAAAGGTGTGAGTTGGCCAGAACT...
Mmu_rhesus AAGATCTTTGAAAGGTGTGAGTTGGCCAGAACT...
Ssc_squirrelM AAGGTCTTCGAAAGGTGTGAGTTGGCCAGAACT...
Cja_marmoset AAGGTCTTTGAAAGGTGTGAGTTGGCCAGAACT...

PAML要求输入的Phylip格式,其物种名和后面的序列之间至少间隔两个空格(是为了允许物种名的属名和种名之间有一个空格)。

3. 输入文件(3/3):codeml命令的配置文件 codeml.ctl

软件难点在于配置文件的理解,常用示例如下:

*输入输出参数:
      seqfile = input.txt    *设置输入的多序列比对文件路径。
     treefile = input.trees  *设置输入的树文件路径。
      outfile = mlc          *设置输出文件路径。
        noisy = 9            *设置输出到屏幕上的信息量等级:0,1,2,3,9。
      verbose = 0            *设置输出到结果文件中的信息量等级:0,精简模式结果(推荐);1,输出详细信息,包含碱基序列;2,输出更多信息。
        getSE = 0            *设置是否计算并获得各参数的标准误:0,不需要;1,需要。
 RateAncestor = 0            *设置是否计算序列中每个位点的替换率:0,不需要;1,需要。当设置成1时,会在结果文件rates中给出各个位点的碱基替换率;同时也会进行祖先序列的重构,在结果文件rst中有体现。

*数据使用说明参数:
      runmode = 0      *设置程序获取进化树拓扑结构的方法:0,直接从输入文件中获取进化树拓扑结构(推荐);1,程序以输入文件中的多分叉树作为起始树,并使用星状分解法搜索最佳树;2,程序直接以星状树作为起始树,并使用星状分解法搜索最佳树;3,使用逐步添加法搜索最佳树;4,使用简约法获取起始树,再进行邻近分枝交换寻找最优树;5,以输入文件中的树作为起始树,再进行邻近分枝交换寻找最优树;-2,对密码子序列进行两两比较并使用ML方法计算DnDs,或对蛋白序列进行两两比较计算ML距离,而不进行其它参数(枝长和omega等)的计算。
* fix_blength = 1      *设置程序处理输入树文件中枝长数据的方法:0,忽略输入树文件中的枝长信息;-1,使用一个随机起始进行进行计算;1,以输入的枝长信息作为初始值进行ML迭代分析,此时需要注意输入的枝长信息是碱基替换率,而不是碱基替换数;2,不需要使用ML迭代计算枝长,直接使用输入的枝长信息,需要注意,若枝长信息和序列信息不吻合可能导致程序崩溃;3,让ML计算出的枝长和输入的枝长呈正比。
      seqtype = 1      *设置输入的多序列比对数据的类型:1,密码子数据;2,氨基酸数据;3,输入数据虽然为密码子序列,但先转换为氨基酸序列后再进行分析。
    CodonFreq = 2      *设置密码子频率的计算算法:0,表示除三种终止密码子频率为0,其余密码值频率全部为1/61;1,程序分别计算三个密码子位点的四种碱基频率,再算四种碱基频率在三个位点上的算术平均值,计算61种密码子频率时,使用该均值进行计算,这61种密码子频率之和不等于1,然后对其数据标准化,使其和为1,得到所有密码子的频率;2,程序分别计算三个密码子位点的四种碱基频率,计算61种密码子频率时,使用相应位点的碱基频率进行计算,这61种密码子频率之和不等于1,然后对其数据标准化,使其和为1,得到所有密码子的频率;3,直接使用观测到的各密码子的总的频数/所有密码值的总数,得到所有密码子的频率。一般选择第三种方法,设置CodonFre的值为2。第一种方法让所有密码子频率相等,不符合实际轻卡;第二种方法没有考虑到三种位点各碱基频率的差异,得到的密码子频率不准确;第三种方法能比较准确计算出各密码子的频率;第四种方法由输入数据得到真实的各密码子频率,但有时候有些非终止密码子的频率为0,可能不利于后续的计算。
    cleandata = 1      *设置是否移除不明确字符(N、?、W、R和Y等)或含以后gap的列后再进行数据分析:0,不需要,但在序列两两比较的时候,还是会去除后进行比较;1,需要。
*       ndata = 1      *设置输入的多序列比对的数据个数。
        clock = 0      *设置进化树中各分支上的变异速率是否一致,服从分子钟理论:0,变异速率不一致,不服从分子钟理论(当输入数据中物种相差较远时,各分支变异速率不一致);1,所有分枝具有相同的变异速率,全局上服从分子钟理论;2,进化树局部符和分子钟理论,程序认为除了指定的分支具有不同的进化速率,其它分枝都具有相同的进化速率,这要求输入的进化树信息中使用#来指定分支;3,对多基因数据进行联合分析。
        Mgene = 0      *设置是否有多个基因的多序列比对信息输入,以及多各基因之间的参数是否一致:0,输入的多序列比对文件中仅包含一个基因时,或多个基因具有相同的Kappa和Pi参数;1,输入文件包含多个基因,这些基因之间是相互独立的(这些基因之间具有不同的Kappa和Pi值,且其进化树的枝长也不相关);2,输入文件包含多个基因,这些基因具有相同的Kappa值,不同的Pi值;3,输入文件包含多个基因,这些基因具有相同的Pi值,不同的Kappa值;4,输入文件包含多个基因,这些基因具有不同的Kappa值和不同的Pi值。当值是2、3或4时,多个基因的进化树枝长虽然长度不一样,但是呈正比关系。这点和参数值等于1是不一样的。
        icode = 0      *设置遗传密码。其值1-10和NCBI的1-11遗传密码规则对应:0,表示通用的遗传密码。
   Small_Diff = .5e-6  *设置一个很小的值,一般位于1e-8到1e-5之间。推荐检测设置不同的值在比较其结果,需要该参数值对结果没什么影响。

*位点替换模型参数:
        model = 1            *若输入数据是密码子序列,该参数用于设置branch models,即进化树各分枝的omega值的分布:0,进化树上所有分枝的omega值一致;1,对每个分枝单独进行omega计算;2,设置多类omega值,根据树文件中对分枝的编号信息来确定类别,具有相同编号的分枝具有相同的omega值,没有编号的分枝具有相同的omega值,程序分别计算各编号和没有编号的omega值。
                             *若输入数据是蛋白序列,或数据是密码子序列且seqtype值是3时,该参数用于设置氨基酸替换模型:0,Poisson;1,氨基酸替换率和氨基酸的观测频率成正比;2,从aaRatefile参数指定的文件路径中读取氨基酸替换率信息,这些信息是根据经验获得,并记录到后缀为.dat的配置文件中。这些经验模型(Empirical Models)文件位于PAML软件安装目录中的dat目录下,例如,Dayhoff(dayhoff.dat)、WAG(wag.dat)、LG(lg.dat)、mtMAN(mtman.dat)和mtREV24(mtREV24.dat)等;
   aaRatefile = dat/wag.dat  *当对蛋白数据进行分析,且model = 2时,该参数生效,用于设置氨基酸替换模型。
       aaDist = 0            *设置氨基酸之间的距离。
      NSsites = 0            *输入数据时密码子序列时生效,用于设置site model,即序列各位点的omega值的分布:0,所有位点具有相同的omega值;1,各位点上的omega值小于1或等于1(服从中性进化neutral);2,各位点上的omega值小于1、等于1或大于1(选择性进化selection);3,discrete;4,freq;5:gamma;6,2gamma;7,beta;8,beta&w;9,betaγ10,beta&gamma+1;11,beta&normal>1;12,0&2normal>1;13,3normal>0。
                             *可以一次输入多个模型进行计算并比较,其结果输出的rst文件中。
    fix_alpha = 1            *序列中不同的位点具有不同的碱基替换率,服从discrete-gamma分布,该模型通过alpha(shape参数)和ncatG参数控制。该参数设置是否给定一个alpha值:0,使用ML方法对alpha值进行计算;1,使用下一个参数设置一个固定的alpha值。
	                     *对于密码子序列,当NSsites参数值不为0或model不为0时,推荐设置fix_alpha = 1且alpha = 0,即不设置alpha值,认为位点间的变异速率一致,否则程序报错。若设置了alpha值,则程序认为不同密码子位点的变异速率不均匀,且同时所有位点的omega值一致,当然各分枝的omega值也会一致,这时要求NSsites和model参数值都设置为0(这一般不是我们需要的分析,它不能进行正选择分析了)。
        alpha = 0            *设置一个固定的alpha值或初始alpha值(推荐设置为0.5)。该值小于1,表示只有少数热点位置的替换率较高;该值越小,表示位点替换率在各位点上越不均匀;若设置fix_alpha = 1且alpha = 0则表示所有位点的替换率是恒定一致的。
       Malpha = 0            *当输入的多序列比对结果中有多基因时,设置这些基因间的alpha值是否相等:0,分别对每个基因单独计算alpha值;1,所有基因的alpha值保持一致。
        ncatG = 5            *序列中不同位点的变异速率服从GAMMA分布的,ncatG是其一个参数,一般设置为5,4,8或10,且序列条数越多,该值设置越大。
                             *对于密码子序列,当NSites设置为3时,ncatG设置为3;当NSites设置为4时,ncatG设置为5;当NSites值设置>=5时,ncatG值设置为10。
    fix_kappa = 0            *设置是否给定一个Kappa值:0,通过ML迭代来估算Kappa值;1,使用下一个参数设置一个固定的Kappa值。
        kappa = 2            *设置一个固定的Kappa值,或一个初始的Kappa值。
    fix_omega = 0            *设置是否给定一个omega值:0,通过ML迭代来估算Kappa值;1,使用下一个参数设置一个固定的omega值。 
        omega = .4           *设置一个固定的omega值,或一个初始的omega值。
       method = 0            *设置评估枝长的ML迭代算法:0,使用PAML的老算法同时计算所有枝长,在clock = 0下有效;1,PAML新加入的算法,一次对一个枝长进行计算,该算法仅在clock参数值为1,2或3下工作。

4. 运行codeml进行分析

运行程序的命令:

codeml codeml.ctl

5. 正选择基因分析

对基因进行正选择分析,可以按如下步骤进行:

(1) 收集该基因在多个物种中的CDS和protein序列(可以使用orthoMCL分析结果)。
(2) 对该基因的蛋白序列进行多序列比对,再根据CDS序列转换成Codon序列比对结果。这个步骤需要自己编写一些程序来实现。
(3) 根据Codon序列的比对结果构建系统发育树。推荐使用RAxML软件使用ML算法对所有的单拷贝同源基因进行物种树构建。然后,输入物种树信息、Codon多序列比对信息,使用codeml进行omega计算和选择压模型检验。
(4) 快速过滤非正选择基因:设置runmode = -2运行codeml命令对两两序列比较,使用ML方法计算dNdS,若存在omega值大于1,则认为该基因可能属于正选择基因,进行后续分析。
(5) 正选择基因的初步鉴定方法:设置model = 0、NSsites = 1 2运行codeml命令(site models)分别对正选择模型(M2a)和近中性进化模型(M1a)进行检测,若两种模型的似然值相差较大,通过自由度为2的卡方检验,可以确定该基因是否为正选择基因。从程序结果中(例如:lnL(ntime: 11 np: 16): -870.867266 +0.000000)可以找到M2a模型的似然值l1和M1a模型的似然值l0,再计算2Δl = 2(l1-l0),通过命令“chi2 2 2Δl”根据卡方检验算出p值。此外,也可以设置NSsites = 7 8进行计算,则表示分别对M8(beta & omega)和M7(beta)模型进行检测,和前者类似,若两种模型的似然值相差较大,通过自由度为2的卡方检验,可以确定该基因是否为正选择基因。根据PAML作者所说,M2a和M1a的比较,比M7-M8的比较更严格。即若想得到更多的正选择基因,可以使用M7_VS_M8分析。此外,根据M2和M8模型的BEB(Bayes empirical Bayes)分析结果,还可以得到在多序列比对结果中第一条序列上的正选择位点。根据PAML作者所说,这种site models之间的比较,能用于检测是否在序列上不同的位点具有不同的omega值,而不是用于正选择检测(We suggest that The M0-M3 comparison should be used as a test of variable ω among sites rather than a test of positive selection)。
(6) 指定分化枝上的正选择基因鉴定方法:经过上一步初步鉴定后,设置model = 2、NSsites = 2、fix_omega = 0、omega = 2.0运行codeml命令(branch-site model A);再设置model = 2、NSsites = 2、fix_omega = 1、omega = 1运行codeml命令(modified branch-site model A / null model)。进行这两种模型分析时,要求输入的树文件中对目标分化枝进行标注。对这两种模型进行LRT分析,计算2Δl = 2(l1-l0),注意是前者的似然值减后者(null model)的似然值;再使用自由度为1的卡方检验,通过命令“chi2 1 2Δl”计算出的值再除以2,即得到p值。
(7) 目标分化枝上的快速进化基因(Rapidly evolving gene)鉴定方法:设置model = 0、NSsites = 0运行codeml命令,再设置model = 2、NSsites = 0运行codeml命令(这时要求输入的树文件中对目标分化枝进行标注)。然后比较两种运行模式下的似然值,通过自由度为1的卡方检验,可以鉴定该基因在目标分化枝的omega值和背景差异显著。

1 thought on “使用PAML软件利用密码子序列进行正选择分析

发表评论

电子邮件地址不会被公开。 必填项已用*标注

此站点使用Akismet来减少垃圾评论。了解我们如何处理您的评论数据