Published at 2022-11-07 11:11
Author:liuyy
View:2827
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来评估它们。
参数名称 | 参数解释 |
---|---|
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是用来使程序进行两个额外的分析。 |
参数名称 | 参数解释 |
---|---|
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中有体现。 |
以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
以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重建祖先序列