最近在学习构建进化树,尝试了用PHYML构建ML树和用MrBayes构建贝叶斯树,因为构建ML树时用的是网站的PHYML进行的,操作简便,所以就未做整理,只详细整理了一下贝叶斯树的构建步骤及相关的知识。算是学习笔记,有错误和不足之处或其它我还没有get到的功能还望道友指明并告知。
一、比对序列格式转换
将fasta格式的序列先在MEGA软件里进行比对,然后导出比对结果为PAUP format(图1),用记事本打开导出的结果,将里面的“missing=?”删掉,将里面的"datatype=nucleotide“改为"datatype=dna”,如果是蛋白序列就是“datatype=protein”。否则在用MrBayes程序打开时会报错。修改好后应是如下所示的样子:
至此就获得了可用的“.nexus”格式文件。注意:要将文件与MrBayes程序放在同一目录下,否则会无法打开。(nexus格式:ntax=类群数,nchar=序列长度,datatype=数据类型)
二、MrBayes软件下载安装
直接在其官网下载就可以了(http://nbisweden.github.io/MrBayes/download.html),注意根据自己的电脑系统及配置下载安装相应的版本。(似乎是不用安装的,解压后双击应用程序就可以直接用了)
三、建树
1、运行MrBayes,输入“EXE filename.nex”命令,打开第一步中准备好的“filename.nexus”格式文件。执行成功,执行窗口会输出文件的简单信息。
2、参数设置
lset 设置替换模型参数
输入“Lset nst=6 rates=invgamma”命令,该命令设置进化模型为with gamma-distributed rate variation across sites和a proportion of invariable sites的GTR模型。模型可根据需要更改,不过一般无须更改。
Nst:核酸替代模型,1=JC69,2=F81,6=GTR,可以尝试使用三个模型运行,选择最优的结果
Rates:指定序列上每个位点的替换速率,equal表示替换速率都是一致的,gamma表示 用 gamma来确定序列上的替换速率
prset 设置模型的先验信息,一般不确定就设置以下默认参数,依次输入以下命令:
“prset statefreqpr=fixed(0.2612, 0.2180,0.2563, 0.2645)” 设置GTR模型中核苷酸平衡频率的先验概率
“prset revmatpr=fixed(1.0000, 3.6695, 0.4065,0.4065, 9.2409, 1.0000)” 设置GTR模型中替换速率的先验分布
“prset shapepr=fixed(0.8448)” 设置速率分布的尺度
“prset pinvarpr=fixed(0.5703)” 设置不变位点的比例
mcmc 设置抽样信息,输入如下命令
“mcmc Ngen=10000,Samplefreq=10,Printfreq=100”
Ngen:设置总抽样次数
Samplefreq:设置抽样频率
Printfreq:设置打印频率
注意:此处需要根据界面返回的分裂频率分支频率的标准偏差值判断抽样次数是否足够,分裂频率分支频率的标准偏差需要小于0.01,否则需要增加抽样次数。即:如果分裂频率分支频率split frequencies的标准偏差standard deviation在100,000代generations以后低于0.01,则当程序询问:“Continue the analysis? (yes/no)”,回答no;如果高于0.01,回答yes并输入抽样次数,直到该值低于0.01。(如下图2所示)
总结参数结果,输入如下命令:
“sump burnin=250”(在此为1000个样品,即任何相当于你取样的25%的值)。
程序会输出一个关于样品(sample)的替代模型参数的总结表,包括mean,mode和95 % credibility interval of each parameter,要保证所有参数PSRF(the potential scale reduction factor)的值接近1.0,如果不接近,分析时间要延长。
总结系统树,输入如下命令:
“sumt burnin=250”,然后就等待结果出现。
程序会输出一个具有每一个分支的posterior probabilities的树以及一个具有平均枝长mean branch lengths的树。所有输出文件均保存在与MrBayes程序同一目录下。生成“ . t” 文件和“ . parts” 、 “ . con” 及“ . trprobs”文件各一个。其中“ . con” 文件就是含有两棵共有树 ( consensus trees)的文件。
参数设置知识补充:
1、 lset 设置的参数 --------使用help lset查看 lset设置的参数
参数Nucmodel:指核酸的类型;其可选项中4by4指不区分序列上的位点;condon指使用密码子模型(序列上每个位点的替换速率会依据密码子模型来推断);doublet通常用于具有协同进化效应的序列;一般情况下都可以使用4by4,编码序列最好使用condon。
参数code:指密码子编码的规律;其可选项中universal指通用密码子使用规律;metmt用于推测线粒体内的基因;叶绿体用mycoplasma。
参数ploidy:指物种是单倍体还是二倍体。
2、 prset 设置模型的先验信息----------使用help prset查看相关参数及其说明
参数tratiopr:指定转换和颠换的比例,可以使用fixed来指定也可以通过beta分布来模拟产生。
参数revmatpr:指定GTR模型里替换速率的先验分布。
参数aamodelpr:指定氨基酸替换模型中先验参数的分布。
参数statefreqpr:指定GTR模型中核苷酸平衡频率的先验概率。
参数shapepr:设置速率分布的尺度参数。
3、 mcmc 设置抽样信息----------使用help mcmc查看相关参数
参数ngen:指的是总抽样次数。
参数nruns:制定独立分析的次数。如果为2,表明程序从两个独立的树形开始抽样,分析完成后综合两个分析结果。
参数nchain:设置每次分析时运行的chain的数量。
参数samplefreq:指定从总的样本数中抽样的频率。这个一般和参数ngen配合使用,以保证最后用于分析的样本量足够。举例:samplefreq设置为100,000,ngen设置为1000,这样100,000个随机样本中每隔1000个抽出一个样本,最后一共可以得到1000个样本。
参数burninfrac:该参数控制用于分析的样本的数量。在mcmc抽样初期的数据往往不可靠,需要去掉。burninfrac控制去掉的比例,如为0.25则表示样本的前25%的数据被去除,因此最终用来分析的总样本数为1000*(1-0.25)=750。
运行之后看方差,当方差<0.01时即可认为抽样达到平稳,不需要进行更多的抽样分析,否则应继续增加抽样次数以达到平稳为止。然后在屏幕输出结果中找到chain swap information,当其显示的四条链之间的交换频率在0.1~0.8之间可以认为结果是合理的,可以进行下一步分析,否则需要重新设置参数,包括足够长的ngen,适当降低temp等。
MrBayes的高级功能:
1、可以直接在序列文件后输入相关设置参数,然后让程序自动运行就可以了,但不建议将“sump”和“sumt”这两个命令写入,因为这两个命令具有诊断的作用。
2、partition功能,如果所分析的序列不均一,如既有编码区又有非编码区,或者想把编码区分为密码子第一、第二和第三位碱基单独分析的话就需要使用该功能。需要在序列文件中增加如下内容:
begin mrbayes;
charset pos1=1./3;
charset pos1=1./3;
charset pos1=1./3;
partition region=3:pos1,pos2,pos3;
Iset applyto=(all) nst=6 rates=gamma;
prset applyto=(all) ratepr=variable;
end;
charset用于设置变量并赋值;1./3指从第一个位开始,每隔三个位点取出一个值,并把这些值用变量pos1表示,代表密码子的第一位。其它类推。partition一行是告诉程序序列分为三部分,prset一行是用来指定三个部分的参数是独立估计的。
如果序列分为编码区和非编码区可以按一下写法:
charset noncoding=1-200;
charset coding=200-600;
partition region=2: noncoding, coding;
prset ratepr=variable;
3、指定外群,在一组序列中可以指定外群,如果不指定,则以序列文件的第一个物种作为外群。
指定外群设置命令为:outgroup 5或者outgroup my_taxon 5(数字指要指定物种在序列中的位置)。