前言：Mrbayes构建贝叶斯系统发育树的具体步骤在闵农高芳銮教授的课件（见文尾）已经说得较为详细，高老师不仅学识渊博，更愿意将自己所学分享给他人，这种可贵的品质让人敬佩。这里依据高老师的课件和软件说明书标记点重要知识。
参考文献：《MrBayes version 3.2 Manual: Tutorials and Model Summaries》
Manual_MrBayes_v3.2（说明书）.pdf

1. 常用脚本模板

BEGIN MRBAYES;
outgroup name1 #外群名1
outgroup name2 #外群名2
lset nst=6 rates=invgamma ngammacat=4; #进化模型参数
Prset statefreqpr=dirichlet(1,1,1,1); #进化模型参数
mcmc ngen=2000000 printfreq=1000 nruns=2 diagnfreq=5000 samplefreq=100 nchains=4 temp=0.1 burninfrac=0.25 checkpoint=yes savebrlens=yes;
sump burnin=5000;
sumt burnin=5000 contype=allcompat showtreeprobs=yes;
END;

上述模板代码如下：（复制可用）

  1. BEGIN MRBAYES;
  2. outgroup  name1
  3. outgroup  name2
  4. lset nst=6 rates=invgamma ngammacat=4;
  5. Prset statefreqpr=dirichlet(1,1,1,1);
  6. mcmc ngen=2000000 printfreq=1000 nruns=2 diagnfreq=5000 samplefreq=100 nchains=4 temp=0.1 burninfrac=0.25 checkpoint=yes savebrlens=yes;
  7. sump burnin=5000;
  8. sumt burnin=5000 contype=allcompat showtreeprobs=yes;
  9. END;

（1）参数说明

①红色部分根据自己的模型或者实际情况去调整。设置外群可以设置多个，也可以不设置，本示例设置了两个。本文的模型示例是GTR+G+I模型。
舍弃样本比例（burninfrac）一般是设置0.25或者0.1。
contype是生成树的类型，allcompat表示“adds all compatible groups to such a tree”，如果设置成contype=Halfcompat，’Halfcompat’ results in a 50% majority rule tree。

模型与代码参考表：
核苷酸替换模型速查表(for Mrbayes block)_修正.docx

②关于 ngammacat，即number of gamma categories，有些建树软件里简写作Ncat。这个值的默认值是4。关于这个值的界定有点难，参考了高老师的课件里提到的方法。
Mrbayes使用要点 - 图1
（图1）
Mrbayes使用要点 - 图2
（图2,点击图片放大更清晰）

Mrbayes使用要点 - 图3 （图3）
计算方法：（from 高老师）
（i）如果得到模型+G或（+G+I），在替换选择选择的结果中就有一个alpha值（或叫shape值），这个值就与+G对应，记下这个值备用（如果模型后没有+G或+G+I），那么那个选项是unform rate，说明不用考虑nCat值这个问题；

（ii）根据这个alpha值（或叫shape值）确定已绘制的参考曲线（图2左下角），看看最接近哪个，比如alpha值等0.43，那么最接受图2左下角红色那条Alpha线，再看该Alpha线与横坐标的交点值，该交点值的大小即为nCat值的大小。比如alpha值等0.22，这是介于红和黑的Alpha线之间，那么把ngammacat设置为4没问题，如果大致设置为2或3也够了，nCat =4 是默认的。

（iii） 如果是GTR或HKY后面不带G或G+I，那个Rate就是uniform ，就没有上面要设定nCat的问题了。

（iv）图3说的是(ii)的原理，即根据等分法，确定曲线中的x轴上的数值，简单理解：大约需要多少个长方形把这个曲线的面积大致求出来，这个个数就是对应的nCat数。
（v）Mrbayes运算的时候，有时候数据不收敛，可以由默认值的4改为 nCat=16，但是运算量会加大。在一些找最佳进化模型的软件里，得到的结果，有时候将GTR+G(nCat=4)模型就写成下图这样的：
Mrbayes使用要点 - 图4
③如果不输入diagnfreq(诊断频率），默认值是5000，我们一般是不输入diagnfreq，放在这里是为了让大家知道有这个东西。
在Prset里还可以设置两个参数shape和pinv，shap值只有模型为+G或+G+I时模型软件才会给出，pinv这个值只有模型为+I时模型软件才会给出。在mrbayes的原版说明书里提到了过度参数化，mrbayes本身就可以计算GAMMA的shape值（或叫alpha值）和pinv值，所以我们可以不指定shape值和pinv值，让mrbayes软件自己去估算它们。

④紫色部分基本可以保存不变。showtreeprobs=yes 这个要写进脚本里，最后的树枝上才会有后验概率数值。temp默认值是0.1。nruns=2表示执行两次独立的搜索，默认值也是2。 sumt，summary tree的缩写，可以设置多个子参数。sump，summary probability的缩写，用于估算后验概率。

⑤黑色部分也是需要调整的，总代数（ngen）、多少代生成一个样本（samplefreq）、舍弃样本数是变化的，公式：burnin= ngen ÷ samplefreq × burninfrac
其中burnin表示老化值， mcmc链在收敛前存在一段波动较为明显的过程，对枝长以及后验概率的估算影响较大，需要舍弃。
数据量比较小的话，总代数设置为100万代或者200万代够用了。

⑥nchains值
该值一般设置是4，即三个热链一个冷链。设置为n，即有n-1个热链heated chains被使用。默认n＝4，表示bayes 会使用3 个热链和1 个”cold” chain。根据相关资料，heating 对于大于50 个类群（序列）的分析是很重要的。增加热链数量对于分析大的困难的数据集可能有帮助。但分析时间也会随着链的增加成比例增加。
当热链设置大于超过3时，即 nchains值大于4时， mcmc里的参数Nswaps要增大，默认值是1.
关于这里，有一段原文描述如下：

  1.  By default, MrBayes uses Metropolis coupling to improve the MCMC sampling of the targetdistribution. The Swapfreq, Nswaps, Nchains, and Temp settings together control the Metropolis coupling behavior. When Nchains is set to 1, no heating is used. When Nchains is set to a value n larger than 1, then n-1 heated chains are used. By default, Nchains is set to 4, meaning that MrBayes will use 3 heated chains and one "cold" chain.
  2. In our experience, heating is essential for some data sets but it is not needed for others. Adding more than three heated chains may be helpful in analyzing large and difficult data sets. The time complexity of the analysis is directly proportional to the number of chains used (unless MrBayes runs out of physical RAM memory, in which case the analysis will suddenly become much slower), but the cold and heated chains can be distributed among processors in a cluster of computers and among cores in multicore processors using the MPI version of the program, greatly speeding up the calculations.

MrBayes uses an incremental heating scheme, in which chain i is heated by raising its posterior probability to the power 1/(1 + i”), where “ is the heating coefficient controlled by the Temp parameter. Every Swapfreq generation, two chains are picked at random and an attempt is made to swap their states. For many analyses, the default settings should work nicely. If you are running many more than three heated chains, however, you may want to increase the number of swaps (Nswaps) that are tried each time the chain stops for swapping. If the frequency of swapping between chains that are adjacent in temperature is low, you may want to decrease the Temp parameter.


特别说明！！！你设置的线程数要小于或等于你在mybayes脚本里设置的nchains数（马尔科夫链数），且两者为倍数关系，比如你设置了nchains=4，那么一次独立搜索最多可以开4个线程(如果执行两次独立搜索，最多可以开8线程）；如果设置nchains=8，执行一次独立搜索最多可以开8个线程。软件mybayes默认nchains值是4。
对于大的困难的数据集，nchains设置越大越好，但是nchains越大也会增加分析时间，同时，nchains增大的同时Nswaps也要相应增大，但是Nswaps却没有针对nchains值的参考的范围。综合考虑，还是设置nchains=4。

⑦MrBayes会保留了上一代的参数值并将使用那些作为下一次分析的起始值，除非使用重置值：

  1. mcmc starttrees=random startvals=reset


⑧mrbayes其他的参数默认就行，默认的参数相对而言是比较合适的且适合大多数的分析，除非你真正知道这个参数的含义。

（2）序列长度

如果用来建树的序列太长，mrbayes运行报错，需要将序列变成interleave形式(交错排列式）的连续.nex文件。

2. linux上安装多线程版本Mrbayes

（1）目前最新版本是mrbayes3.2.6版本(截止至2018年11月16日），先下载源码编译包（http://mrbayes.sourceforge.net/download.php，点击 source code下载），要MrBayes-3.2.6.tar.gz格式的。不要下载MrBayes-develop.zip，这个是单线程版本。
怎么编译安装多线程版本的mrbayes？将MrBayes-3.2.6.tar.gz解压后，有个src文件夹，里面有个名字为CompileInstructions.txt的文件，介绍了如何进行软件的安装，其中就介绍如何安装多线程版本。
请注意：安装多线程版本mrbayes需要mpi等环境，所以最好按照CompileInstructions.txt文件里的步骤，先把需要的环境安装了（文件里有详细名称），再次强调一定要在联网环境下弄这些。

（2）在网上有提到MrBayes-3.2.1版本安装多线程版后，用多线程版跑数据，树的分枝不对，需要打补丁。
故在使用目前最新的MrBayes 3.2.6版本前，我先用一个小数据做了测试（小数据是十几条1000bp左右DNA序 列，建树前先用jmodeltest找的模型），分别开单线程、4线程、8线程，Nswaps都取默认值1；发现三个不同线程得到的树的拓扑结构和关键节点的后验概率基本一致，其中单线程和8线程的拓扑结构更相似。另外再用另一建树软件phyML对这个用来测试的小数据DNA序列同样进行建树分析（和前面的用同一核苷酸模型），发现phyML生成的树和MrBayes 3.2.6的单线程、4线程、8线程建的树拓扑结构也基本一致。
以上表明，目前的MrBayes 3.2.6版本可能消除了之前3.3.1版本的bug，安装的多线程是可以用的。

3. 最后收敛里部分参数的意义

  1. （1）the total tree length：IL
  2. （2）the six reversible substitution rates：r(A<->C)，r(A<->G)，r(A<->T)，r(C<->G)，r(C<->T)，r(G<->T)
  3. （3）the four stationary state frequencies：pi(A), pi(C), pi(G), pi(T)
  4. （4）the shape of the gamma distribution of ratevariation across sites: alpha
  5. （5）the proportion of invariable sites: pinvar

Note： that the six rate parameters of the GTR model are given as proportions of the ratesum.

最后收敛里ESS值低怎么办？
增加抽样频率；整合多个独立链结果；增加链的长度。

该软件详细的使用方法，请参考高芳銮教授的课件。
附件：系统发育分析之BI篇（By Raindy）.pdf

若有错误之处，欢迎及时反馈。