PAML软件和用法简介

Published at 2022-11-07 11:11

Author：liuyy

View：2670

PAML简介

PAML (Phylogenetic Analysis by Maximum Likelihood)，是杨子恒教授开发的一款利用DNA或Protein数据使用最大似然法进行系统发育分析的软件。

PAML手册：(http://abacus.gene.ucl.ac.uk/software/pamlDOC.pdf)

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）描述了统计和计算的详细内容。分析的例子可以从软件包中找到。

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运行。

参数设置方式与解释

baseml控制文件的关键参数设置

参数名称      参数解释
seqfile = brown.nuc 指定输入序列数据文件
treefile = brown.trees 指定输入数据的树文件
outfile = mlb 设置输出文件路径。
model = 4 指定核酸替代模型。模型0,1,…,8分别表示模型JC69， K80， F81， F84， HKY85， T92， TN93， REV (也称为 GTR)， 和 UNREST。注释参看Yang（1994 JME 39:105-111）。两个模型可以分别运行。model = 9是REV模型的特例，而model = 10是无限制模型的特例。格式如下面所示不需要过多解释。当model = 9或10时，同一行里包含的主要是指定模型的额外信息。中括号中的数字是自由参数比值的数字。例如，REV应该是5而UNREST应该是11。接在这些数字后面的是圆括号中的数字。输出文件中的比值参数将按照这里的顺序摆列。括号中没有提到的比值为1。当model = 9时，你只能指定TC或CT，但不能同时指定。当model = 10时，TC和CT是不同的。
nhomo = 0 nhomo仅针对baseml，并涉及同一替代模型中的频率参数。nhomo = 1适合一个均匀模型，但是通过最大似然迭代估算频率参数（πT，πC和πA；πG在频率总和为1中不是一个自由参数。）。这适用于F81，F84，HKY85，T92（这种情况下πGC是一个参数），TN93，或REV模型。通常（nhomo = 0）这些通过观察到的频率的平均值来估算。在两种情况下（nhomo = 0 和 1），你需要为碱基频率计算三个（或对T92计算一个）自由参数。 nhomo = 2 适用K80，F84和HKY85模型对树中的每个分枝估算一个转换/颠换比率（k）。 选项 nhomo = 3，4，和5仅用于F84，HKY85 或 T92。
fix_blength = 1 设置程序处理输入树文件中枝长数据的方法：0，忽略输入树文件中的枝长信息；-1，使用一个随机起始进行进行计算；1，以输入的枝长信息作为初始值进行ML迭代分析，此时需要注意输入的枝长信息是碱基替换率，而不是碱基替换数；2，不需要使用ML迭代计算枝长，直接使用输入的枝长信息，需要注意，若枝长信息和序列信息不吻合可能导致程序崩溃；3，让ML计算出的枝长和输入的枝长呈正比。
RateAncestor = 0 RateAncestor = 0或1。通常使用0。值1是用来使程序进行两个额外的分析。

codeml控制文件的关键参数设置

参数名称 参数解释
seqfile = input.txt 指定输入序列数据文件
treefile = input.trees 指定输入数据的树文件
outfile = mlc 设置输出文件路径。
seqtype = 1 设置输入的多序列比对数据的类型：1，密码子数据；2，氨基酸数据；3，输入数据虽然为密码子序列，但先转换为氨基酸序列后再进行分析。
model = 0 若输入数据是密码子序列，该参数用于设置branch models，即进化树各分枝的omega值的分布：0，进化树上所有分枝的omega值一致；1，对每个分枝单独进行omega计算；2，设置多类omega值，根据树文件中对分枝的编号信息来确定类别，具有相同编号的分枝具有相同的omega值，没有编号的分枝具有相同的omega值，程序分别计算各编号和没有编号的omega值。
aaRatefile = dat/jones.dat 当对蛋白数据进行分析，且model = 2时，该参数生效，用于设置氨基酸替换模型。
RateAncestor = 1 设置是否计算序列中每个位点的替换率：0，不需要；1，需要。当设置成1时，会在结果文件rates中给出各个位点的碱基替换率；同时也会进行祖先序列的重构，在结果文件rst中有体现。

baseml示例

以PAML软件自带中的数据为例，编辑baseml.ctl控制文件（设置seqfile和treefile，其他参数保持默认）。

seqfile = /path/to/paml4.9j/examples/brown.nuc
treefile = /path/to/paml4.9j/examples/brown.trees
outfile = mlb       * main result file

在将控制文件baseml.ctl作为输入，执行baseml程序。

(base) [user@server ~]# baseml baseml.ctl

mlb文件即为最后输出的结果。内容如下：

BASEML (in paml version 4.9j, February 2020)  /path/to/paml4.9j/examples/brown.nuc  HKY85 dGamma (ncatG=5)
Frequencies..
                                    T      C      A      G
Human                          0.25698 0.33073 0.30503 0.10726 GC = 0.438
Chimpanzee                     0.26480 0.32402 0.31061 0.10056 GC = 0.425
Gorilla                        0.25810 0.32514 0.31173 0.10503 GC = 0.430
Orangutan                      0.23687 0.34525 0.31508 0.10279 GC = 0.448
Gibbon                         0.24916 0.31955 0.31732 0.11397 GC = 0.434

Homogeneity statistic: X2 = 0.00429 G = 0.00428

Average                        0.25318 0.32894 0.31196 0.10592

# constant sites:    613 (68.49%)
ln Lmax (unconstrained) = -2476.965133

Distances:HKY85 (kappa)  (alpha set at 0.50)
This matrix is not used in later m.l. analysis.

Human
Chimpanzee         0.1227(42.4634)
Gorilla            0.1508(31.9110)  0.1533(28.5860)
Orangutan          0.2862(11.3980)  0.3125(13.1537)  0.3088(13.5463)
Gibbon             0.3311( 9.5815)  0.3473(10.5068)  0.3470(10.2059)  0.3428( 8.2174)

TREE #  1:  (((1, 2), 3), 4, 5);  MP score: 357.00
lnL(ntime:  7  np:  9):  -2621.554340    0.000000
   6..7     7..8     8..1     8..2     7..3     6..4     6..5
 0.116030 0.033500 0.059487 0.074846 0.078037 0.296862 0.449334 22.975306 0.229124

tree length =   1.10810

(((Human, Chimpanzee), Gorilla), Orangutan, Gibbon);

(((Human: 0.059487, Chimpanzee: 0.074846): 0.033500, Gorilla: 0.078037): 0.116030, Orangutan: 0.296862, Gibbon: 0.449334);

Detailed output identifying parameters

Parameters (kappa) in the rate matrix (HKY85) (Yang 1994 J Mol Evol 39:105-111):
 22.97531

alpha (gamma, K=5) =  0.22912
rate:   0.00048  0.01948  0.16074  0.75101  4.06828
freq:   0.20000  0.20000  0.20000  0.20000  0.20000

codeml示例

计算选择压力(dN/dS)

以PAML软件自带中的数据为例，设置控制文件lysozymeSmall.ctl中的seqfile和treefile，设置参数model = 0，将结果文件命名为mlc_small_model0，其他参数保持默认。

seqfile = /path/to/paml4.9j/examples/lysozyme/lysozymeSmall.txt
treefile = /path/to/paml4.9j/examples/lysozyme/lysozymeSmall.trees
outfile = mlc_small_model0

在将控制文件lysozymeSmall.ctl作为输入，执行condeml程序。

(base) [user@server ~]# codeml lysozymeSmall.ctl

再次修改控制文件，将输出文件更改为mlc_small_model2，并将参数model设置为model = 2，其他参数保持默认。

再次执行condeml程序。

(base) [user@server ~]# codeml lysozymeSmall.ctl

两次运行codeml程序的结果分别记录在mlc_small_model0和mlc_small_model2两个文件中。

打开mlc_small_model0（零假设）文件，找到以下内容。

lnL(ntime: 11  np: 13):   -906.017440      0.000000

打开 mlc_small_model2（备择假设） 文件，找到以下内容。

lnL(ntime: 11  np: 14):   -903.076551      0.000000

用R计算likelihood ratio test，将lnL二者差值的二倍 ((906.017440-903.076551)x2) 与自由度为(14-13)(df=1)的卡方分布进行比较。

> pchisq((906.017440-903.076551)*2,1,lower.tail=FALSE)
[1] 0.01529837

P = 0.01529837 < 0.05，说明相比于 mlc_small_model0，mlc_small_model2更符合数据，选择备择假设。

PAMLX是PAML的图形界面版，利用PAMLX计算选择压力，可参考**codeml的理论和实践

视频教程可参考*利用codeml计算dN/dS*

重建祖先序列

以PAML软件自带中的数据为例，设置codeml.ctl控制文件中的seqfile和treefile，其他参数保持默认。

seqfile = /path/to/paml4.9j/examples/stewart.aa    * sequence data filename
treefile = /path/to/paml4.9j/examples/stewart.trees    * tree structure file name
outfile = mlc           * main result file name

在将控制文件codeml.ctl作为输入，执行condeml程序。

(base) [user@server ~]# codeml condeml.ctl

打开rst文件找到以下内容

tree with node labels for Rod Page's TreeView
(((1_Langur, 2_Baboon) 9 , 3_Human) 8 , 4_Rat, (5_Cow, 6_Horse) 10 ) 7 ;

Nodes 7 to 10 are ancestral

List of extant and reconstructed sequences

    10    128

Langur            KIFERCELAR TLKKLGLDGY KGVSLANWVC LAKWESGYNT EATNYNPGDE STDYGIFQIN SRYWCNNGKP GAVDACHISC SALLQNNIAD AVACAKRVVS DQGIRAWVAW RNHCQNKDVS QYVKGCGV
Baboon            KIFERCELAR TLKRLGLDGY RGISLANWVC LAKWESDYNT QATNYNPGDQ STDYGIFQIN SHYWCNDGKP GAVNACHISC NALLQDNITD AVACAKRVVS DQGIRAWVAW RNHCQNRDVS QYVQGCGV
Human             KVFERCELAR TLKRLGMDGY RGISLANWMC LAKWESGYNT RATNYNAGDR STDYGIFQIN SRYWCNDGKP GAVNACHLSC SALLQDNIAD AVACAKRVVR DQGIRAWVAW RNRCQNRDVR QYVQGCGV
Rat               KTYERCEFAR TLKRNGMSGY YGVSLADWVC LAQHESNYNT QARNYDPGDQ STDYGIFQIN SRYWCNDGKP RAKNACGIPC SALLQDDITQ AIQCAKRVVR DQGIRAWVAW QRHCKNRDLS GYIRNCGV
Cow               KVFERCELAR TLKKLGLDGY KGVSLANWLC LTKWESSYNT KATNYNPSSE STDYGIFQIN SKWWCNDGKP NAVDGCHVSC SELMENDIAK AVACAKKIVS EQGITAWVAW KSHCRDHDVS SYVEGCTL
Horse             KVFSKCELAH KLKAQEMDGF GGYSLANWVC MAEYESNFNT RAFNGKNANG SSDYGLFQLN NKWWCKDNKR SSSNACNIMC SKLLDENIDD DISCAKRVVR DKGMSAWKAW VKHCKDKDLS EYLASCNL
node #7           KVFERCELAR TLKRLGMDGY RGISLANWVC LAKWESNYNT QATNYNPGDQ STDYGIFQIN SRYWCNDGKP GAVNACHISC SALLQDNIAD AVACAKRVVR DQGIRAWVAW RNHCQNRDVS QYVQGCGV
node #8           KVFERCELAR TLKRLGMDGY RGISLANWVC LAKWESGYNT QATNYNPGDQ STDYGIFQIN SRYWCNDGKP GAVNACHISC SALLQDNIAD AVACAKRVVR DQGIRAWVAW RNHCQNRDVS QYVQGCGV
node #9           KIFERCELAR TLKRLGLDGY RGISLANWVC LAKWESGYNT QATNYNPGDQ STDYGIFQIN SRYWCNDGKP GAVNACHISC SALLQDNIAD AVACAKRVVS DQGIRAWVAW RNHCQNRDVS QYVQGCGV
node #10          KVFERCELAR TLKRLGMDGY RGISLANWVC LAKWESNYNT QATNYNPGDE STDYGIFQIN SKWWCNDGKP GAVNACHISC SELLEDNIAD AVACAKRVVR DQGITAWVAW RNHCQDRDVS QYVQGCGL

node #7、node #8 、node #9、node #10 即为预测的祖先蛋白序列。

PAMLX是PAML的图形界面版，可参考codeml重建祖先序列