PAML 中文说明.docx
《PAML 中文说明.docx》由会员分享，可在线阅读，更多相关《PAML 中文说明.docx（24页珍藏版）》请在冰豆网上搜索。
PAML中文说明
PAML:
最大似然法分析系统发育
PhylogeneticAnalysisbyMaximumLikelyhood
版本：
4.3（2009年9月）
ZihengYang
马向辉翻译
1、概述
PAML（forPhylogeneticAnalysisbyMaximumLikelihood）是一个用最大似然法分析蛋白质或DNA序列系统发育的一个程序包。
1.1PAML文件：
除了这个手册以外，以下资源也需要注意：
PAML网站：
http:
//abacus.gene.ucl.ac.uk/software/PAML.html。
在这个网站上有PAML的下载以及编译程序;
PAMLFAQ页面：
http:
//abacus.gene.ucl.ac.uk/software/pamlFAQs.pdf;
PAML讨论群：
http:
//www.rannala.org/phpBB2/，在这里你可以提出你的问题，或者提出你发现的漏洞。
1.2PAML可以做些什么？
PAML的最新版本包含一下几个程序模块：
baseml,basemlg,codeml,evolver,pamp,yn00,mcmctree,以及chi2。
其中最常用的模块的介绍可以参考杨子恒教授2007年发表的文章。
模块运行中用到的计算、统计方法在杨子恒教授的书中有详细的介绍。
模块的主要作用包括：
计算以及检测系统发育树（baseml和codeml）;计算复杂的碱基替代或者氨基酸替代模型中的参数，如不同位点间不同速率的模型或多个基因或者位点的综合分析模型（baseml和codeml）;用似然比例检测比较几个模型（baseml，codeml以及chi2）;用全局分子钟或者局部分子钟估算分歧时间（baseml和codeml）;用最大似然法重建祖先氨基酸、核苷酸序列以及密码子模型（baseml和codeml）;用蒙特卡洛模拟生成氨基酸、密码子或者核酸序列（evolver）;估算同义替代、非同义替代的速率，检测DNA的蛋白编码区的正选择（yn00和codeml）;综合贝叶斯法以及化石校正估算物种分歧时间（mcmctree）。
PAML的优势在于它整合了各种复杂的替代模型。
在baseml和codeml中建树的算法相对简单，所以较少的物种（如<10个）可以用这两个软件分析，对于大量物种的建树分析，最好还是用其他的程序去分析树结构，例如phylip、paup或者myBayes。
当然，你可以用其他的软件构树，然后作为用户树用baseml或codeml验证。
baseml和codeml：
baseml程序用于最大似然法分析核苷酸序列;codeml程序则是由两个旧程序组合而成：
codonml和aaml。
其中前者是基于Goldman和Yang在1994年提出的编码蛋白质的核酸序列的密码子替代模型，而后者主要用于氨基酸序列的替代模型。
现在，这两种序列可以在codeml.ctl中通过seqtype定义，其中1表示密码子序列，2表示氨基酸序列。
在这个手册里面，我将使用codonml和aaml来分别表示codeml中的seqtype=1和seqtype=2。
这三个程序（baseml，codonml和aaml用相同的最大似然算法对于模型进行拟合，而三者之间主要的不同点在于，三个程序中对于序列进化的马尔科夫模型中“位点”的定义：
在baseml中一个位点表示一个核苷酸，在codonml中一个位点表示一个密码子，aaml中一个位点表示一个氨基酸。
马尔科夫过程模型常常用于描述核苷酸、氨基酸序列之间或者密码子之间的替代。
对于不同的位点，这种替代既可以是恒定的，也可以是可变的。
evolver：
这个程序可以在特定的核苷酸、氨基酸、密码子替代模型下模拟序列的产生。
它还可以用于其他的一些操作，如产生随机树、计算树间的距离。
basemlg：
这个程序主要用于执行Yang在1993年提出的gamma模型的运算。
在计算6或7个物种以上的数据时，这段程序运算非常的慢，而且较难执行。
而baseml程序中的不连续的gamma模型则可以弥补这一不足。
mcmctree：
这个程序用于计算物种的分歧时间，使用的模型是Yang和Rannala在2006和2007年提出的。
pamp：
这个程序用于执行Yang和Kumar在1996年提出的简约分析。
yn00：
这段程序用于计算蛋白质编码的DNA的同义突变和非同义突变的速率，运算主要基于2000年Yang和Nielsen提出的方法。
chi2：
用于似然比例运算。
它可以计算chi平方的临界值，这个值可以用于比较统计值和真实值之间的差异，据此可以检测在显著性水平为5%或者1%时的差异是否显著。
运行这个程序可以直接输入“chi2”并回车。
另外，这个程序还可以计算p值。
这时候需要输入“chi2p”并回车。
1.3PAML不能做什么？
你可能希望一个系统发育软件可以做很多事情，但是其中有许多是PAML不能完成的。
下面是一个不完全的列举，希望你不要因此而浪费时间。
1）序列排列：
你可以用一些程序自动的排列核苷酸或氨基酸序列（如Clustal或TreeAlign），也可以在其他一些软件（如BioEdit和GeneDoc）的帮助下手动排列序列。
自动排列序列还远没有达到成熟的地步，所以手动调整是必须的。
如果你要进行数千个基因组水平的排列，你必须进行排列质量的控制，例如计算序列间的距离，并一次作为序列排列是否可靠的标准。
对于编码序列，可以首先对蛋白序列进行排列，然后根据蛋白序列排列核苷酸序列。
需要注意的是，在baseml和codeml中序列排列时产生的gap会被识别为缺省数据（如果cleandata=1）。
如果cleandata设置为1，那么所有的不明确的数据以及gap数据都会被删除。
2）基因预测：
程序codonml只能用于编码序列的分析，所以codonml在运行事是假设序列是预先排列好的外显子，并且序列长度为3的倍数，序列中的第一个核苷酸会被识别为密码子位置1。
内含子、居间序列以及其他非编码区域必须事先删除，编码序列也必须事先排列完毕。
程序也不能处理那些直接从GeneBank里面下载的数据，就算CDS信息已经确定也不行。
这段程序不能用于编码区的预测。
3）在大量数据中寻找构树信息：
如前所述，你可以通过其他的软件得到一个树或者几个备选树，然后用他们作为用户树去拟合某个模型。
2、PAML程序的编译以及应用
PAML使用老式且简单的命令行界面。
你可以从PAML的网站上下载档案文件（一般来说，这个档案文件的名字时“PAML*.*.tar.gz”），然后解压缩到本地磁盘上。
这个文件适用于所有的操作系统平台，对于Windows用户来说，只需解压即可，对于UNIX或者MACOSX用户，则需要在运行程序之前编译一下。
2.1对于Windows用户
Windows的可执行文件放在如下所述的文件夹中：
1）到PAML的网站上（http:
//abacus.gene.ucl.ac.uk/software/paml.html）下载最近升级的档案文件并存放到本地硬盘上。
然后将档案文件解压缩（如用WinZip）到一个特定的文件夹中，例如（D:
\software\paml\），记住文件夹的名字。
2）进入命令行模式。
可以通过以下方式：
“开始-程序-附件”或者“开始-运行”然后在命令框中输入“cmd”然后点击确定。
你可以通过右键点击标题栏改变窗口的字体、颜色、大小等。
3）将目录改到PAML的文件夹下。
例如你可以输入：
d:
cd\software\pamldir
4）注意，Windows命令以及文件名不区分大小写。
“src”文件夹中存放的时一些源文件，“examples”文件夹中存放着各式各样的实例文件，“bin”文件夹中存放着Windows的可执行文件。
你可以用WindowsExplorer浏览这些文件。
如果需要在当前目录下通过默认的控制文件“baseml.ctl”运行baseml程序的话，你可以在命令行中输入这些：
bin\basemlD:
\software\paml4\bin\baseml
这可以让baseml程序在当前文件夹中读取默认的控制文件“baseml.ctl”，并根据其中设置的参数进行分析。
这时，你就可以输出“baseml.ctl”的拷贝，并用文本文件编辑器浏览相应的序列以及树文件了。
同样，对于codeml程序的控制文件“codeml.ctl”也可以使用相同的操作。
接下来你就可以准备你自己的序列数据文件以及树文件了。
控制文件以及其他输出文件都是一些简单的文本文件。
一个普遍存在的问题是由于UNIX和Windows对于回车符和换行符的使用和识别。
如果你使用MSWord准备输入文件的话，需要把这样存储文件“另存为-纯文本”，不要存储为Word文档。
我收集了一些注意事项，这些注意事项在附加文件B的“OvercomingWindowsAnnoyances”详细列出。
如果你坚持使用双击，你可以打开WindowsExplorer，然后拷贝所有的可执行文件，然后把它们粘贴到包含控制文件的文件夹中，双击可执行文件即可。
2.2对于UNIX和MacOSX用户（略去）
2.3运行程序
如上所述，你可以在命令行模式中输入程序名来运行一个程序，但是你应该知道，你的序列文件、树文件、控制文件相对于你当前的工作目录的位置。
如果不熟练的话，你可以把所有的可执行文件拷贝到数据文件夹中。
对于codeml程序，可能会需要一些数据文件，如：
grantham.dat，dayhoff.dat,jones.dat,wag.dat,mtREV24.dat或者mtmam.dat，所以你最好还是一起拷贝这些文件。
程序们运行的结果将会以特定的文件名存放，如rub、lnf、rst或者rates。
你不要用这些文件名命名自己的文件，因为这些文件会被覆盖的。
2.4数据格式实例
examples文件夹中包含许多数据的实例。
这些数据在出处的论文中是用于检测新方法的，我把这些数据文件放到这里是为了你们可以重现论文中的结果。
序列的排列、控制文件以及详细的readme文件也包含在内。
这些都是为了帮助你们发现程序中的bug。
如果你对于某个特定的分析及其结果感兴趣的话，你可以找到相应论文，用其中描述的方法分析实例数据以重现那些发表过的结果。
这时尤其重要的，因为手册中所述的只是程序中用到的各类变量及其意义，但是并未清楚地描述如何针对特定的数据设置控制文件中的参数。
1）examples\HIVNSsites\文件夹：
这个文件夹中包含了Yang在2000发表的论文中使用的HIV-1病毒的envV3区域的序列数据。
这些数据是为了阐述文章中的NSsites模型，即不同氨基酸位点的不同ω速率的模型。
这个模型在之后发表的文章中（Yang&Swanson）称为“random-sites”模型，因为有这样一个前提：
我们不清楚哪些位点可能是高度保守的，哪些位点是经历正选择的。
这个模型优势也称作“fishing-expedition”模型。
文件夹中的数据是2000年论文中的第十组数据，分析结果则列在那篇论文的TABLE12中。
请阅读此文件夹中的readme文件。
2）examples\lysin\文件夹：
这个文件夹中存放着25个鲍鱼物种的细胞溶解酶基因，这些基因在两个文章中分析报道（Yang,Swanson&Vacquier（2000a）andYangandSwanson（2002）），这些数据用了两种模型进行了分析：
“random-sites”模型（asinYang,Swanson&Vacquier（2000a））和“fixedsites”模型（asin（YangandSwanson2002））。
在2002年的论文中，我们根据结构信息把细胞溶解酶中的氨基酸分成了两种：
包埋的和暴露的，根据这个分别计算不同位置氨基酸的ω比例。
这种分析证实了一个假设：
位于蛋白表面的氨基酸是受正选择的。
可以参考文件夹中的readme文件。
3）examples\lysozyme\文件夹：
这个文件夹包含灵长类的溶菌酶c基因，这个酶基因在1997年曾被Messier和Stewart分析过，但是在1998年Yang对这些数据进行了重新分析。
这是为了阐述这样一个模型：
系统发育树上不同的枝有不同的ω比例。
这对于检测进化枝中的正选择是有用的。
这些模型有时被人们称作分支模型或者分支特异性模型。
在1998年Yang发表的论文中的大数据和小数据都包括在文件夹中。
这些模型需要用户为发育树上的分支分类，在readme文件中包含树文件并非常详细解释了树的结构。
请参考随后的“Treefileandrepresentationsoftreetopology”。
溶菌酶数据还在Yang和Nielse的2002年的论文中分析“branch-site”模型，这个模型允许ω在进化枝和氨基酸位点之间存在差异。
请参考readme文件中关于这个模型的使用。
4）examples\MouseLemurs\文件夹：
这个文件夹中包括mtDNA的排列后数据，这些数据在YangandYoder在2003年发表的论文中用于估计mouselemurs的分歧时间。
这些数据用于阐述在全局和局部分子钟模型下，用最大似然法估计物种的分歧时间。
在这个文章中描述的最复杂的模型是同时应用了多重的进化速率（multiplecalibrationnodes）计算多个基因之间的不同，此外还计算了不同进化枝之间的不同速率。
readme文件解释了2004年Yang论文中的adhoc比例平滑程序。
5）examples\mtCDNA\文件夹：
这个文件夹中存放了Yang、Nielsen和Hasegawa在1998年发表的论文中使用的数据，这些数据包括猿类动物线粒体DNA编码的12个蛋白编码基因，对于这些序列的分析基于数种不同的氨基酸和密码子替代模型。
根据这个论文的说法，这些数据是小数据（“small”dataset），分析这些数据不仅仅采用密码子替代的机理模型（“mechanistic”model），还使用了经验模型（empiricalmodel）。
这个模型可以用于检测保守的和突变率较高的氨基酸位点的替代速率是否相等。
详细叙述请参考文件夹中的readme文件。
6）examples\TipDate\文件夹：
这个文件夹中包括了Rambaut在2000年发表的论文中的数据，这些数据用于描述他的TipDate模型。
readme文件中阐述了如何用baseml程序拟合TipDate模型——一个全局的分子钟模型。
局部的分子中模型也一样可以使用，参考examples\MouseLemurs\文件夹中的介绍。
注意我在序列名字中使用@表示序列确定的时间。
这里的文件可以被Rambaut的TipDate程序直接读取，但是如果要被baseml程序读取的话，则需要编辑一下（插入@符号）。
7）网站上下载到得文档包中还有一些其他的文件：
brown.nuc和brown.tree文件：
Brown等人在1982年报道的线粒体DNA的一个895bp的片段，这个片段在1994年由Yang等人拿来检测位点特异的突变速率模型。
mtprim9.nuc和9s.trees：
9个灵长类动物的线粒体中的一个888个碱基的DNA排列文件（Hayasakaetal.1988），这些序列还用来检测不连续的gamma模型（Yang（1994a））和自动不连续gamma模型（Yang（1995））。
abglobin.nuc和abglobin.tree：
α-和β-珠蛋白基因，这些数据用来描述密码子模型（GoldmanandYang（1994））。
abglobin.aa是序列排列后翻译成氨基酸的序列。
stewart.aa和stewart.tree：
六个哺乳动物的溶菌酶蛋白序列（Stewartetal.1987），用于预测祖先氨基酸序列（Yangetal.1995a）。
3、数据格式
3.1序列数据格式
参考下载的文档包中的一些实例数据文件（扩展名为.nuc、.aa和.nex）。
你可以把你的数据文件存储为任何一种格式，这样PAML程序们就应该可以读取这些文件。
比较合适的文件格式就是PHYLIP格式，这个格式是JoeFelsenstei开发的PHYLIP软件包的格式（Felsenstein2005）。
PAML程序们对于NEXUS文件格式的支持比较有限（PAUP程序和MacClade程序）。
对于这种格式的文件，PAML文件可以读取序列数据和树，但是命令块则忽略掉了。
PAML没有办法处理注释的部分，因此请避免使用这些。
sequential格式和interleaved格式
以下是PHYLIP程序（Felsenstein2005）的常用格式，第一行包含物种数量以及序列长度（可能随后还有一个可选性质），对于密码子序列（codeml中seqtype=1），序列的长度表示碱基的数目而不是密码子的数目。
序列中允许的选项包括I、S、P、C和G。
序列可以是interleaved格式的（选项I，例如abglobin.nuc），或者sequential格式（选项S，例如brown.nuc）。
默认的选项是S。
选项G用于分析多个基因数据，随后将会介绍。
下面就是一个sequential格式的数据实例（图1），这个实例包含4个长度为60的核苷酸（20个密码子）：
图1.sequential格式的数据
物种/序列的名字
在物种/序列的名字中不要出现一下字符：
“,:
#（）$=”，因为这会使程序感到很痛苦。
@字符可以用于序列的名字以定义序列确定的日期，例如virusl@1984，这时@会成为名字的一部分，即序列在1984年确定。
物种名字的最大字符数目在主体程序baseml.c和codeml.c的开始即定义，在PHYLIP软件中，物种的名字必须精确地为10个字符，我觉得这太严格了。
所以我用的默认的字符数目为30。
为了解决名字字符数目差异这一问题，PAML考虑在物种名字后面加上两段连续的空格符，因此物种的名字就不必是精确的30个（或者10个）。
为了符合这个规则，请不要在物种的名字中间加上两段连续的空格。
例如上述的序列数据也可以是下面的格式（图2）。
图2.图1格式的一个可变格式
如果你想让一个文件既可以被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以及它们的小写字符）都会被识别为氨基酸。
不明确的字符（未确定的核苷酸或氨基酸）也是允许的。
三个特定的字符（“.”“-”“?
”）会被认为是：
点代表当前DNA序列中的核苷酸与第一条DNA相应位置的核苷酸相同，破折号意味着序列排列时产生的gap，而问好表示不确定的核苷酸或者氨基酸。
序列中的一些非字母字符，例如<>’”$%&^{}[]0123456789，将会被程序忽略掉，因此可以用这些字符作为一些特殊的标识。
对于几个序列，每行的字符数可以不等长，你可以把所有的序列放在同一行中。
在baseml和codeml程序中，对于不明确字符和gap的处理依赖于它们的控制文件中对于变量cleandata的定义。
在最大自然法分析中，如果变量cleandata设置为1，则意味着如果一个DNA序列中包含不明确字符，那么所有序列中的相应位置的核苷酸或氨基酸都会被删除；如果cleandata=0则表示“成对删除”，即只有在序列对中含有“missingcharacter”的位点被删除。
在PAML的各个程序中，没有对于插入或者缺失突变的模型，所以在序列排列时的gap会被认为是不确定氨基酸或核苷酸（即问号）。
请注意，对于一个密码子序列，移除任何一个核苷酸意味着移除整个密码子。
注意事项可以放在序列文件的末尾，这会被程序忽略掉。
选项G：
这个选项是为了组合分析异种数据（如多个基因）而设的。
这些序列必须连接起来，选项是为了制定基因或者密码子的来源。
这个选项有三种格式：
第一中格式会用一组序列文件的一部分来说明（图3），实例数据来源于线粒体基因组的一段895bp的片段（Brownetal.（1982）），这段序列编码两个蛋白（ND4和ND5）的尾端部分，两者之间还有三个tRNA基因。
因此序列中的位点就分为四类：
三个密码子区和tRNA编码区。
文件的第一行包含选项字符“G”，第二行以字符“G”开头，随后是位点分类的数目。
Thefollowinglinescontainthesitemarks,oneforeachsiteinthesequence（oreachcodoninthecaseofcodonml）.Thesitemarkspecifieswhichclasseachsiteisfrom.Iftherearegclasses,themarksshouldbe1,2,...,g,andifg>9,themarksneedtobeseparatedbyspaces.Thetotalnumberofmarksmustbeequaltothetotalnumberofsitesineachsequence.
图3.第一种格式实例
如果数据是多个基因连接起来的序列，第二中格式就是比较有用的了，下面是这种格式的一个实例（图4）。
这个序列包括四个DNA，总长度为1000，其中第一个基因长100个核苷酸，第二个基因长200个核苷酸，第三个基因长300个核苷酸，第四个基因长400个核苷酸。
如果所有基因的长度放在一行，那么这一行必须以G开头，例如图4中第二行（这能让程序确定数据中使用的是第一种格式还是第二种格式）。
所有的基因长度的和应该等于组合序列的长度（对于baseml或basemlg程序是核苷酸，aaml程序是氨基酸，codonml程序是密码子）。
图4.第二种格式实例
第三种格式只能用于编码蛋白质的DNA序列（对于baseml程序）。
在每个数据文件第一行用字符“GC”表示（第一、第二种格式用字符“G”）。
程序在分析序列时，会把密码子的三个位置不同对待，但是需要序列的长度严格等于3的倍数。
图5.第三种格式实例
在分析序列是密码子序列时，选项G的作用（codeml程序运行时，seqtype=1）：
格式与baseml程序的同样操作相似，但是请注意序列的长度为核苷酸的数目，而基因的长度是密码子的数目。
这时常常让人出错的原因之一。
接下来是一个实例：
图6.
图6表示，这组数据包含5个序列，每个序列有300个核苷酸（100个密码子），每个序列包括2个基因，第一个基因含40个密码子，第二个基因含60个密码子。
3.2位点模式计算（Sitepatterncounts）
序列的排列可以以位点模式输入，计数时也以位点模式进行。
这种格式的特点是在输入的数据文件第一行包含一个选项字符“P”。
图7是一个实例，在这个实力中包含3个序列，8个位点模式，字符“P”表示数据是位点模式而不是位点。
对于“interleavedformat”，选项“P”和选项“I”的用法相同，而对于“sequentialforma
升级会员 