参考下载的文档包中的一些实例数据文件(扩展名为.nuc、.aa和.nex)。你可以把你的数据文件存储为任何一种格式,这样PAML程序们就应该可以读取这些文件。比较合适的文件格式就是PHYLIP格式,这个格式是Joe Felsenstei开发的PHYLIP软件包的格式(Felsenstein 2005)。PAML程序们对于NEXUS文件格式的支持比较有限(PAUP程序和MacClade程序)。对于这种格式的文件,PAML文件可以读取序列数据和树,但是命令块则忽略掉了。PAML没有办法处理注释的部分,因此请避免使用这些。
sequential格式和interleaved格式
以下是PHYLIP程序(Felsenstein 2005)的常用格式,第一行包含物种数量以及序列长度(可能随后还有一个可选性质),对于密码子序列(codeml中seqtype=1),序列的长度表示碱基的数目而不是密码子的数目。序列中允许的选项包括 I、S、P、C和G。序列可以是interleaved格式的(选项I,例如abglobin.nuc),或者sequential格式(选项S,例如brown.nuc)。默认的选项是S。选项G用于分析多个基因数据,随后将会介绍。下面就是一个sequential格式的数据实例(图1),这个实例包含4个长度为60的核苷酸(20个密码子):
物种/序列的名字
在物种/序列的名字中不要出现一下字符:“, : # ( ) $ =”,因为这会使程序感到很痛苦。@字符可以用于序列的名字以定义序列确定的日期,例如virusl@1984,这时@会成为名字的一部分,即序列在1984年确定。物种名字的最大字符数目在主体程序baseml.c和codeml.c的开始即定义,在PHYLIP软件中,物种的名字必须精确地为10个字符,我觉得这太严格了。所以我用的默认的字符数目为30。为了解决名字字符数目差异这一问题,PAML考虑在物种名字后面加上两段连续的空格符,因此物种的名字就不必是精确的30个(或者10个)。为了符合这个规则,请不要在物种的名字中间加上两段连续的空格。例如上述的序列数据也可以是下面的格式(图2)。
如果你想让一个文件既可以被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,而问好表示不确定的核苷酸或者氨基酸。序列中的一些非字母字符,例如< > ’ ” $ % & ^ { } [ ] 0 1 2 3 4 5 6 7 8 9,将会被程序忽略掉,因此可以用这些字符作为一些特殊的标识。对于几个序列,每行的字符数可以不等长,你可以把所有的序列放在同一行中。
在baseml和codeml程序中,对于不明确字符和gap的处理依赖于它们的控制文件中对于变量cleandata的定义。在最大自然法分析中,如果变量cleandata设置为1,则意味着如果一个DNA序列中包含不明确字符,那么所有序列中的相应位置的核苷酸或氨基酸都会被删除;如果cleandata=0则表示“成对删除”,即只有在序列对中含有“missing character”的位点被删除。
在PAML的各个程序中,没有对于插入或者缺失突变的模型,所以在序列排列时的gap会被认为是不确定氨基酸或核苷酸(即问号)。请注意,对于一个密码子序列,移除任何一个核苷酸意味着移除整个密码子。
注意事项可以放在序列文件的末尾,这会被程序忽略掉。
选项G:这个选项是为了组合分析异种数据(如多个基因)而设的。这些序列必须连接起来,选项是为了制定基因或者密码子的来源。
这个选项有三种格式:第一中格式会用一组序列文件的一部分来说明(图3),实例数据来源于线粒体基因组的一段895bp的片段(Brown et al. (1982)),这段序列编码两个蛋白(ND4和ND5)的尾端部分,两者之间还有三个tRNA基因。因此序列中的位点就分为四类:三个密码子区和tRNA编码区。文件的第一行包含选项字符“G”,第二行以字符“G”开头,随后是位点分类的数目。The following lines contain the site marks, one for each site in the sequence (or each codon in the case of codonml). The site mark specifies which class each site is from. If there are g classes, the marks should be 1, 2, ..., g, and if g > 9, the marks need to be separated by spaces. The total number of marks must be equal to the total number of sites in each sequence.
如果数据是多个基因连接起来的序列,第二中格式就是比较有用的了,下面是这种格式的一个实例(图4)。这个序列包括四个DNA,总长度为1000,其中第一个基因长100个核苷酸,第二个基因长200个核苷酸,第三个基因长300个核苷酸,第四个基因长400个核苷酸。如果所有基因的长度放在一行,那么这一行必须以G开头,例如图4中第二行(这能让程序确定数据中使用的是第一种格式还是第二种格式)。所有的基因长度的和应该等于组合序列的长度(对于baseml或basemlg程序是核苷酸,aaml程序是氨基酸,codonml程序是密码子)。
第三种格式只能用于编码蛋白质的DNA序列(对于baseml程序)。在每个数据文件第一行用字符“GC”表示(第一、第二种格式用字符“G”)。程序在分析序列时,会把密码子的三个位置不同对待,但是需要序列的长度严格等于3的倍数。
在分析序列是密码子序列时,选项G的作用(codeml程序运行时,seqtype=1):格式与baseml程序的同样操作相似,但是请注意序列的长度为核苷酸的数目,而基因的长度是密码子的数目。这时常常让人出错的原因之一。接下来是一个实例:
图6表示,这组数据包含5个序列,每个序列有300个核苷酸(100个密码子),每个序列包括2个基因,第一个基因含40个密码子,第二个基因含60个密码子。
序列的排列可以以位点模式输入,计数时也以位点模式进行。这种格式的特点是在输入的数据文件第一行包含一个选项字符“P”。图7是一个实例,在这个实力中包含3个序列,8个位点模式,字符“P”表示数据是位点模式而不是位点。对于“interleaved format”,选项“P”和选项“I”的用法相同,而对于“sequential format”,选项“P”与选项“S”的用法相同。图7中,数字8表示位点数量有8中模式。例如在图7中,所有的三个物种都包含G,在所有的三个物种的200位点都含有T,以此类推。总的来说,每个物种都包含100+200+40+……+14=44=440个位点。
上面的实例用于baseml和basemlg程序分析核苷酸序列。对于多个基因,需要选项P和选项G一起使用。
图7所示,为一个包含三个物种、10个位点模式、两个基因的数据。其中前四种模式是第一个基因,后面6个模式是第二个基因,一共10种模式。对于第一部分,有40个位点在三个物种之间都是碱基A,在第二部分种则有61个位点在三个物种中都是碱基A。
同样的格式也适用与蛋白序列(对于codeml程序,seqtype=2)。将图7中的碱基换成氨基酸即可。
对于密码子序列(codeml中,seqtype=1),格式如图8所示。图8中有三个物种,9种位点模式,其中6个位点是第一种位点模式(在所有三个物种中都是密码子GTG),9*3=27,即程序需要你在第一行中使用密码子位点模式数量的三倍。奇怪吧?这是为了把这种格式与“sequential format”和“interleaved format”的格式要求一致,在后面的两个程序中,序列长度定义为核苷酸的数量而不是密码子的数量。开始的时候,Yang这么做是为了同一个文件可以同时被baseml(碱基序列分析)和codonml(密码子序列分析)两个程序同时使用。
对于多基因的密码子位点模式,见图9。图9所示包含3个物种,9个核苷酸位点模式,其中前四个模式是基因1,后四个模式是基因2。
此外,选项变量P可以与选项变量I或者S一起使用。例如,PI表示位点模式的数据列出格式是“interleaved format”,而PS表示位点模式的数据列出格式是“sequential format”。而P不加S和I表示使用默认的序列格式。如果所有的序列模式都在同一行,那么这种格式既与I也与S格式相同,因此就没有比较再加上I或者S了。
如果运行baseml和codeml程序读取“sequential format”或者“interleaved format”格式的序列时,输出结果将会包含一个特定的这些特定格式的print-out。找到含有“Printing out site pattern counts”的行即可。You can move this block into a new file and later on read that file instead, if it takes a long time to pack sites into patterns. 请注意一下关于P格式的一些限制。首先,在这个选项下,有些输出是无法完成的,这些输出包括祖先序列的重建、以及当前速率下(对于序列或者序列模式)后代序列的估计,这些输出本来是可以通过设置RateAncestor=1来完成。第二,一些运算是需要序列长度的,我把序列的长度设置成了位点模式频率之和。如果位点模式的频率不是位点之和而是位点模式的概率,那么一些与序列长度相关的运算则会出现错误。这些运算包括MLEs的SEs、codonml程序中S位点和N位点数量等。
这个选项的潜在的用途:有时,我会用evolver程序模拟一段很长的序列(如>1M位点),把位点收集然后处理成位点模式需要几分钟到几个小时不等,而最大似然法分析只需要用几秒时间,而我想用不同的模型分析同样的数据,这时候就比较让人气愤了。同样,再分析一个大型基因组数据(>100M)时候也会出现这种状况。这时,你可以把模式从输出的数据文件中拷贝出来放到一个新文件中,下次运行程序时候就可以读取新文件了,这样就可以避免计算位点模式了。另外一种情况是,为了计算某个模型下位点模型的概率,随后读取概率去分析概率看是否能从一个错误的系统发育树推断出正确的树。这样,你可以检查树的重建方法是否相同。对于这种情况,你可以参考这些文献:
Debry (1992) and Yang (1994c). (I need to enable the code for printing site pattern probabilities.)