Beast系统发育关系和分化时间的估算(上)

以下内容参考了同研究组师弟整理的笔记、和网络资源:Beast v1.8简易手册-修订版 · 语雀 (yuque.com),根据自己的实际操作整理。

Beast软件是采用贝叶斯演化分析的原理,用来估算系统发育关系和分化时间的软件,使用前提是需要有化石证据或其他的先验信息进行时间校准。

整体流程用到的软件安装都很简单,就是一般的windows系统下安装软件的正常操作,不再赘述。

Beats官网 https://www.beast2.org/

需要的软件包括:

Beast:该软件包包括Beats程序、BEAUTi、TreeAnnotator和其他的应用程序

Tracer:该程序是用来探究Beast结果的,对连续型参数的分布进行量化和可视化,提供可以诊断的信息。

Figtree:用来呈现系统发育树,尤其是通过Beast获取的树。

下面用我自己的一个实例展示一下整个流程。

数据的介绍

我用到的数据是线粒体基因组中的13个蛋白质编码基因和2个rRNA基因(12S和16S)。

因为是多基因,又有蛋白编码基因(编码核苷酸的密码子对应的三个碱基的演化速率存在差异),所以我选择对数据进行分区估算最优核苷酸替换模型。在进行多序列比对时,对两个rRNA基因分别采用mafft alignment,而对13个蛋白质编码基因分别translation align,保证不破坏编码同一核苷酸的三个碱基。之后,将15个比对后的多重序列串联在一起(整个过程在Geneious里操作)。

通过查阅文献,我找到了这里用到的类群的化石证据(包括年代信息),以及前人研究发现的几个类群的分化时间。

除了查阅文献,可以去这两个网站寻找化石证据和分化时间信息:

化石查询网站:https://fossilcalibrations.org/

分化时间查询网站:http://timetree.org/

第零步 估算分区模型

这一步不是必要的一步,所以叫第零步。这一步里我用PartitionFinder这个软件估算分区和每个分区最优的核苷酸替换模型。具体操作可以参考我的这篇文章:

https://www.jianshu.com/p/1a0214ef79a7

第一步 用BEAUTi建立BEAST XML 输入文件

1 载入比对序列的.nex文件(load NEXUS format alignment)

这一步的目的是将我们的序列信息和模型参数信息整合在一起,建立一个输入BEAST进行运算的文件。

首先要把序列文件导出为.nex文件(Geneious里操作)。

打开.nex文件,将分区的信息加入进去,这些信息在PartitionFinder输出的结果文件analysis的best_schem.txt里:

best_scheme.txt文件中展示分区和每个分区的最有替换模型的信息
best_scheme.txt文件中可以直接加入进.nex文件的分区信息

我们打开导出的序列数据的.nex文件,将上述分区信息加入进去,整体的效果如下:

打开BEAUTi,上方工具栏:File-Import Data,导入.nex文件:

可以看到,是按照给出的分区信息被拆分成12个部分。

2 定义校准节点(defining the calibration node)

主要参数的设置选项

我们可以看到主窗口上有很多个选键,我们从左向右依次设置。

点击Taxa,这里用来设置分类单元子集,之后你就可以为设置的分类单元的共同最近祖先(MRCA,most recent common ancester)添加分化时间校准信息。

点击左下角的+,添加分类单元,双击taxon set一列的行名,可以修改名称。刚打开窗口,所有分类群都处于左边的excluded taxa一栏,选择属于该分类单元的分类群,点击右箭头,把它们导入右栏included taxa。

Taxa

我们看到在taxon set右边有一个Mono?复选框,如果你确定你设置的分类单元是一个单系类群,可以在复选框里打勾。在之后及进行的MCMC分析中,这一类群就会保持为单系类群。

Mono?

3 设置核苷酸替换模型(setting the substitution model)

因为我这里选择了对数据进行分区,因此需要对之前partitionfinder运算得到的12个subset分别进行模型的设置。点击Partitions,全选后,点击窗口上方Unlink subst. Models,

点击Sites,左面substitution model下面是所有的subset,一个个分别选中后,根据partitionfinder的结果文件内容,在右面设置模型参数:

这里要修改的参数是Substitution Model和Site Heterogeneity Model。很简单,如果最优模型是GTR+I+G+X,Substitution Model选择GTR,Site Heterogeneity Model选择Gamma+Invariant Sites。因为我之前已经对编码氨基酸的三个碱基进行了分区,这里就不用考虑partition into codon positions了,选择off。

4 设置分子钟模型(setting the clock model)

点击Clocks,clock type下拉菜单中选择uncorrelated relaxed clock (注:这里的选择根据个人数据情况)。

5 设置树先验(tree prior)

树先验用于设置种群或物种规模随时间变化的模型。在tree prior的下拉菜单中选择Yule Process。这是一个简单的物种形成模型,更适用于分析来自不同物种的序列。

6 设置先验(Priors)

Priors用于设置模型参数的先验。我们需要基于已知的化石证据或参考文献信息,为一些分化时间设置一个先验分布,点击prior一栏的Using Tree Prior,在弹出的对话框中,设置Prior Distribution为Normal。在mean中填55,stdev中填0.5,意思是最近共同祖先距今55百万年,标准差为0.5个百万年。这里需要根据自己的先验信息设置。

对于一些我们设置的分类群子集,我们可能没有先验信息,就不用设置,可以根据已有的信息进行估算。

7 设置MCMC选项(setting the MCMC options)

Length of chain: this is the number of steps the MCMC will make in the chain before finishing. How long this should be depends on the size of the data set, the complexity of the model and the quality of answer required.

length of chain 指的是MCMC分析运行的总代数,可以通过增加代数提高ESS(effective sample size)值,使参数收敛。默认是10,000,000。

一个粗略的计算方法是总代数=3000*序列条数的平方。这里有20条序列,就是1,200,000代

Log parameters every:for the log file, the value should be set relative to the total length of the chain. Sampling too often will result in very large files with little extra benefit in terms of the precision of the analysis. Sample too infrequently and the log file will not contain much information about the distribution of the parameters.

log parameters every 指的是采样的频率,导致太多无用或意义不大的信息产生,采样不够频繁又不能包含足够多的参数分布信息。一般要根据length of chain的数值来设置,要保证样本数量(样本数量=length of chain/log parameter every)至少有10,000个。

8 生成BEAST XML文件(generating the BEAST XML file)

点击右下角的Generate BEAST File,生成一个.XML文件,然后通过BEAST运算该文件。

打开beast,点击BEAST XML File选项的Choose File,选择生成的XML文件,点击Run。

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

推荐阅读更多精彩内容