[TOC]
Bayesian的系统推断相较于ML推断，Bayesian推断会考虑进先验概率，考虑进先前研究的结果。更为重要的是，Bayesian推断是后续做分化时间分析的重要前提，Bayesian推断可以根据化石或地理记录以及分子进化速率推断出物种之间的分化时间。
Preparation
软件：
- BEAST2分子进化；分子时间计算经典软件；https://www.beast2.org
- TracerBayesian贝叶斯树的后验计算；Tracer
- FigTree系统树的可视化软件；FigTree
利用BEAST2中的BEAUti, BEAST2, TreeAnnotator做Bayesian推断；其次会利用Tracer做辅助推断；作后利用BEAST2中的bModelTest包做自动化模型的选择。
BEAST2学习网址：
数据： 同样利用16s_filtered.nex和rag1_filtered.nex。假设这两个基因的历史与物种分化历史一致；即我们可以利用这两个串联基因的基因树gene tree代表其中41个物种的物种树species tree
一些基本术语：
- Markov Chain Monte Carlo(
马尔科夫链蒙特卡洛方法)简称MCMC，是一个抽样方法，用于解决难以直接抽样的分布的随机抽样模拟问题 - Bayesian Inference:
概率统计中的应用所观察到的现象对有关概率分布的主观判断（即先验概率）进行修正的标准方法。重点就是加入一个先验概率的模型推断。
利用BEAST2做Bayesian系统树推断
- 打开软件BEAUti，导入16S and RAG1比对后文件；根据前面对核基因RAG1发现三联CDS不同碱基的模型不同，因此对RAG1分开处理，点击左下角split后会分开三个序列。
- 点Link tree，表示我们对这两个基因做串联处理，理论上建树是更能代表物种树的分化情况。点Link Clock Models，代表串联的基因做相同的分子模型对待。某一个支系中会以相同的分子进化速率相对待。而实际上在进化上，物种的进化速率是由物种特定的代谢速率和迭代时间所决定的（例如寿命长，两代之间时间长的物种分子速率就较慢）
- 跳过Tip Dates选项，此选项是针对迭代速度快的物种，例如virus所设定。对于Site model位点模型根据前面PAUP计算得到的核酸替换模型，Subst Model选择GTR，
Gamma Category Count选择4，Substitution Rate点击estimate自动计算。选择所有四个序列，设置相同的sitemodel设置。 - 对于ClockModel， 选择relaxed clock log normal，此clock model是大多数情况所适用的模型。设置clock rate为系统估算estimate。
- 对于Priors先验模型来说，选择Birth Death Model。对于系统默认的Yule模型来说，其会假设物种进化历史过程中不存在物种灭绝时间。Birth DeathModel适用于大多数情况。Nevertheless，实际运用过程中发现这部分模型筛选对结果影响较小。另外需要指定外类群OUTGROUP。在此41个物种中，选择”Danioxxrerioxx(Zebrafish)”作为外类群，其余的即为要计算的类群，标为Euteleosteomorpha。
- 随后，对于标为Euteleosteomorpha和Zebrafish的分化时间，我们根据前人文献中确定的时间或者http://timetree.org/网站的时间进行先验分化时间的确定分为250Ma以前。选择Log Normal设置先验密度。M代表先验密度平均值设为10，S代表标准差SD设为0.5，选择Mean In Real Space，设置Offset为240；设置Use Originate；设置Set monophyletic
- 随后根据ML树的结果设置其它为Monophyletic的物种支系；设置的多会减少BEAST计算的时间。
- 对于MCMC设置来说，抽样是越多越好，100million以上，会跑相当长时间。
tracelog选项填输出文件的名字combined.trees，log every设置为10000次，表示每10000次输出日志文件。最后file里保存为combined.xml文件。最后BEAST里运行。
利用拓展包bModelTest做核酸模型的自动选择。
前面我们在Site Model里根据前面PAUP的模型计算选择了GTR模型及相关参数，而BEAST2中也同样可以利用相关软件做自动计算。
- 在file里的manage packages里可以下载相关的拓展包。下载bModelTest，随后即可在SiteModel里选择BEAST model test并应用于所有序列。其它设置相同。
- 在Linux下可根据BEAUti保存的xml配置文件直接输入命令行进行计算。
1 | beast combined_bmodeltest.xml |
利用Tracer估算MCMC先验迭代次数。
Bayesian推断中需要大量的MCMC迭代，可以利用Tracer去做迭代次数的预估，进而减少计算量（类似的有R中的coda）。
- 计算ESS(Effective sample sizes)，ESS值应该大于200。对Trace的图进行探索确定标准。
Summarize the posterior tree distribution统计后验分布
- 利用figtree打开BEAST生成的树文件trees，在此树中，支长代表时间，大部分支聚集在一起也是代表其之间分化的时间较短。并且这样一颗树也只是所有树中的一颗，并不能代表所有的
- 利用BEAST2中的TreeAnnotator进行树的整合，老化率（Burnin percentage）设置为10，根据Tracer plots中stationary after the first 10% of the chain。选择Maximum clade credibility tree。Node height选择Mean Heights。
- 结果在figtree中的node labels中选择posterior后验概率即为每个支的支持率。可以看到有些支系小于0.6是不可靠的。
- 选择Scale axis即代表时间尺度，可以看到各个支系的大致分化时间。选择Node Bars中height_95%_HPD(highest posterior density)即代表分化时间的95%置信区间。
结果显示African和 Neotrocial cichlid fishes的分化时间大致为~35million时间。但是这个结果仍然值得质疑，因为Neo类群仅包括2个物种accs(markers)，并且此系统树是基于single root node而非实际大量的化石证据。后续可以继续根据化石记录来进一步推断时间。