分类: PAML, 分子进化

「PAML」分子进化程序包PAML使用入门1-codeml

原文发布日期:20120929

前阵子为了算一些基因的Selection可算费了老劲,我的课题其实非常偏进化但是这方面的知识太欠缺,恶补了一阵子还是觉得很多东西一知半解。随手拿MEGA建个分子进化树非常容易,几分钟就能完成,但是如何选择模型?进化树是否具有生物学意义?如何解读一颗进化树?进化树的拓扑结构是否可重复?比如之前我一直找不到一个最最基本问题的答案:一颗进化树下面的Bar以及其数字代表什么?是否可以简单理解为序列的差异?后来在Nei的那本经典教材「Molecular Evolution and Phylogenetics」里面找到了部分答案,至少在距离模型里这个Bar的长度代表了序列的分歧度。

PAML程序包是由一个分子进化领域的大牛杨子恒教授开发的,他现在在伦敦大学,是英国皇家学会的会员,在计算分子进化领域的工作非常厉害。不过要是没点统计基础,想了解这个领域很多模型的算法非常难,我想对于大多数人来说,能够利用好PAML这个工具并能够解释结果就够了。这次就来简单介绍下PAML程序包中CODEML程序的使用。

感谢Boss Kuang的指导,小雨对于用CHi2计算LRT时候的指导,以及小麦童鞋在学习PAML遇到瓶颈时候的帮助。

「预备知识」

①、命令行的基本操作。

②、MEGA,Clustal, Treeview等程序的使用。

接下来进入正题,从最常用的codeml程序开始说起

①去官网下载PAML程序包:http://abacus.gene.ucl.ac.uk/software/paml.html, Windows版本应该不用解析,Unix/ Linux and Mac OSX请参考官网步骤解析程序。

②codeml程序运行需要以下4个文件,最好将他们放在一个目录下:codeml程序,codemlde.ctl配置文件,比对好的序列文件(Phlip或者PAML格式),以及tree文件。

③ 准备比对好的序列文件,序列要满足以下几个要求:

A.Phylips格式或者PAML格式

B.序列数需要时3的整数倍(密码子)

C.在序列中不能出现终止密码(最后一个如果是终止密码请手动删除)

具体制作步骤:把要进行Selection分析的序列保存成Fasta格式,用MEGA打开,用Clustal或者Muscle进行比对,Coden alignment以后手工删掉所有终止密码(Gap可以选择删除或者不删除,反正codeml计算的时候是忽略Gaps的)—从MEGA导出为比对好的Fasta文件。

④将比对好的Fasta文件转化成Phylips格式或者PAML格式,这种格式其实很简单,第一行是序列号和序列长度,下面就是序列名跟这序列把Fasta格式序列名前面的“>”号,Phylips和PAML格式的区别也仅仅在Phylips的序列名不能超过10个字母,而PAML的序列名没有限制。从这种格式可以看出来,手动将比对好的Fasta文件改成Phlips或者PAML格式是最简单的方法;替代方案是使用一个在线工具,提交比对好的Fasta文件就可以得到转化好的PAML文件:http://app3.titan.uio.no/biotools/tool.php?app=fasta-paml。最后说一下,Mac用户最好拿Windows完成步骤③④,因为在Mac下制作出来的序列文件经常会出问题,Linux没有测试过。PAML运行不了的原因几乎都是因为序列格式有问题。

最后,将转化好的序列格式保存成自己想要的名字就好,这里我存的是n1.txt

⑤制作tree文件,这里的tree文件最简单的方法就是直接用Clustal制作,推荐一个在线的Clustal工具:http://www.genome.jp/tools/clustalw/, 将序列提交上去,点Execute Multiple Alignment,然后将结果里的dnd文件下载,改名为自己想要的文件+后缀名“.trees”,我这里保存的名字为n1.trees。

⑥接下来要配置codeml.ctl文件,具体配置选项的功能和配置可以参考这个链接: http://csbl.bmb.uga.edu/yinyb/codeml.html,比较重要的几个有

seqfile = n1.txt 自己制作的比对好的序列名称

treefile = n1.trees 树文件名称

outfile = mcl 结果文件名称

runmode = 0 0: user tree; 1: semi-automatic; 2: automatic

seqtype = 1 1:codons; 2:AAs; 3:codons–>AAs

model = 0

models for codons:

0:one, 1:b, 2:2 or more dN/dS ratios for branches

models for AAs or codon-translated AAs:

0:poisson, 1:proportional, 2:Empirical, 3:Empirical+F

6:FromCodon, 7:AAClasses, 8:REVaa 0, 9:REVaa(nr=189)

NSsites = 0 1 2 7 8 0:one w;1:neutral;2:selection; 3:discrete;4:freqs;

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

其中NSsites就是模型选择,我一般都是选择0,1,2,7,8,这样在接下去可以对结果进行LRT分析来看结果统计上的显著性。

⑦保证这四个文件:codeml, codeml.ctl, n1.txt, n1.trees在一个文件夹,cd到目标目录,然后运行codeml程序。等命令行跳半天以后目录里就会出现很多个新的文件,其中mcl就是主要文件的结果。如果数据量大的话可能会很慢,请考虑选择服务器进行运算。

⑧目标目录会出现很多结果文件,其中MCL里面是主要的结果。

⑨打开mcl文件,看M8模型下面的这个部分,第一列是序列的位置和承受选择的氨基酸,第二列是W>1(W>1: positive selection, W<1: purifying selection)的概率,如果P>95%就会被标记上,证明该位点收到positive selection。

这样的结果其实还不可信,接下来还应该用通过比较M7,M8, 计算LRT来看结果统计上的可信度,至于如何算嘛,请听下回分解。

「参考资料」

①、「计算分子进化」杨子恒

②、「分子进化与系统发育」Masatoshi Nei

③、paml中文手册:http://blog.csdn.net/winglyx/article/details/6731787

④、codeml配置文件说明:http://csbl.bmb.uga.edu/yinyb/codeml.html

⑤、codeml使用流程介绍:http://fhqdddddd.blog.163.com/blog/static/1869915420111091276909/