经验大分享：PAML简介

PAML：

用最大似然法分析系统发育

目录

1.简介

　　PAML能做什么

　　PAML不能做什么

　　如何运行PAML中的程序

2.数据文件格式

　　序列文件格式

　　树文件格式

3.Baseml

　　核酸替代模型

　　控制文件

4.Basemlg

5.Codeml

　　密码子替代模型

　　氨基酸替代模型

　　控制文件

　　密码子序列（seqtype =1）

　　氨基酸序列（seqtype = 2）

6.Evolver

7.Yn00

8.Mcmctree

　　简　　介

　　化石校准

　　控制文件

MCMCtree 和 Multidivtime 的差异

　　近似似然计算

9.其他注释

　　分析大数据集和迭代算法

　　树形检索算法

　　生成引导数据集

　　迭代过程中的磨损记录

　　指定初始值

　　微调迭代算法

　　源程序中可调节参数

PAML与其他系统发育软件的共用

PHYLIP

PAUP、MacClade 和 MrBayes

Clustal

MEGA

MOLPHY

TreeView

1.简介

PAML能做什么

当前的PAML软件包有以下程序：baseml、basemlg、codeml、evolver、pamp、yn00、mcmctree和chi2. 可以用这个软件包进行的分析包括：

系统进化树的比较和检验（baseml 和 codeml）；

复杂替代模型的参量估计，包括sites之间的变量比率的模型和多重基因或者site划分的组合分析模型（baseml 和 codeml）；

应用模型比较假设的似然比率检验（baseml、codeml、chi2）；

基于全局和局部的分子钟的分歧时间的估计（baseml 和 codeml）；

用核酸、氨基酸和密码子模型进行的祖先序列似然重建（先验贝叶斯）（baseml 和 codeml）；

由蒙特卡洛仿真产生核酸、密码子、氨基酸序列数据集（evolver）；

蛋白质编码序列中正向选择的同义与非同义替代的估计和检测（yn00 和 codeml）；

在化石校准中的分歧时间合并的不确定性的贝叶斯估计（mcmctree）.

PAML的长处在于它收集了复杂的替代模型.树形检索算法在baseml和codeml中的应用时相当粗糙的，因此除了非常小的数据集（<10 species），最好用别的软件包来推断拓扑树.比如：phylip，paup，mrBayes.

Baseml 和 codeml

Baseml主要用于核酸序列的最大似然分析. Codeml 是由codonml 和 aaml合并而来，其中codonml 是编码蛋白序列的密码子替代模型的应用，而aaml是氨基酸序列模型的应用.现在这两个软件通过控制文件codeml.ctl中的序列类型来区分，1为密码子序列，2为氨基酸序列.baseml、codonml和aaml使用相似的算法来适用最大似然模型，主要区别在于马尔科夫模型中的进化单元，在序列中被称为“site”，对这三个软件来说分别是是一个核苷酸，一个密码子或者一个氨基酸.

Evolver

这个软件可以在核酸、密码子和氨基酸替代模型下用来仿真序列.他还有一些其他选项，比如得到一个随机树，和计算树之间的分歧距离.

Basemlg

这个软件应用（连续）gamma模型.它非常的慢并且对于大于6或7个物种的数据是不能使用的.在baseml中的离散-gamma模型是可以使用的.

Macmctree

这个软件应用贝叶斯MCMC算法来估算物种分歧时间.

Pamp

基于简约法的分析.

Yn00

通过编码蛋白DNA序列之间的比较来评估同义和非同义替代的比率（ds 和dn）.

Chi2

用来进行似然比率检验.它计算chi平方的临界值，可以用来和检验统计量（由真实数据计算得来）进行比较来确定这个实验实在5%或者1%水平.通过软件的名字“chi2”来适用该软件，当你输入检验统计量和d.f时它同时开还可计算P值，命令为“chi2 p”.

PAML不能做什么

序列比对.

基因预测.

在大量数据集里进行树形检索.

2.数据文件格式

PAML的原始数据格式是PHYLIP软件包里用到的PHYLIP格式.PAML对PAUP和MacClade使用的NEXUS格式是有限制的.只有序列数据和树是可识别的，其他的会被忽略掉.PAML不能识别学列数据文件里的注释语句块，所以请避免使用.

序列和交叉格式

PHYLIP格式：

4 60

sequence 1

AAGCTTCACCGGCGCAGTCATTCTCATAAT

CGCCCACGGACTTACATCCTCATTACTATT

sequence 2

AAGCTTCACCGGCGCAATTATCCTCATAAT

CGCCCACGGACTTACATCCTCATTATTATT

sequence 3

AAGCTTCACCGGCGCAGTTGTTCTTATAAT

TGCCCACGGACTTACATCATCATTATTATT

sequence 4

AAGCTTCACCGGCGCAACCACCCTCATGAT

TGCCCATGGACTCACATCCTCCCTACTGTT

第一行包含了物种的个数以及序列的长度.对编码序列（codeml seqtype = 1），序列的长度是指核酸的个数而不是密码子的个数.这个序列文件中只允许以下选项：I，S，P，C，G.序列可以是层叠格式（选项I），或者连续的格式（选项S）.默认的选项是S.选项G用于多重基因数据的联合分析（如上示例）.物种/序列名称.不能使用下面的特殊符号：“，：#（）￥=”，因为他们有特殊的用途，以免混淆，符号@可以为序列名字的一部分或者结尾用来区分数据说明，例如：virusl@1984.名字的字符个数在PHYLIP中的baseml.c和codeml.c的开头进行设置，为10个字符，而PAML默认为30.物种名字的结尾应该有两个连续的空格.//代码效果参考：http://www.zidongmutanji.com/zsjx/19471.html

4 60

sequence 1 AAGCTTCACCGGCGCAGTCATTCTCATAAT

CGCCCACGGACTTACATCCTCATTACTATT

sequence 2 AAGCTTCACCGGCGCAATTATCCTCATAAT

CGCCCACGGACTTACATCCTCATTATTATT

sequence 3 AAGCTTCACC GGCGCAGTTG TTCTTATAAT

TGCCCACGGACTTACATCATCATTATTATT

sequence 4 AAGCTTCACCGGCGCAACCACCCTCATGAT

TGCCCATGGACTCACATCCTCCCTACTGTT

在一个序列中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 或者他们的小写形式）被识别为氨基酸.没有定义为核酸或者氨基酸符号的字符也是允许出现在序列里的.特殊的符号“.”“，”“-”“？”被识别为：“.”与第一条序列相同的字符，“，”为比对gap，“？”为没有定义的字符.非字母符号“> < ! ' " $ % ^ 【 】 ( ) { } 0 1 2 3 4 5 6 7 8 9在序列中直接被忽略或者作为一个指示符号. 模糊字符在控制文件里根据变量cleandata决定怎么处理. 在最大似然分析之前，如果cleandata=1，那么含有模糊字符的所有序列的该位点都被移除.如果cleandata=0，那么模糊字符和比对gaps都被视为模糊字符.在计算序列对之间的距离时（输出中的下三角距离矩阵），cleandata=1意味着“完全删除”，cleandata=0意味着“配对序列删除”.

在PAML软件中没有插入和删除模型.因此一个比对gap被视为一个模糊字符“？”.对密码子序列删除某个核酸以为这删除整个密码子.

说明部分可以放在序列文件的后面，而且会被软件忽略掉.

选项G：这个选项是为了不均匀数据集合的联合分析，比如多基因数据或者三个密码子位置的数据. 序列必须连接在一起并且这个选项用来识别基因或者密码子位置的来源.这个选项有三种格式.第一种格式如下：

5 895 G

G 4

3

123123123123123123123123123123123123123123123123123123123123

123123123123123123123123123123123123123123123123123123123123

123123123123123123123123123123123123123123123123123123123123

123123123123123123123123123123123123123123123123123123123123

123123123123123123123123123123123123123123123123123123123123

123123123123123123123123123123123123123123123123123123123123

123123123123123123123123123123123123123123123123123123123123

1231231231231231231231231231231231231

444444444444444444444444444444444444444444444444444444444444

444444444444444444444444444444444444444444444444444444444444

444444444444444444444444444444444444444444444444444444444444

444444444444444444

123123123123123123123123123123123123123123123123123123123123

123123123123123123123123123123123123123123123123123123123123

123123123123123123123123123123123123123123123123123123123123

12312312312312312312312312312312312312312312312312312312312

Human

AAGCTTCACCGGCGCAGTCATTCTCATAATCGCCCACGGACTTACATCCTCATTACTATT

CTGCCTAGCAAACTCAAACTACGAACGCACTCACAGTCGCATCATAATC........

Chimpanzee

.........

这组数据是线粒体上的一个895bp的片段，它的两端分别编码两个蛋白的部分（ND4和ND5），中间编码三个tRNAs，序列的sites自然分成了4个class：三个编码位点，和tRNA编码区域.第一行包含符号G.第二行以G开头，后跟着一个class数.下一行包含了site标记，用来区分site来源.如果有g个classes，标记为1,2,……g，如果g>9，标记将以空格隔开.

第二种可用格式是多重基因的连续序列，如：

5 1000 G

G 4 100 200 300 400

Sequence 1

TCGATAGATAGGTTTTAGGGGGGGGGGTAAAAAAAAA.......

这个序列来自4个基因的1000个核苷酸，有连续的4个基因片段100,200,300,400链接起来.基因的长度必须跟在G后.

第三种格式适合蛋白编码DNA序列（baseml）.你可以在第一行使用GC来代替G.在核酸分析上将会完全不同，它假设序列的长度刚好是3的倍数.

5 855 GC

human GTG CTG TCT CCT ...

选项G对编码序列（codeml seqtype=1）.这个格式和baseml的格式相似，但注意序列的长度是核酸的数目，而基因的长度是密码子的数目：

5 300 G

G2 40 60

这个集合有5条序列，每条有300个核苷酸（100 codons），分成两个基因，第一个含有40个密码子，第二个有60个密码子.

Site patterns计数

序列的比对同样可以是site patterns的输入.这个格式用数据文件第一行的选项“P”区分如：

3 8 P

human GTACTGCC

rabbit GTACTACT

rat GTACAGAC

100 200 40 50 11 12 13 14

这里有三条序列，八个site，“P”说明这是一个site patterns 而不是sites，“P”的用法和层叠格式的“T”以及连续序列的“S”（默认）用法一样.比对结果下面的8个数字是上面site的8个模式.例如，在100 sites，三个位点都含有G，而在200 sites所有位点都有T，等等.总共有100+200+40+50+11+12+13+14=440sites.

这个例子适用于baseml和basemlg，基于核酸分析的软件，为了区分多基因（site 划分），可以G和P一起使用：

3 10 PG

G 2 4 6

human GTTA CATGTC

rabbit GTCA CATATT

rat GTTA CAAGTC

100 200 40 50 120 61 12 13 54 12

这里有10个site模式和2个基因，前四个模式是第一个基因的后面6个是第二个基因的.

同样的格式适用于蛋白序列（codeml seqtype=2），用氨基酸代替上述的核酸.

对编码序列（codeml seqtype=1），格式如下：

3 27 P G

human GTG CTG TCT CCT GCC GAC AAG ACC

rabbit ... ... ... G.C ... ... ... T..

rat ... ... ... ..C ..T ... ... ...

6 1 1 1 1 4 3 1 1

有三个物种，九个模式，六个sites为第一个模式（三个物种中都有GTG密码子）.注意到27 = 93. 这个软件需要你用密码子site patterns的三倍.这与连续序列或者层叠序列格式一致，序列的长度表示核苷酸的个数而不是密码子的个数.

区分多基因密码子site patterns，如下：

3 27 P G

G 2 4 5

human GTG CTG TCT CCT GCC GAC AAG ACC

rabbit ... ... ... G.C ... ... ... T..

rat ... ... ... ..C ..T ... ... ...

6 1 1 1 1 4 3 1 1

此外，选项P可以和I或S一起使用. PI意味着site patterns用层叠的格式列出而PS 以为以连续序列的格式列出（默认）.

这些选项的限制：

有一些内容将会没有输出，包括祖先序列重建和sites 比率的估计，你可以通过使用RateAncestor=1来得到序列. 其次，一些计算需要序列的长度，设置为site pattern的频率.如果site pattern频率不是sites 的个数而是site pattern 概率，计算序列的长度将不准确.这些计算包括：MLEs 的 SEs，codonml 中sites S和Ｎ数量的计算.

这些选项可能的使用. 有时我用 evolver来仿真非常长的序列（>1M sites)并且它可以花几分钟或几个小时把setes 转换成 patterns，那是非常气愤的，当我想用一样的数据用最大似然迭代花几秒来使用多重模型时。一个相似的情况是长序列（>100Mb sites)大基因数据的分析.这种情况，你可以只运行baseml 或者 codeml一次，然后从输出文件中把pattern 数量复制出来. 下一次，你运行这个软件你可以读这个新文件. 这样，软件就会跳过site patterns的计算. 另外一种情形是在模型下计算site pattern 概率，然后读概率来分析，用一个错误来看这个正确的树是否依然被恢复，这样你可以检查这个重构树的方法是否一致.

树文件格式 和 拓扑树的表示法

当runmode=0 或 1时，将会用到树结构的文件.这个文件的名字在适当的控制文件里指出.拓扑树一般用圆括号记法.

4 5 // 4 species, 5 trees

(1,2,3,4); // the star tree

((1,2),3,4); // species 1 and 2 are clustered together

((1,2),3,4); // Commas are needed with more than 9 species ((human,chimpanzee),gorilla,orangutan); ((human:.1,chimpanzee:.2):.05,gorilla:.3,orangutan:.5);

如果树有分支的长度，baseml和codeml在最大似然迭代中允许在树中使用分支长度作为起点值.

你使用的是有根树还是无根树依赖于模型，例如，是否已经假设了分子钟.没有分子钟（clock=0），将会用无根树，如（（1,2）3,4）或（1,2，（3,4））.对分子钟或局部分子钟模型，将会用有根树，这两个树是不同的而且都不同于（（（1,2），3），4）.

在PAML中，一个有根数在根部有一个分歧点，而无根树则在根部有分叉点.

由PAUP和MacClade产生的树文件. PAML软件对这两个软件生成的树文件只有有限的兼容性. 首先“【 U】”符号记法用来识别一个无根树将被忽略. 对于有PAML得到可接受的无根树，你需要手动修改树文件使得在根部有一个分歧点，例如，通过修改“（（（1,2），3,），4）” 为 “（（1,2），3,4）”.其次，“转化”关键字同样被PAML忽略，并且假设树中序列的顺序和序列文件中的顺序一样.

分支或节点标签. baseml和codeml 中应用的一些模型允许一些书的分支组，设为不同的参数. 例如，在baseml或codeml中的局部分子钟模型（clock=2 或3）中，你可以有3个分支组的比率，慢，中等，和高比率. 同样特有的编码模型（codonml 中model =2 或3）允许不同的分支组有不同的ws，导致所谓的“双比率”和“三比率”. 所有的这些模型需要在树分支或节点处加上标签. 分支标签和分支长度一样用符号“#”而不是“：”.分支标签是由0开始的连续的整数（默认）.例如，下面的树：

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

分支标签适用于不同分支适用不同的w. 人和长臂猿的内部的祖先分支比率为w1，其他的分支（默认标签 #0）为w0. 这适用于小数据集合.

对一个非常大的树，你可能想在一个进化支内标记所有的分支. 为此，你可以使用进化支标签$. $是对于Δ ，看起来想一个很好的进化支标记但是在大多的关键字上会丢失.因此$2

等同于在进化支#2内标记所有的分支/节点. 下面是两个等价的树：

(((rabbit, rat) $1, human), goat_cow, marsupial);

(((rabbit #1, rat #1) #1, human), goat_cow, marsupial);

这里是考虑到网状进化支标签的规则. #号优先于$，并且靠近tips的进化标签优先于进化标签对靠近与根部的祖先节点.因此下面的两个数是等价的.在第一个树中$1首先应用到胎盘哺乳动物的所有的进化支，然后$2应用到兔子的所有进化支.

((((rabbit, rat) $2, human #3), goat_cow) $1, marsupial);

((((rabbit #2, rat #2) #2, human #3) #1, goat_cow #1) #1, marsupial);

注意到在这个规则下，无论你是否包含标签$0，它会不同. 例如：

((a, b) $0, (c, d)) $1;

左边标记三个分支（a，b，和 祖先分支a和b）就像#0，而别的分支为#1.然而下面将会把所有分支标记为#1.

((a, b), (c, d)) $1;

用TreeView方法创建树文件和读取树是非常方便的，然而，你需要把标签放在一个引号里面，如：

((((rabbit '#2', rat '#2') '#2', human '#3') '#1', goat_cow '#1') '#1', marsupial);

这样，就可以被TreeView和TreeView X（baseml/codeml)读取. 注意到TreeView和TreeView X不能接受tips或者tip 分支的标签，并且可能打乱序列名字的标签.

分歧时间符号@.化石检验信息用符号@指出. 这个在baseml和codeml中的分子钟和局部分子中模型中用到. 如：人-黑猩猩分歧时间约在5MY.

((gorilla, (human, chimpanzee) '@0.05'), orangutan);

拓扑树的分支表现：PAML中用到的第二种表示拓扑树的方法是计算它的分支，每一个用一个起点和终点节点表示. 这种表示同样用在估计分支长度的输出文件中，但你可以在树文件中使用，如树（（1,2），3,4）可以通过它的5个分支表示：

5

5 6 6 1 6 2 5 3 5 4

在这个树中，节点用连续的自然数做索引，1,2，……，s 表示已知数据中第s条已知的序列，和数据中的顺序一致．大于s 标签的是内部节点，这个序列是未知的.因此才上面的树中，节点5是祖先到6,3,4节点而6是祖先到1,2节点.

这个种符号记法能很方便的区分数据剧烈中的祖先到其他序列. 例如，下面的树有5条序列4个分支，5为1,2,3,4共同的祖先：

4

5 1 5 2 5 3 5 4

这种记法并没有在baseml和codeml中的所有模型中应用，所以使用时要注意.

3.Baseml

核苷酸替代模型

PAML中用到的模型包括JC69（Jukes and Cantor 1969），K80（Kimura 1980），F81（Felsenstein 1981），F84（Felsenstein，DNAML 软件1984，PHYLIP 版本2.6），HKY85（Hasegawa 等1984；Hasegawa等1985），Tamur（1992），Tamura 和Nei（1993），和REV，又名GTR（general-time-reversible）.比率矩阵如下：

控制文件

Baseml默认的控制文件是baseml.ctl，（下图）.

注意到等号的两边都需要空格，并且注释行以开头，没有用到的选项可以删除，变量的顺序是没有影响的.

seqfile = brown.nuc sequence data file name

outfile = mlb main result file

treefile = brown.trees tree structure file name

noisy = 3 0,1,2,3: how much rubbish on the screen

verbose = 0 1: detailed output, 0: concise output

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

3: StepwiseAddition; (4,5):PerturbationNNI

model = 5 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85

5:T92, 6:TN93, 7:REV, 8:UNREST, 9:REVu; 10:UNRESTu

Mgene = 0 0:rates, 1:separate; 2:diff pi, 3:diff kapa, 4:all diff

ndata = 1 number of data sets

clock = 0 0:no clock, 1:clock; 2:local clock; 3:CombinedAnalysis

fix_kappa = 0 0: estimate kappa; 1: fix kappa at value below

kappa = 2.5 initial or fixed kappa

fix_alpha = 1 0: estimate alpha; 1: fix alpha at value below

alpha = 0. initial or fixed alpha, 0:infinity (constant rate)

Malpha = 0 1: different alpha's for genes, 0: one alpha

ncatG = 5 # of categories in the dG, AdG, or nparK models of rates

fix_rho = 1 0: estimate rho; 1: fix rho at value below

rho = 0. initial or fixed rho, 0:no correlation

nparK = 0 rate-class models. 1:rK, 2:rKfK, 3:rKMK(1/K), 4:rKMK

nhomo = 0 0 1: homogeneous, 2: kappa for branches, 3: N1, 4: N2

getSE = 0 0: don't want them, 1: want S.E.s of estimates

RateAncestor = 0 (0,1,2): rates (alpha>0) or ancestral states

Small_Diff = 1e-6

cleandata = 1 remove sites with ambiguity data (1:yes, 0:no)?

icode = 0 (RateAncestor=1 for coding genes, "GC" in data)

fix_blength = 0 0: ignore, -1: random, 1: initial, 2: fixed

method = 0 0: simultaneous; 1: one branch at a time

seqfile，outfile，和treefile 指定序列数据文件，主要输出文件和树结构文件.在文件名里不能有空格，一般避免使用特殊字符.

Noisy控制在屏幕输出的多少. 如果选择的模型包含大量的计算，你可以选择一个大的值给noisy来避免孤独.verbose控制输出到结果文件的多少.

runmode=0意味这个在拓扑结构树文件里详细说明估计值，而runmode=1 或 2意味着用星-分解算法进行启发式树形检索，而runmode=2算法从星状树开始，当用runmode=1时，软件将会从树结构文件里读一个多叉树并估计与之一致的最好的树的分支.runmode=3意味着阶梯式的增加.runmode=4意味着对由最大简约法得到的树NNI摄动，runmode=5意味着对树结构文件NNI摄动.树形检索工作的不理想，所以尽量使用runmode=0.对于小的数据集合可用阶梯增加算法.

Model指定核苷酸替代模型.模型0,1，...，8代表模型JC69，K80，F81，F84，HKY85，T92，TN93，REV（GTR），和UMREST.最近补充了两个模型model=9指定REV模型的特殊情况，而model=10指定无约束模型的情形. 格式如下：

model = 10 【0】 / JC69 /

model = 10 【1 (TC CT AG GA)】 / K80 /

model = 10 【11 (TA) (TG) (CT) (CA) (CG) (AT) (AC) (AG) (GT) (GC) (GA) 】 / unrest /

model = 10 【5 (AC CA) (AG GA) (AT TA) (CG GC) (CT TC)】 / SYM /

model = 9 【2 (TA TG CA CG) (AG)】 / TN93 */

当指定model=9或10的时候你最好注释额外的信息.中括号里面的数字是自由比率参数的个数.例如，对REV为5或UNREST为11.数字等于圆括号对的个数. 输出文件的比率参数也为这个数.没有提到的为比率1.当model=9，你指定TC或CT，但不要两个都指定. 当model=10，TC和CT是不同的.

Mgene在序列数据文件中和选项G结合使用，用于联合分析.这样的模型也叫分隔模型.如果在数据文件中没有用到G则值为0.

这些模型总结在表1. 最简单的模型假设在基因中完全同源可以用来把不同的基因连接程一条序列，不适用G选项（并指定Mgene=0）.最常用的一个模型等价于分开分析，并且可以把这个模型用到每个基因，同样也可以在合并的数据文件中指定Mgene=1和G.不同基因的最大似然的和最常用模型的似然值都做考虑. 考虑到不同源的多基因模型，通过指定Mgene和G.Mgene=0意味着假设不同基因间相同模式替代的不同的替代比率.Mgene=2意味着对不同基因但是相同比率参数不同频率的参数（k in the K80，F84，F80，HKY85模型 或者TN93和REV模型中的比率参数）.……