PHYLIP使用.docx

1、PHYLIP使用 PHYLIP使用下面介绍几个软件的使用。首先是PHYLIP。其是多个软件的压缩包，下载后双击则自动解压。当你解压后就挥发现PHYLIP的功能极其强大，主要包括五个方面的功能软件：i，DNA和蛋白质序列数据的分析软件。ii，序列数据转变成距离数据后，对距离数据分析的软件。 iii，对基因频率和连续的元素分析的软件。iv，把序列的每个碱基/氨基酸独立看待（碱基/氨基酸只有0和1的状态）时，对序列进行分析的软件。v，按照DOLLO简约性算法对序列进行分析的软件。vi，绘制和修改进化树的软件。在此，我主要对前两种功能软件进行说明。 我们现在有几个序列如下：Mo3ATGTATTTCGT

2、ACATTACTGCCAGCCACCATGAATATTGCACGGTACCATMo5ATGTATTTCGTACATTACTGCCAGCCACCATGAATATTGTACGGTACCATMo6ATGTATTTCGTACATTACTGCCAGCCACCATGAATATTGTACGGTACCATMo7ATGTATTTCGTACATTACTGCCAGCCACCATGAATATTGTACAGTACCATMo8ATGTATTTCGTACATTACTGCCAGCCACCATGAATATTGTACAGTACCATMo9ATGTATCTCGTACATTACTGCCAGCCACCATGAATATTGTACGGTA

3、CCATMo12ATGTATTTCGTACATTACTG CCAGCCACCATGAATATTGTACGGTACCATMo13ATGTATCTCGTACATTACTGCCAGCCACCATGAATATTGTACGGTACCAT要对这8个序列进行进化树分析，按照上面的步骤，首先用CLUSTALX排列序列，输出格式为 *.PHY。用记事本打开如下图：图中的8和50分别表示8个序列和每个序列有50个碱基。然后，打开软件SEQBOOT，如下图：按路径输入刚才生成的 *.PHY文件，并在Random number seed (must be odd) ?的下面输入一个4N+1的数字后，屏幕显示如下：图中

4、的D、J、R、I、O、1、2代表可选择的选项，键入这些字母，程序的条件就会发生改变。D选项无须改变。J选项有三种条件可以选择，分别是Bootstrap、Jackknife和Permute。文章上面提到用Bootstraping法对进化树进行评估，所谓Bootstraping法就是从整个序列的碱基（氨基酸）中任意选取一半，剩下的一半序列随机补齐组成一个新的序列。这样，一个序列就可以变成了许多序列。一个多序列组也就可以变成许多个多序列组。根据某种算法（最大简约性法、最大可能性法、除权配对法或邻位相连法）每个多序列组都可以生成一个进化树。将生成的许多进化树进行比较，按照多数规则（majority-r

5、ule）我们就会得到一个最“逼真”的进化树。Jackknife则是另外一种随机选取序列的方法。它与Bootstrap法的区别是不将剩下的一半序列补齐，只生成一个缩短了一半的新序列。Permute是另外一种取样方法，其目的与Bootstrap和Jackknife法不同，这里不再介绍。R选项让使用者输入republicate的数目。所谓republicate就是用Bootstrap法生成的一个多序列组。根据多序列中所含的序列的数目的不同可以选取不同的republicate。当我们设置好条件后，键入Y按回车。得到一个文件outfile Outfile用记事本打开如下：这个文件包括了100个repub

6、licate。打开DNAPARS（最大简约性法）或DNAML（最大可能性法）软件。将刚才生成的outfile文件更名后输入。如下图：选项O是让使用者设定一个序列作为outgroup。一般选择一个亲缘关系与所分析序列组很接近的序列作为outgroup（本例子不选outgroup），outgroup选择的好坏将直接影响到最后的进化树的好坏。选项M是输入刚才设置的republicate的数目。设置好条件后，键入Y按回车。生成两个文件outfile和treefile。Outfile打开如下图：该文件包括了227个进化树。Treefile可以用TREEVIEW软件打开同样包含了这227个进化树。打开CO

7、NSENSE软件，将刚才生成的treefile文件更名后输入。如下图：键入Y按回车。生成两个文件outfile和treefile。Treefile用TREEVIEW打开，如下图：Outfile打开如下图：我们看出两个树是同样的。但在outfile的树上的数字表示该枝条的Bootstrap支持率（除以100.6）。到现在，8个序列的进化树分析（最大简约法）已经完成。如果要用邻位相连法对这8个序列进行分析的话，也首先执行SEQBOOT软件将这8个序列变成100个republicate。然后，打开DNADIST软件，把SEQBOOT生成的文件输入，如下图：选项D有四种距离模式可以选择，分别是Kimu

8、ra 2-parameter、Jin/Nei、Maximum-likelihood和Jukes-Cantor。选项T一般键入一个15-30之间的数字。选项M键入100。运行后生成文件如下图：这个文件包含了与输入文件相同的100个republicate，只不过每个republicate是以两两序列的进化距离来表示。文件中的每个republicate都省略了第一排的Mo3 Mo5Mo6Mo7Mo8Mo9Mo12Mo13。以这个输出文件为输入文件，执行NEIGHBOR软件。如下图：选项M键入100。生成两个文件outfile和treefile用记事本和TREEVIEW打开后，发现这两个文件都含有10

9、0个进化树。再将treefile文件更名后输入CONSENSE软件，又得到两个文件outfile和treefile，这就是最后的结果。以上是对DNA序列的分析，如果要对蛋白质序列进行分析，PROTDIST、PROTPARS等软件。其他软件的用法可以参照PHYLIP的documents。PAML软件的一些简单的具体的使用操作(2008-11-14 15:36:42) 转载标签： 杂谈1. 首先用Clustal X进行序列比对：要保证：保证核苷酸序列是三的倍数，没有终止密码子，核苷酸序列的第一位是密码子的第一位。假设序列名为cox1.fas2. 使用DAMBE软件进行转换成PML格式。软件使用截图

10、： 打开要换换的文件，然后“file” “save and convert sequence format”，在保存类型中选择“Yangs PAML”。那么此时的序列名为“cox1.PML”3 这样就可以得到文件“*.PML”，然后就直接把后缀改成“*.nuc”。那么此时的序列名为“cox1.nuc” 这样就完成了文件格式的转换。4 打开PAML软件的文件夹，找到文件名是“bin”的文件夹，打开之后，找到程序“codeml.exe”,把该程序复制到D盘的根目录下。（这一步并不是必要的，只是要把用到的几个程序放在同一个目录下）5 在你使用ClustalX进行序列比对的时候，会生成一棵进化树，适用

11、treeview软件可以打开，你需要的是把文件的后缀名改称“*.trees”。即树的文件名是“cox1.trees”，这就完成了树的格式的转换。6 然后再PAML4的文件夹中找到一个后缀是“*.ctl”的文件，把文件名改成“cox1.ctl”，复制到和“codeml.exe”相同的地方。7 要对codeml.ctl文件中的各个选项的值进行修改，具体内容如下：8 seqfile = cox1.nuc 按你自己的文件名进行修改，就可以了，9 treefile = cox1.treesoutfile = mlc * main result file name，noisy = 9 * 0,1,2,3,

12、9: how much rubbish on the screen，verbose = 0 * 0: concise; 1: detailed, 2: too muchrunmode = 0seqtype = 1 * 1:codons; 2:AAs; 3:codons-AAsCodonFreq = 2 * 0:1/61 each, 1:F1X4, 2:F3X4, 3:codon tableclock = 0aaDist = 0 * 0:equal, +:geometric; -:linear, 1-6:G1974,Miyata,c,p,v,aaaRatefile = wag.dat * onl

13、y used for aa seqs with model=empirical(_F) * dayhoff.dat, jones.dat, wag.dat, mtmam.dat, or your ownmodel = 0，这是使用的最简单的模型， * models for codons: * 0:one, 1:b, 2:2 or more dN/dS ratios for branches* models for AAs or codon-translated AAs:* 0:poisson, 1:proportional, 2:Empirical, 3:Empirical+F 28 * 6:

14、FromCodon, 7:AAClasses, 8:REVaa_0, 9:REVaa(nr=189)NSsites = 0 3 1 2 7 8 ，依次选取了6个模型。也可以选其中的两个，但必须是0和3，1和2，7和8。相互配对icode = 4 * 0:universal code; 1:mammalian mt; 2-10:see below如果是核基因的话就选0。fix_kappa = 0kappa = 5fix_omega = 0omega = 0.2getSE = 0RateAncestor = 0Small_Diff = .5e-6cleandata = 1 * remove sit

15、es with ambiguity data (1:yes, 0:no)?method = 0 * 0: simultaneous; 1: one branch at a time10 最后保证在同一个文件夹内同时具有：三个文件“codeml.exe”,“cox1.nuc”“cox1.trees”，这时候你双击codeml.exe，就可以运行程序。如果不能正确运行的话，你可以通过运行cmd，在dos情况下，运行codeml.exe，这样会有错误提示，知道你错在哪里了。常见命令解释： 1. Baseml.ctl的命令说明：2. noisy用来控制输出结果的多少，如果模型适用的运算比较多的话，no

16、isy的值可以选择的比较大，verbose可以控制结果文件中结果的多少。3. runmode = 0 表明在树的结构文件中估算树的拓扑结构。这个选项是我们通常情况下选择的，基本上可以满足我们的需要。4. Runmode = 1 or 2 表明通过星状-分解算法来进行启发式搜索树。Runmode = 2 这种算法是从星状树开始搜索，而runmode = 1 则表明软件读取多歧树是从树的结构文件中，并且同过比较去估计最佳二歧树。 5. runmode = 3 表明是逐步增加的。6. runmode = 4 通过简约法来搜索具有NNI perturbation的起始树。7. runmode = 5

17、表明从树的结构文件中来读取NNI perturbation with 起始树。8. Model 0，1，8 分别代表以下模型：JC69, K80, F81, F84, HKY85, T92, TN93, REV( also know as GTR), and UNREST。9. Mgene 用于和序列数据文件中的option G进行联合，用于多个基因和多个位点的联合分析。如果不使用option G的话，则选择0。10. ndata 用于指定文件中的分隔的数据集的数目。它的变化被用于模拟，你可以使用evolver来产生200个复制数据集，这是设置ndata = 200，然后用baseml进行分析

18、。11. clock 用于指定谱系之间速率恒定或变化的模型。Clock = 0，意味着整棵树中，不同分支之间不存在clock现象；Clock = 1，意味着global clock，所有的分支具有相同的进化速率；clock = 2 意味着local clock，所有分支之间的进化速率被分成几个部分；clock = 3意味着多个基因或多重分隔数据，允许分支的进化速率以不同的方式变化。；12. Codeml.ctl的使用说明：13. CodonFreq 用于平衡密码子替换模型中的密码子使用频率。Codonfreq = 0 说明每种密码子的使用频率是相同的；codonfreq = 1 说明是从平均核

19、苷酸频率中计算出来的；codonfreq = 2说明是从三个密码子位置的平均核苷酸频率得来的；codonfreq = 3则使用了三个参数。Codonfreq = 0，1，2和3 所代表的模型中使用的参数的数目分别为：0，3，9，和60。14. aadist用于指定氨基酸距离是否是相同的（= 0），还是属于Granthams matrix（= 1）。15. runmode = -2 执行ML方法来推测蛋白序列两两之间的dn和ds。16. model 用于估计各个分支之间的w值。 Model = 0，表明所有的谱系具有一个w比率（one w ratio）；model = 1，每一支具有一个速率（f

20、ree-ratio）；model = 2 表明速率的任意数字。17. NSsites 主要是用于指定模型允许dn/ds（w）在不同的位点之间变化。NSsites = m 表明对应于 model = m。变化的 ncatG被用来指定在一些特定的模型下的w分布的类型的数目。NcatG的值被用于执行一下分析：paper are 3 for M3 (discrete), 5 for M4 (freq), 10 for the continuous distributions (M5 gamma, M6: 2gamma, M7: beta, M8:beta&w, M9:beta&gamma, M10:

21、beta&gamma+1, M11:beta&normal1, and M12:0&2normal1, M13:3normal0). This means M8 will have 11 site classes (10 from the beta distribution plus 1 additional class) 。通过NSsites可以同时执行多个模型，例如：NSsites = 0 1 2 3 7 8，的意思就是同时执行M0，M1，M2a，M3，M7，和M8。作者建议：使用M1a和M2a来重建LRT，使用M7和M8来重建LRT，使用M2a和M8来鉴别受到正选择的位点。18. ico

22、de 用来更改所选序列的遗传密码子，以期得到更加准确的结果。19. RateAncestor = 1 表明你想重建原始序列，如果 RateAncestor = 0 说明你将避免这个计算。不过使用效果并不明显，还需要进一步研究如何使用。 PAML使用中最重要的就是模型的选择： PAML中所有的模型都在baseml和codeml这两个程序中使用。这两个程序是最大似然程序，它们使用数值优化算法来最大化对数似然值。这些模型最大的用途就是适用likelihood ratio test（似然比率检验）来检测有趣的生物学假设。这些模型是在Baseml中使用的，软件中常用的数学模型有：JC69（Jukes a

23、nd Cantor 1969），K80（Kimura 1980），F81（Felsenstein 1981），F84（Felsenstein 1984），HKY85（Hasegawa 1984，1985），Tamura（1992），Tamura and Nei（1993），and REV，also know as GTR for general-time-reversible（Yang 1994）。 模型的一般遵循以下假设： 1. 在不同的谱系中替换是独立发生的。 2. 在不同的位点中替换也是独立发生的。 3. 替换的过程我们通过时间均匀马尔科夫过程（time-homogeneous Mark

24、ov process）。常用的两种检测方法：1. Maximum likelihood estimates（MLEs）：观测到的数据X的概率（probability），当做为一个未知参数的函数的时候，就叫做似然函数（likelihood function）：L (：X) = f（X）。根据似然规则（likelihood principle），似然函数包括数据中关于参数所有的信息。参数的最佳点估计（optimal point estimate）可以通过最大化似然L的值或l（；X）的似然对数进行估计。并且，似然曲线可以为未确定的点估计提供信息。2. Likelihood ratio tests （

25、LRTs）：假设一个简单模型或无效模型（simpler or null model）有一个参数 p0，更通用的模型或可选择的模型（general or alternative model）有一个参数p1，两个模型的最佳似然值分别为l0和l1。那么对数似然值差异（log likelihood difference）的两倍是：2l = 2（l1 - l0），如果无效模型（null model）成立的话，那么对数似然值差异的二倍将与自由度是d.f. = p1- p0的卡方分布具有渐进关系（asymptotically）。因此，对数似然值差异的二倍的检验统计可以通过比较卡方分布来检验无效模型（null

26、 model）是否拒绝备择模型（alternative model）。 所谓Likelihood ratio test（似然比率检测）是用来检验两个模型的。离散伽玛模型（discrete-gamma model）允许不同位点具有不同的变化速率。 Baseml中有核苷酸替换模型，Codeml中有不同位点替换速率变化的模型。1. 作者在Codeml中进行比较的两个模型比较有：M1a（Nearly Neutral）和M2a（Positive Selection）；M7（beta）和M8（beta&）。2. 作者认为M3对于正选择的LRT检测并不是十分适合，并不推荐适用M3模型。 使用似然比率检测可以

27、验证正选择（Testing positive selection using the likelihood ratio test）。作者推荐使用二到三种LRT来验证正选择。第一个检测是比较M1a和M2a，第二个检测是比较M7和M8。 Gamma分布中形状参数所表示的含义： 1. 1，大多数位点的替换速率在1附近，但有少数位点具有比较高或比较低的替换速率。曲线形状为 bell-shaped 2. ，表明所有的位点具有一个相同速率。 3. 1，表明大部分位点的替换速率比较低，或接近于不变，可是有一些位点具有比较高的替换速率。曲线形状为 L-shape。 PAML的一个重要功能就是检测基因是否受到正

28、选择，即适应性选择。但是现在用于估计适应性选择的方法，忽略了氨基酸的化学性质，这样得出的结果是不准确的，作者表示，直接通过dn和ds的比较来确定受到什么样的选择压力，是不准确的。 PAML中的无效模型（null model）是指不允许任何位点的值大于1。PAML中null model是不允许w值大于1，如果null model成立，则w小于1，基因受到负选择；如果null model不成立，则w大于1，基因受到正选择。研究表明，通过比较两个点模型，而得到的结果尤其可靠。 Ancestral reconstruction 为探索数据提供了一个直观的方法，他被用于大量的数据分析，例如，评估不同谱系

29、中的选择压力。但是由于这种方法的简单和直观，会产生很多错误。大多数重建原始序列的工作都忽略了这样一个事实，即使用假数据（pseudo-data）代替真实观察到的数据（real observed data），并且仅仅使用处于最佳状态（optimal character states）的一些特征，而忽略未处于最佳状态（suboptimal states）的一些特征，从而产生一些系统上的偏差（systematic biases）。如果数据中的ds区域饱和的话，那么会导致我们低估ds，从而使dn/ds的值偏高，即 的值偏高。进化距离估计中有关序列间隔的处理进化距离估计中，排列时的间隔导致了某些复杂的问

30、题。同时，由于实验上的原因，也可能出现丧失信息的位点。在距离估计中，一般忽略这些位点，可用两种不同的方法来进行处理。一种方法是从数据分析中删除这些位点，称为完全删除（complete deletion）。一般来说，这种方法较好，因为DNA或氨基酸序列的不同区段往往具有不同的演变规律。然而在所研究的序列中，间隔不大或者是随机分布的，则可计算每个配对序列间的距离，并只忽略两个配对序列间的那些间隔。这种方法称为成对删除（pairwise deletion）。这个过程在PAML中可以通过cleandata = 0 或1。使用PAML进行数据分析的时候，所选用的序列越多，则得到的结果越可靠。一般应该大于

31、17条序列，此外影响LRT检验结果好坏的因素还有：序列长度（sequence length）、序列分异度（sequence divergence）和正选择的强度(the strength of positive selection )(Anisimova, Bielawski et al. 2001; Opazo, Palma et al. 2005)。 还有使用PAML进行序列分析的时候必须检验序列是否发生过重排，可以使用的软件由PLATO 2.0，HYPHY等软件(Pond, Frost et al. 2005; Kosakovsky Pond, Posada et al. 2006; P

32、etersen, Bollback et al. 2007)。如果是检测序列受到的选择压力的话，那么序列大于50codons即可，但是这些序列对于系统进化分析可能得到的结果不可靠(Pie 2006)。 物种的基因只有在受到正选择作用才能不断的适应环境的变化，所以正选择在物种进化中起到了非常重要的作用(Vallender and Lahn 2004)。 PAML中常用的模型有以下几个： 1. M0，所有的谱系具有相同的0值。 2. M3，discrete，它的位点具有三个离散类（discrete classes），并且具有不同的0。 3. M1a，nearly neutral model，允许两个位点分类，0 0 1 或1。 4. M2a，selection，具有一个额外的位点分类1。 5. M7，be