PAML操作.docx

PAML操作.docx

《PAML操作.docx》由会员分享，可在线阅读，更多相关《PAML操作.docx（7页珍藏版）》请在冰豆网上搜索。

PAML操作

PAML软件的一些简单的具体的使用操作

1.首先用ClustalX进行序列比对：

要保证：

保证核苷酸序列是三的倍数，没有终止密码子，核苷酸序列的第一位是密码子的第一位。

假设序列名为cox1.fas

2.使用DAMBE软件进行转换成PML格式。

打开要换换的文件，然后“file”“saveandconvertsequenceformat”，在保存类型中选择“Yang’sPAML”。

那么此时的序列名为“cox1.PML”，这样就可以得到文件“*.PML”，然后就直接把后缀改成“*.nuc”。

那么此时的序列名为“cox1.nuc”这样就完成了文件格式的转换。

3.打开PAML软件的文件夹，找到文件名是“bin”的文件夹，打开之后，找到程序“codeml.exe”,把该程序复制到D盘的根目录下。

（这一步并不是必要的，只是要把用到的几个程序放在同一个目录下）

4.在你使用ClustalX进行序列比对的时候，会生成一棵进化树，适用treeview软件可以打开，你需要的是把文件的后缀名改称“*.trees”。

即树的文件名是“cox1.trees”，这就完成了树的格式的转换。

5.然后再PAML4的文件夹中找到一个后缀是“*.ctl”的文件，把文件名改成“cox1.ctl”，复制到和“codeml.exe”相同的地方。

6.要对codeml.ctl文件中的各个选项的值进行修改，具体内容如下：

  seqfile=cox1.nuc按你自己的文件名进行修改，就可以了，

  treefile=cox1.trees

  outfile=mlc*mainresultfilename，

  noisy=9*0,1,2,3,9:

howmuchrubbishonthescreen，

  verbose=0*0:

concise;1:

detailed,2:

toomuch

  runmode=0

  seqtype=1*1:

codons;2:

AAs;3:

codons-->AAs

  CodonFreq=2*0:

1/61each,1:

F1X4,2:

F3X4,3:

codontable

  clock=0

  aaDist=0*0:

equal,+:

geometric;-:

linear,1-6:

G1974,Miyata,c,p,v,a

  aaRatefile=wag.dat*onlyusedforaaseqswithmodel=empirical（_F）*

               dayhoff.dat,jones.dat,wag.dat,am.dat,oryourown

  model=0，这是使用的最简单的模型，*modelsforcodons:

*0:

one,1:

b,2:

2ormoredN/dSratiosforbranches

  *modelsforAAsorcodon-translatedAAs:

  *0:

poisson,1:

proportional,2:

Empirical,3:

Empirical+F28*6:

FromCodon,7:

AAClasses,

    8:

REVaa_0,9:

REVaa（nr=189）

  NSsites=031278，依次选取了6个模型。

也可以选其中的两个，但必须是0和3，1和2，7和8。

相互配对

       icode=4*0:

universalcode;1:

mammalianmt;2-10:

seebelow如果是核基因的话就选0。

       fix_kappa=0

       kappa=5

       fix_omega=0

       omega=0.2

       getSE=0

       RateAncestor=0

       Small_Diff=.5e-6

       cleandata=1*removesiteswithambiguitydata（1:

yes,0:

no）?

       method=0*0:

simultaneous;1:

onebranchatatime

   最后保证在同一个文件夹内同时具有：

三个文件“codeml.exe”,“cox1.nuc”“cox1.trees”，这时候你双击codeml.exe，就可以运行程序。

如果不能正确运行的话，你可以通过运行cmd，在dos情况下，运行codeml.exe，这样会有错误提示，知道你错在哪里了。

常见命令解释：

1.Baseml.ctl的命令说明：

 noisy用来控制输出结果的多少，如果模型适用的运算比较多的话，noisy的值可以选择的比较大，verbose可以控制结果文件中结果的多少。

 runmode=0表明在树的结构文件中估算树的拓扑结构。

这个选项是我们通常情况下选择的，基本上可以满足我们的需要。

 Runmode=1or2表明通过星状-分解算法来进行启发式搜索树。

Runmode=2这种算法是从星状树开始搜索，而runmode=1则表明软件读取多歧树是从树的结构文件中，并且同过比较去估计最佳二歧树。

runmode=3表明是逐步增加的。

runmode=4通过简约法来搜索具有NNIperturbation的起始树。

runmode=5表明从树的结构文件中来读取NNIperturbationwith起始树。

 Model0，1，…，8分别代表以下模型：

JC69,K80,F81,F84,HKY85,T92,TN93,REV（alsoknowasGTR）,andUNREST。

 Mgene用于和序列数据文件中的optionG进行联合，用于多个基因和多个位点的联合分析。

如果不使用optionG的话，则选择0。

 ndata用于指定文件中的分隔的数据集的数目。

它的变化被用于模拟，你可以使用evolver来产生200个复制数据集，这是设置ndata=200，然后用baseml进行分析。

 clock用于指定谱系之间速率恒定或变化的模型。

Clock=0，意味着整棵树中，不同分支之间不存在clock现象；Clock=1，意味着globalclock，所有的分支具有相同的进化速率；clock=2意味着localclock，所有分支之间的进化速率被分成几个部分；clock=3意味着多个基因或多重分隔数据，允许分支的进化速率以不同的方式变化。

；；

2.Codeml.ctl的使用说明：

 CodonFreq用于平衡密码子替换模型中的密码子使用频率。

Codonfreq=0说明每种密码子的使用频率是相同的；codonfreq=1说明是从平均核苷酸频率中计算出来的；codonfreq=2说明是从三个密码子位置的平均核苷酸频率得来的；codonfreq=3则使用了三个参数。

Codonfreq=0，1，2和3所代表的模型中使用的参数的数目分别为：

0，3，9，60。

 aadist用于指定氨基酸距离是否是相同的（=0），还是属于Grantham’smatrix（=1）。

 runmode=-2 如果是核酸序列执行ML方法来推测蛋白编码序列两两之间的dn和ds，如果是氨基酸序列，则是计算ML距离的。

 model用于估计各个分支之间的w值。

Model=0，表明所有的谱系具有一个w比率（onewratio）；model=1，每一支具有一个速率（free-ratio）；model=2表明速率的任意数字。

 NSsites主要是用于指定模型允许dn/ds（w）在不同的位点之间变化。

NSsites=m表明对应于model=m。

变化的ncatG被用来指定在一些特定的模型下的w分布的类型的数目。

NcatG的值被用于执行一下分析：

paperare3forM3（discrete）,5forM4（freq）,10forthecontinuousdistributions（M5gamma,M6:

2gamma,M7:

beta,M8:

beta&w,M9:

beta&gamma,M10:

beta&gamma+1,M11:

beta&normal>1,andM12:

0&2normal>1,M13:

3normal>0）.ThismeansM8willhave11siteclasses（10fromthebetadistributionplus1additionalclass）。

通过NSsites可以同时执行多个模型，例如：

NSsites=012378，的意思就是同时执行M0，M1，M2a，M3，M7，和M8。

整理建议：

使用M1a和M2a来重建LRT，使用M7和M8来重建LRT，使用M2a和M8来鉴别受到正选择的位点。

icode用来更改所选序列的遗传密码子，以期得到更加准确的结果。

RateAncestor=1表明你想重建原始序列，如果RateAncestor=0说明你将避免这个计算。

不过使用效果并不明显，还需要进一步研究如何使用。

PAML使用中最重要的就是模型的选择：

PAML中所有的模型都在baseml和codeml这两个程序中使用。

这两个程序是最大似然程序，它们使用数值优化算法来最大化对数似然值。

这些模型最大的用途就是适用likelihoodratiotest（似然比率检验）来检测有趣的生物学假设。

这些模型是在Baseml中使用的，软件中常用的数学模型有：

JC69（JukesandCantor1969），K80（Kimura1980），F81（Felsenstein1981），F84（Felsenstein1984），HKY85（Hasegawa1984，1985），Tamura（1992），TamuraandNei（1993），andREV，alsoknowasGTRforgeneral-time-reversible（Yang1994）。

模型的一般遵循以下假设：

1.在不同的谱系中替换是独立发生的。

2.在不同的位点中替换也是独立发生的。

3.替换的过程我们通过时间均匀马尔科夫过程（time-homogeneousMarkovprocess）。

常用的两种检测方法：

1.Maximumlikelihoodestimates（MLEs）：

观测到的数据X的概率（probability），当做为一个未知参数θ的函数的时候，就叫做似然函数（likelihoodfunction）：

L（θ：

∣X）=f（θ∣X）。

根据似然规则（likelihoodprinciple），似然函数包括数据中关于参数θ所有的信息。

参数θ的最佳点估计（optimalpointestimate）可以通过最大化似然L的θ值或l（θ；X）的似然对数进行估计。

并且，似然曲线可以为未确定的点估计提供信息。

2.Likelihoodratiotests（LRTs）：

假设一个简单模型或无效模型（simplerornullmodel）有一个参数p0，更通用的模型或可选择的模型（generaloralternativemodel）有一个参数p1，两个模型的最佳似然值分别为l0和l1。

那么对数似然值差异（loglikelihooddifference）的两倍是：

2△l=2（l1-l0），如果无效模型（nullmodel）成立的话，那么对数似然值差异的二倍将与自由度是d.f.=p1-p0的卡方分布具有渐进关系（asymptotically）。

因此，对数似然值差异的二倍的检验统计可以通过比较卡方分布来检验无效模型（nullmodel）是否拒绝备择模型（alternativemodel）。

所谓Likelihoodratiotest（似然比率检测）是用来检验两个模型的。

离散伽玛模型（discrete-gammamodel）允许不同位点具有不同的变化速率。

Baseml中有核苷酸替换模型，Codeml中有不同位点替换速率变化的模型。

1.整理在Codeml中进行比较的两个模型比较有：

M1a（NearlyNeutral）和M2a（PositiveSelection）；M7（beta）和M8（beta&ω）。

2.整理认为M3对于正选择的LRT检测并不是十分适合，并不推荐适用M3模型。

使用似然比率检测可以验证正选择（Testingpositiveselectionusingthelikelihoodratiotest）。

整理推荐使用二到三种LRT来验证正选择。

第一个检测是比较M1a和M2a，

第二个检测是比较M7和M8。

Gamma分布中形状参数所表示的含义：

1.α＞1，大多数位点的替换速率在1附近，但有少数位点具有比较高或比较低的替换速率。

曲线形状为bell-shaped.

2.α→∞，表明所有的位点具有一个相同速率。

3.α≤1，表明大部分位点的替换速率比较低，或接近于不变，可是有一些位点具有比较高的替换速率。

曲线形状为L-shape。

PAML的一个重要功能就是检测基因是否受到正选择，即适应性选择。

但是现在用于估计适应性选择的方法，忽略了氨基酸的化学性质，这样得出的结果是不准确的，整理表示，直接通过dn和ds的比较来确定受到什么样的选择压力，是不准确的。

PAML中的无效模型（nullmodel）是指不允许任何位点的ω值大于1，如果nullmodel成立，则w小于1，基因受到负选择；如果nullmodel不成立，则w大于1，基因受到正选择。

研究表明，通过比较两个点模型，而得到的结果尤其可靠。

Ancestralreconstruction为探索数据提供了一个直观的方法，他被用于大量的数据分析，例如，评估不同谱系中的选择压力。

但是由于这种方法的简单和直观，会产生很多错误。

大多数重建原始序列的工作都忽略了这样一个事实，即使用假数据（pseudo-data）代替真实观察到的数据（realobserveddata），并且仅仅使用处于最佳状态（optimalcharacterstates）的一些特征，而忽略未处于最佳状态（suboptimalstates）的一些特征，从而产生一些系统上的偏差（systematicbiases）。

如果数据中的ds区域饱和的话，那么会导致我们低估ds，从而使dn/ds的值偏高，即ω的值偏高。

进化距离估计中有关序列间隔的处理进化距离估计中，排列时的间隔导致了某些复杂的问题。

同时，由于实验上的原因，也可能出现丧失信息的位点。

在距离估计中，一般忽略这些位点，可用两种不同的方法来进行处理。

一种方法是从数据分析中删除这些位点，称为完全删除（completedeletion）。

一般来说，这种方法较好，因为DNA或氨基酸序列的不同区段往往具有不同的演变规律。

然而在所研究的序列中，间隔不大或者是随机分布的，则可计算每个配对序列间的距离，并只忽略两个配对序列间的那些间隔。

这种方法称为成对删除（pairwisedeletion）。

这个过程在PAML中可以通过cleandata=0或1。

使用PAML进行数据分析的时候，所选用的序列越多，则得到的结果越可靠。

一般应该大于17条序列，此外影响LRT检验结果好坏的因素还有：

序列长度（sequencelength）、序列分异度（sequencedivergence）和正选择的强度（thestrengthofpositiveselection）（Anisimova,Bielawskietal.2001;Opazo,Palmaetal.2005）。

还有使用PAML进行序列分析的时候必须检验序列是否发生过重排，可以使用的软件由PLATO2.0，HYPHY等软件（Pond,Frostetal.2005;KosakovskyPond,Posadaetal.2006;Petersen,Bollbacketal.2007）。

如果是检测序列受到的选择压力的话，那么序列大于50codons即可，但是这些序列对于系统进化分析可能得到的结果不可靠（Pie2006）。

物种的基因只有在受到正选择作用才能不断的适应环境的变化，所以正选择在物种进化中起到了非常重要的作用（VallenderandLahn2004）。

PAML中常用的模型有以下几个：

1.M0，所有的谱系具有相同的ω0值。

2.M3，discrete，它的位点具有三个离散类（discreteclasses），并且具有不同的ω0。

3.M1a，nearlyneutralmodel，允许两个位点分类，0＜ω0＜1或ω＞1。

4.M2a，selection，具有一个额外的位点分类ω＞1。

5.M7，beta，ω＜1，具有额外10个位点分类。

6.M8，ω＞1，具有额外11个位点分类。

其中可以用于检验正选择的模型是M2a，M3和M8，得到了结果之后，我们可以再使用NEB（na?

veempirical）和BEB（bayesempirical）进行验证，通过后验概率，一般应该大于0.95。

PAML软件的模型中需要考虑的问题是：

1.遗传密码子的结构。

2.转换和颠换的比值。

3.密码子不同位置的速率。

PAML软件中的序列处理：

1.如果序列之间的分异度比较明显的话，那么需要4-5条序列；如果有10条序列的话，效果会比较好；如果序列数大于20，那么得到的结果就会比较可靠。

当然这也和序列之间的分异度有比较密切的了解。

2.最佳序列分异取决于序列的数目，如果树比较大的话，则可以容忍更多的变异。

一般情况下，如果总共的ds的距离大于0.5，那么就认为这种方法是合理的。

3.一般情况下，软件可以鉴别出一到两个位点受到较强的选择压力，但是有可能一些位点受到了选择压力，但比较弱，这是LRT方法会告诉你，这个位点存在选择压力，但是对于鉴别出，比较麻烦。

4.一般情况下，使用简单模型和复杂模型得出的结果应该是一致的。

所以，M0作为比较简单的模型，使用它得出的枝长、K以及w值应该和其他复杂模型得出的结果是一致的。

5.如果比对的序列是高度相似或高度分异的话，我们应该执行程序两遍。

6.如果使用NSsites执行多个模型，ncatG的参数也需要重新设置。

7.在计算过程中，序列的饱和性并不是一个主要的问题，一般较高的序列分异会带来更多的问题，不同的序列之中会具有不同的密码子使用偏好性和核苷酸组成。

8.计算出来的w值是不能为负值的。

9.现在已经找到，lnL的值在rst文件中，但是每个位点的w值，在什么地方呢？

即使找到了相关的参数，那么接下来再如何分析呢？

10.在进行序列分析的时候，要把编码序列末端的终止密码子去掉，以防止出现误差。

PAML中用到的进化树的格式：

Atreestructurefileisusedwhenrunmode=0or1.Thefilenameisspecifiedintheappropriatecontrolfile.Thetreetopologyistypicallyspecifiedusingtheparenthesisnotation（括号标序法）,althoughitispossibletouseabranchrepresentation,asdescribedbelow.Parenthesisnotation.Thefirstisthefamiliarparenthesisrepresentation,usedinmostphylogeneticsoftware.Thespeciescanberepresentedusingeithertheirnamesortheirindexescorrespondingtotheorderoftheiroccurrencesinthesequencedatafile.Ifspeciesnamesareused,theyhavetomatchexactlythoseinthesequencedatafile（includingspacesorstrangecharacters）.Branchlengthsareallowed.Thefollowingisapossibletreestructurefileforadatasetoffourspecies（human,chimpanzee,gorilla,andorangutan,occurringinthisorderindatafile）.Thefirsttreeisastartree,whilethenextfourtreesarethesame.Ifthetreehasbranchlengths,basemlandcodemlallowyoutousethebranchlengthsinthetreeasstartingvaluesformaximumlikelihooditeration.Whetheryoushoulduserootedorunrootedtreesdependsonthemodel,forexample,onwhetheramolecularclockisassumed.Withouttheclock（clock=0）,unrootedtreesshouldbeused,suchas（（1,2）,3,4）or（1,2,（3,4））.Withtheclockorlocal-clockmodels,thetreesshouldberootedandthesetwotreesaredifferentandbotharedifferentfrom（（（1,2）,3）,4）.InPAML,arootedtreehasabifurcationattheroot,whileanunrootedtreehasatrifurcationormultifurcationattheroot.TreefilesproducedbyPAUPandMacClade.PAMLprogramshaveonlylimitedcompatibilitywiththetreefilegeneratedbyPAUPorMacClade.Firstthe“[&U]”notationforspecifyinganunrootedtreeisig