用MrBayes构建系统发育树操作步骤整理总结

最近在学习构建进化树,尝试了用PHYML构建ML树和用MrBayes构建贝叶斯树,因为构建ML树时用的是网站的PHYML进行的,操作简便,所以就未做整理,只详细整理了一下贝叶斯树的构建步骤及相关的知识。算是学习笔记,有错误和不足之处或其它我还没有get到的功能还望道友指明并告知。

一、比对序列格式转换

将fasta格式的序列先在MEGA软件里进行比对,然后导出比对结果为PAUP format(图1),用记事本打开导出的结果,将里面的“missing=?”删掉,将里面的"datatype=nucleotide“改为"datatype=dna”,如果是蛋白序列就是“datatype=protein”。否则在用MrBayes程序打开时会报错。修改好后应是如下所示的样子:

图1

至此就获得了可用的“.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所示)

图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(数字指要指定物种在序列中的位置)。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 203,324评论 5 476
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,303评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 150,192评论 0 337
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,555评论 1 273
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,569评论 5 365
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,566评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,927评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,583评论 0 257
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,827评论 1 297
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,590评论 2 320
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,669评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,365评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,941评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,928评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,159评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,880评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,399评论 2 342

推荐阅读更多精彩内容

  • 刘小泽写于18.11.21聚类树是个有意思的东西,可以十分直观地看到每个物种的亲缘关系之前没有涉及到相关的分析,相...
    刘小泽阅读 4,582评论 0 8
  • 国家电网公司企业标准(Q/GDW)- 面向对象的用电信息数据交换协议 - 报批稿:20170802 前言: 排版 ...
    庭说阅读 10,849评论 6 13
  • 机器学习术语表 本术语表中列出了一般的机器学习术语和 TensorFlow 专用术语的定义。 A A/B 测试 (...
    yalesaleng阅读 1,957评论 0 11
  • 专业考题类型管理运行工作负责人一般作业考题内容选项A选项B选项C选项D选项E选项F正确答案 变电单选GYSZ本规程...
    小白兔去钓鱼阅读 8,970评论 0 13
  • 官网 中文版本 好的网站 Content-type: text/htmlBASH Section: User ...
    不排版阅读 4,362评论 0 5