..微笑着看着杯中的花茶一片片撑开.. ..透明的花瓣里水破开的声音很轻微..
分类: 信息化
2016-07-12 16:19:10
By 斑斑 <23920620>
PAML(将最大似然法用于系统发育分析)是一个用最大似然法来对DNA和蛋白质序列进行系统发育分析的软件包。
除了这个手册,请注意以下资源:
PAML站点:http://abacus.gene.ucl.ac.uk/software/PAML.html 含有下载和编译程序的信息。
PAML FAQ页面:http://abacus.gene.ucl.ac.uk/software/pamlFAQs.pdf
PAML讨论组 http://www.rannala.org/phpBB2/,在这里你可以提交bug报告和提问。
PAML软件包目前包含如下程序:baseml,basemlg,codeml,evolver,pamp,yn00,mcmctree,和chi2。Yang(2007)提供了一个简短的关于PAML中常用的模型和方法的综述。有本书(Yang 2006)描述了统计和计算的详细内容。分析的例子可以从软件包中找到。
l 比较和测试系统发育树(baseml 和codeml);
l 复杂替代模型中的参数估算,包括在位点和模型中用于对多个基因或位点划分进行结合分析的变化率模型。(baseml 和 codeml)
l 通过对实施的模型的对照进行似然率假设性检验(baseml,codeml,chi2);
l 基于全局或局部分子钟模型估算分歧时间(baseml 和 codeml);
l 使用核酸,氨基酸和密码子模型对祖先序列似然性(经验贝叶斯)重构(baseml 和 codeml);
l 用Monte Carlo模拟生成核酸,密码子和氨基酸序列数据集(evolver);
l 估算同义和非同义替代率和发现蛋白编码DNA序列中的正选择;
l 贝叶斯估算物种的分歧时间来合并化石尺度上的不确定性。
PAML的长处是它的一系列复杂替代模型。Baseml和codeml中使用的树搜索算法还都很简单,因此除非是非常小的数据如小于10个物种,否则你最好使用其他的软件包,如phylip,paup,或mrbayes来推断树的拓扑结构。你可以使用其他程序获得一系列的树并那它们作为使用树来通过baseml和codeml来评估它们。
Baseml和codeml:Baseml这个程序是使用最大似然法分析核酸序列。Codeml这个程序由两个旧程序合并构成,它们是codonml(Goldman和Yang(1994)的对蛋白编码DNA序列执行密码子替代模型)和aaml(用来执行氨基酸序列模型)。它们俩现在靠控制文件codeml.ctl中的变量seqtype来区分,1对应密码子序列,2对应氨基酸序列。这个文档中,我们使用codonml和aaml意味着codeml中的seqtype分别等于1和2。程序baseml,codeml和aaml使用相似的算法(最大似然法)来配合模型。三个程序主要的不同是马尔科夫模型中的进化单元不同,序列中被称为“位点”的分别是一个核酸,一个密码子或者一个氨基酸。马尔科夫模型用替代率(计算位点间不变或变化)来描述核酸,密码子或氨基酸间的替代。
Evolver:这个程序基于核酸,密码子和氨基酸替代模型模拟序列。它还具备一些其它的选项如生成随机树,计算树之间的分隔距离(Robinson 和 Foulds 1981)。
Basemlg:这个程序执行(连续的)Yang(1993)的gamma模型。当数据超过6或7个物种时它得到速度非常慢并且难执行。作为替代应使用Baseml中的离散gamma模型。
Mcmctree:这个执行的是Yang和Rannala(2006)和Rannala和Yang(2007)的bayesian MCMC算法来估算物种的分歧时间。
Pamp:这个执行的是Yang和Kumar(1996)的基于简约的分析。
Yn00:这个执行的是Yang和Nielsen(2000)用来成对估算蛋白编码DNA序列中同义和非同义替代率(ds 和 dn)的方法。
Chi2:这个用来执行似然率测试。它计算chi平方临界值。你可以将它同真实数据的检验统计相比较来确定检验在5%和1%的水平上是否显著。键入程序名chi2来运行这个程序。当你输入检验统计值和d.f时这个程序也能计算P值。通过键入chi2 p运行。
有些事情你可能期望一个系统发育软件包应该做,但PAML不能做。这儿有部分列表,希望能帮你避免浪费时间。
序列对齐:你应该使用其它的一些程序,如Clustal或者TreeAlign来自动对齐序列或者手动对齐,或者从其它一些如BioEdit和GeneDoc这样的程序中得到帮助。手动调整似乎没有达到可以完全信赖的成熟的阶段,那么如果有可能,你应当总是手动调节。如果你在基因组范围内构建成千上万的对齐,你应该实行一些质量控制,计算一些序列分歧的评估来作为对齐的不可信度的指标。对于编码序列,你应当对齐蛋白序列并且基于蛋白的对齐来构建DNA的对齐。值得注意的是在baseml和codeml(如果cleandata=1)中,对齐的空格是作为缺失的数据来对待的。如果cleandata=1,所有的模糊位点和空格都会被移除。
基因预测:codeml中基于密码子的分析(codeml seqtype=1)假定序列是从外显子开始对齐,序列的长度是3的倍数,序列中的第一个核酸是编码的第一个位置。内含子、间隔和其它非编码区域必须被移除并且编码序列必须在运行程序前对齐好。这个程序不能直接运行从GenBank下载回来的序列,即使这里含有CDS的信息。它也不能预测出编码区域。
在大规模数据集中树搜索:如先前提到的,你应该使用其它的程序来获得一棵树或一些候选树,并使用这些树作为使用树来配合其它软件包中不可用的模型。
PAML使用旧式的命令行界面。你可以从PAML网站下载存档,典型的以PAML*.*tar.gz命名,在你硬盘上解压这些文件。这是一个针对所有平台的文件。可执行文件是针对windows,对于UNIX和MAC OS X,你需要在运行前先编译他们。
软件包中包含windows的可执行文件。
1.
到PAML的站点http://abacus.gene.ucl.ac.uk/software/paml.html下载最后的存档并保存在你的硬盘上。解压,如,使用WinZip,存档在一个文件夹里,如D:\software\paml\.记住文件夹名称。
2.
启动一个“命令行”,从“开始 - 程序 - 附件”。也可以“开始 - 运行”并键入cmd,点击ok。你可以右击标题栏来改变窗口的字体,颜色,大小等。
3.
改变路径到paml的文件夹。如下:
D:
d:\software\paml
dir
4.
注意,windows的命令和文件名是大小写不敏感的。文件夹 src 包含源文件。example各种例子文件。bin 文件夹包含windows的可执行文件。你可以使用windows浏览器来查看这些文件。在当前文件夹下,使用默认的控制文件baseml.ctl运行baseml程序。命令行大概如下。
D:\software\paml4\bin\baseml
Baseml将读取当前目录下的默认控制文件baseml.ctl并根据它的参数来执行分析。现在我们可以输出一个baseml.ctl的副本,并用一个文本编辑器来查看相关的序列和树文件。类似地,你也可以运行codeml并查看它的控制文件codeml.ctl.
接下来你可以准备自己序列文件和树文件。控制文件和其它输入文件都是纯文本文件。一个常见的错误是由于UNIX和windows对回车和换行的不同处理方式造成的。如果你使用MS word来准备这些输入文件,你应该保存为“以换行结尾的文本”或“不带换行的文本”。有时只有两者中的一种有用。不要将文件保存成word文档。我在附录B的“克服windows的困扰”一节中收集了一些注意事项。如果你坚持双击,你可以打开windows浏览器,复制可执行文件到你包含控制文件的文件夹,然后双击运行。
软件包里没有提供UNIX的可执行文件,因此你必须使用软件包src文件夹中的源文件自己编译它们。值得注意的是UNIX命令行和文件或文件夹名称是大小写敏感的。下面是假设你在UNIX提示符下。
1.
到PAML站点http://abacus.gene.ucl.ac.uk/software/paml.html下载最后的文件并保存在硬盘上。使用gzip解压。命令如下:
gzip –d paml4.tar.gz xf paml4.tar
2.
你可以使用 ls 查看文件夹中的文件。删除bin文件夹下的windows可执行文件(.exe文件)。然后切换到src文件夹下用make编译。
ls -lF bin (列出bin文件夹下的.exe文件)
rm –r bin/*.exe
cd src
make
ls -lF
rm *.o
mv baseml basemlg codeml pamp evolver yn00 chi2 ../bin
cd .. bin/codeml
4.
这些命令编译程序并生成可执行文件baseml,basemlg,codeml,pamp,evolver,yn00,和 chi2。你可以用ls命令查看到。然后删除中间的目标文件*.o,然后移动编译好的可执行文件到PAML的主文件夹bin文件夹。用cd命令切换到PAML主文件夹并使用默认的控制文件codeml.ctl运行codeml.你可以输出codeml.ctl的副本并查看它(主要的结果文件mlc)。
如果编译(make命令)文件不成功,你可能必须在使用make命令前打开并编辑makefile文件。例如,你可以将cc改变为gcc 和-fast到-O3或-O4。如果这些都不起作用,查看scr文件夹中的readme.txt文件中的编译指令。你可以复制编译命令到命令行。例如:
cc –o baseml baseml.c tools.c –lm
cc –o codeml codeml.c tools.c -lm
将使用C编译器cc编译 baseml 和 codeml。然而这种情况下代码的优化没有被开启。你需要使用编译器切换到代码优化,如,
cc –o codeml –O3 codeml.c tools.c -lm
最后,如果你当前的文件夹不在你的搜索路径,你必须添加./到可执行文件的名字前面甚至可执行文件在你的当前工作目录时也是如此。即用./codeml 代替 codeml 来执行codeml。
Mac OS X 如同 UNIX,你需要按照上面UNIX的说明来做。打开一个命令终端(应用 – 工具 – 终端)并且从终端编译和运行程序。你 cd 到 paml/src 文件夹并查看 readme.txt 或 makefile 文件。参照上面。如果你键入命令 gcc 或 make 并得到“命令未找到”的错误,你必须从Apple的站点http://developer.apple.com/tools/下载Apple Developer’s Toolkit。在FAQ里有一些MAC OS X或 UNIX下运行程序所要注意的内容。
我已停止发布MAC OS 9以前的老版本执行程序。
如同上面指出的,你可以通过命令行键入它的名字来运行一个程序。你需要知道你的序列文件,树文件和控制文件在哪个文件夹,这与你的工作文件夹有关。如果不熟练,你可以拷贝这些可执行文件到包含你数据的文件夹。根据模型的使用,codeml可能需要一些数据文件诸如 grantham.dat,dayhoff.dat,jones.dat,wag.dat,mtREV24.dat,或 mtmam.dat。因此你可能也要复制这些文件。
这个程序产生一些名字如rub,lnf,rst,或 rates的结果文件。不要将这些名字用于自己的文件,否则他们将被覆盖。
Examples文件夹包含许多例子数据。它们被作为原始文件来测试新的方法,我将他们包含进来以便你可以重复出我们文档中的结果。对齐的序列,控制文件和详细的说明文件也包含在内。它们可以帮助你获得熟悉的输入数据格式和带有解释的结果,也能帮助你发现程序中的bug。如果你对个例的分析感兴趣,那么拷贝一个带有描述方法和分析例子的文件来重复已公布的结果。这非常重要,因为这个手册只写出和描述了程序控制文件中变量的意义,但没有清楚的解释如何针对个例分析设定控制文件。
examples/HIVNSsites/: 这个文件夹包含Yang等人(2000b)的HIV-1 env V3 区域分析的例子文件。这个数据用来演示文档中描述的NSsites模型,也就是,氨基酸位点中变量ω比值的模型。这些模型被Yang和Swanson(2002)称为 “random-sites”模型。这是由于我们在推理中不知道哪些位点是高度保守并承受正选择的。它们也称为“fishing-expedition”模型。这些数据集是Yang等人(2000b)中分析的第十组数据,结果在该文章的表格12中。参看文件夹中的说明文件。
examples/lysin/:这个文件夹包含的是Yang, Swanson & Vacquier (2000a) 和 Yang and Swanson (2002)对25种鲍鱼的精子细胞溶解酶基因的分析数据。这个数据集是为了演示“random-sites” (Yang, Swanson & Vacquier 2000a)模型和“fixed-sites”(Yang and Swanson 2002)模型的。
后面的这篇文章中,我们使用了结构信息来区分细胞溶解酶中“埋入”的和“暴露”的氨基酸位点信息。并指定和估算两个不同分区的ω比值。假设就是暴露在表面的位点有可能承受正选择。参见文件夹中的说明文档。
examples/lysozyme/:这个文件夹包含的是Messier 和 Stewart(1997)及Yang(1998)重分析的灵长类溶菌酶c基因的数据。这是为演示codon模型,它能够对树的不同枝指定不同的ω比值,对沿着种系进行的正选择测试非常有效。这个模型有时也叫做branch模型或branch-specific模型。Yang(1998)的大小数据集都包含在内。这些模型需要用户标出树种的枝,说明文件和所包含的树文件里对格式进行了非常详细的解释。也可以参看“树文件和树拓扑结构的表示方法”一节中关于指定枝和节点的部分。
Yang 和Nielsen(2002)也使用了溶菌酶数据来执行被称为“branch-sites”的模型,它允许ω比值在种系和位点中都可以变化。参见学习如何运行这个模型的说明文件。
examples/MouseLemurs/: 这个文件夹包含的是Yang和Yoder(2003)通过 mtDNA比对来分析和估算老鼠与狐猴的分歧时间的数据。这个数据集用来演示如何通过最大似然法按照全局和局部分子钟估算分歧时间。文档中描述了更复杂的同时使用多个校准节点的模型的使用方法,对多个基因(或位点分区)进行分析来解释各组分枝间的不同和变化的比值。说明文件中解释了输入数据的格式和模型的详细设定方案。说明文件2中解释了Yang(2004)中对特别比值的平滑过程。
examples/mtCDNA/: 这个文件夹包含的是线粒体基因组中12个编码蛋白基因的同侧序列的比对结果。这些是Yang, Nielsen和Hasegawa(1998)用于对7种猿的密码子和氨基酸替代模型的分析数据。这个是文章中提到的小数据集,用来配合氨基酸替代率的mechanistic和empirical模型,也可以用于密码子替代的mechanistic模型。这个模型可用于测试保守的和进化的氨基酸替代率是否相等。详细的参看说明文件。
examples/TipDate/:这个文件夹包含的是Rambaut(2000)在TipDate模型中使用的例子数据。用于时间已知的病毒序列。说明文件解释了如何使用baseml来配合TipDate模型,用在一个不同时间的序列的全局的时钟上。局部时钟模型也可以应用。如何做请参看examples/MouseLemurs/文件夹。注意,我们使用符号@在序列名称前作为序列测定数据的前缀。文件可以使用Ranmbaut的TipDate程序阅读,但是我们软件包中的文件在送入baseml前需要进行一些编辑。
软件包中也包含了一些其它的数据,下面是详细内容。
brown.nuc and brown.trees: Brown 等人(1982)的895-bp的mtDNA数据。它们被Yang(1994)等人和Yang(1994b)用来测试位点中变化比值的模型。
mtprim9.nuc and 9s.trees:由来自9个灵长类的888个对齐的位点构成的线粒体片段。被Yang(1994a)用于测试离散gamma模型和Yang1995测试自动离散gamma模型。
abglobin.nuc and abglobin.trees:连锁的α和β球蛋白基因。被Goldman和Yang(1994)用在对密码子模型的描述中。abglobin.aa是翻译的氨基酸序列的比对数据。
stewart.aa and stewart.trees:6种哺乳动物的溶解酶蛋白序列(Stewart等人,1987)。被Yang等人(1995a)用来测试重建原始氨基酸序列的方法。
参看软件包中示例数据文件(.nuc .aa和.nex)。只要你获得其中一种格式的数据,PAML程序就可以读入它。默认的格式是(Felsenstein 2005)PHYLIP软件包中的PHYLIP格式。PAML有限的支持PAUP和MacClade中使用的NEXUS文件格式。只有序列数据或树可以被读入,注释块被忽略。PAML不处理序列数据块中的注释块,因此请清除它们。
下面是一个PHYLIP(Felsenstein 2005)格式的例子。第一行包含物种的数目和序列长度(可能跟随着选项符)。对于密码子序列(codeml seqtype = 1),序列文件中的序列长度是指核酸的数目而不是密码子的数目。序列文件中的选项只能是I,S,P,C和G。序列既可能是交替格式(选项I,示例数据文件abglobin.nuc),或者连续的格式(选项S,示例数据文件brown.nuc)。默认选项是S,因此你并不需要指定它。选项G用于多基因数据的结合分析,见下面解释。这面时一个连续型示例数据集。它有4条60个核酸(或20个密码子)的序列。
4 60
sequence 1
AAGCTTCACCGGCGCAGTCATTCTCATAAT
CGCCCACGGACTTACATCCTCATTACTATT
sequence 2
AAGCTTCACCGGCGCAATTATCCTCATAAT
CGCCCACGGACTTACATCCTCATTATTATT
sequence 3
AAGCTTCACCGGCGCAGTTGTTCTTATAAT
TGCCCACGGACTTACATCATCATTATTATT
sequence 4
AAGCTTCACCGGCGCAACCACCCTCATGAT
TGCCCATGGACTCACATCCTCCCTACTGTT
物种/序列名称:物种和序列名字中不要使用下面的特殊符号¨, : # ( ) $ =〃,它们用于特殊的目的可能会引起程序的混乱。@符号可以作为序列名称的一部分或结尾来指定序列测试的时间,例如virus1@1984。@符号被认为是名称的一部分表示这个序列是在1984年测试的。物种名称的最大字符数在主程序baseml.c和codeml.c的开头指定的。在PHYLIP中,物种名是10个字符,我发现经常受限制,因此我将默认值设为30。为解决这个问题,PAML中认定两个连续的空格作为物种名的结束符,以避免物种名要30个(或10个)字符。为遵守这一规则,你不能再一个物种名中使用两个连续的空格。例如上面的数据也可以用下面的格式。
4 60
sequence 1 AAGCTTCACCGGCGCAGTCATTCTCATAAT
CGCCCACGGACTTACATCCTCATTACTATT
sequence 2 AAGCTTCACCGGCGCAATTATCCTCATAAT
CGCCCACGGACTTACATCCTCATTATTATT
sequence 3 AAGCTTCACC GGCGCAGTTG TTCTTATAAT
TGCCCACGGACTTACATCATCATTATTATT
sequence 4 AAGCTTCACCGGCGCAACCACCCTCATGAT
TGCCCATGGACTCACATCCTCCCTACTGTT
如果你想使文件同时被PHYLIP和PAML识别,你需要限定名称的字符数为10并且至少用两个空格来分隔名字和序列。
在一个序列中,T,C,A,G,U,t,c,a,g,u都被识别为核酸序列(在baseml,basemlg和codonml中),标准的单字符编码(A, R, N, D, C, Q, E, G, H, I, L, K, M, F, P, S, T, W, Y, V及小写形式)被识别为氨基。模糊字符(未定的核酸或氨基酸)也允许使用。三个特殊字符“.”“-”“?”做如下解释:圆点表示在第一个序列中相同的字符,破折号表示一个比对的空位,问号表示未确定的核酸或氨基酸。序列中的非字母符号,诸如>