查看“︁MrBayes”︁的源代码

MrBayes是贝叶斯推断法(Bayesian Inference)构建系统发育的经典软件，[https://nbisweden.github.io/MrBayes/download.html 下载地址点击此处。] 

前期数据的准备与处理
*先对数据进行align。可以使用的软件有MAFFT，ClustalW。
*对序列保守区的选择，可以手动，也可以不删除。可以使用的软件有GBlocks。 使用DAMBE的Xia et al.的方法对序列进行饱和度检验。
*将fas等格式文件转换至nexus文件。可以手动，但是建议使用 [https://github.com/MalacoLab/SeqCombGo SeqCombGo] 这个软件转换，其它如DAMBE或者MEGA等写得过于繁琐。
*格式如下最好，名称必须是连续的，一般推荐用_来连接，“-” “:” “.”等不要出现在名称中，并且不要出现中文的各种符号，不然会报错。
*使用 ModelTest-NG 基于BIC数值对模型进行选择。BIC数值越低越适合。
==参数解释==
===lset===
设置模型<pre>lset applyto=(1) nucmodel=4by4 nst=2 rates=propinv coding=all;</pre><code>nucmodel</code> 指定 DNA 模型。<code>4by4</code> 标准情况 <code>Doublet</code> 核糖体 DNA 配对茎区 <code>Codon</code> 密码子编码区

<code>nst</code> 替换模型。<code>1</code> JC/F81 <code>2</code> HKY <code>6</code> GTR

<code>rates</code> 分布。Equal/Gamma/LNorm/Propinv/Invgamma/Adgamma/Kmixture。一般就 +G gamma; +I +G invgamma; +I propinv; 无 equal。

<code>coding</code> character 采样方法。<code>all</code> 全部；<code>variable</code> 只对变化的特征进行采样。
===prset===
<pre>prset applyto=(all) statefreqpr=dirichlet(1,1,1,1) ratepr=variable;</pre><code>revmatpr</code> 6个 GTR 概率矩阵的模型替换率。例如如果是 GTR，<code>modeltes-ng</code> 的 <code>Subst. Rates:</code>可以获得 <code>revmatpr=dirichlet(0.6646,13.5759,1.5451,9.1131,30.5330,1.0000)</code>

<code>statefreqpr</code> GTR 概率矩阵的固定核苷酸频率。如模型为 JC/SYM 模型时，<code>statefreqpr=fixed(equal)</code>；例如如果是 GTR，<code>modeltes-ng</code> 的 <code>Frequencies:</code>可以获得<code>statefreqpr=dirichlet(0.3022,0.1300,0.1544,0.4134)</code>

<code>shapepr</code> GAMMA 速率变化分布的 shape parameter 。如果分布为 <code>gamma</code> 则 <code>modeltes-ng</code> 的 <code>Gamma shape:</code>可以获得 <code>shapepr=fixed(0.5606)</code>

<code>pinvarpr</code> 不变位点的比例。<code>modeltes-ng</code> 的<code>Inv. sites prop:</code> 可以获得 <code>prset pinvarpr=fixed(0.4984)</code>

<code>ratepr</code> 设置让分区以不同的 rates 进化

<code>topologypr</code>对于树的拓扑结构

<code>brelenspr</code> 对于树的枝长

<code>applyto</code> 设置适用范围，比如<code>applyto=(2,3,4,5)</code>
===MCMC===
MCMC 相关<pre>mcmcp savebrlen=yes ngen=2000000 nchains=4 samplefreq=2000 diagnfreq=5000 burninfrac=0.50 temp=0.05 append=yes;
mcmc;</pre><code>mcmc</code> 运行 MCMC

<code>mcmcp</code> MCMC 设置

<code>ngen</code> 运行 MCMC 代数

<code>nchains</code> 是 MCMCMC 的设置，4 表示 1 code chain + 3 heated chain

<code>samplefreq</code> 每多少代采样一次，当<code>ngen</code>大的时候，<code>samplefreq</code>变小

<code>diagnfreq</code> 每多少代 diagnostics test 一次

<code>printfreq</code> 打印在屏幕上的频率

<code>burninfrac</code> 丢弃 cold chain 的前百分之多少。

<code>temp</code> 可以稍短跑一段，看看 log 然后再调整温度。

<code>append</code> 分析过早停止可以续上 ，如果<code>append=yes</code>

<code>savebrlen</code> 报错枝长
===sum parameters===
<pre>sump burninfrac=0.50;</pre><code>sump</code> summarize 所有的 parameters 用 <code>mcmc</code> 一样的 <code>burnin</code>
===sum trees===
<pre>sumt burninfrac=0.50 contype=halfcompat;</pre><code>sumt</code> summarize 所有的 trees 用 <code>mcmc</code> 一样的 <code>burnin</code>

<code>Contype</code> consensus 树类型 <code>halfcompat</code> majority rule consensus tee <code>allcompat</code> 所有兼容的进入这个树
===unlink===
<pre>unlink revmat=(all) pinvar=(all) shape=(all) statefreq=(all);</pre>各个分区不要将这些关联在一起

=== log ===
 log start append filename=any.log;
<code>start</code> 开始记录 <code>append</code> 及时存在也在后面添加

==模型调整==
下述均是 <code>prset</code> 下的设置
===GTR===
GTR 中 <code>prset</code> 下的 <code>statefreqpr</code> 和 <code>revmatpr</code> 是必须要根据 jModelTest 或者 ModelTest-NG 设置的，例如<code>statefreqpr=dirichlet(0.3022,0.1300,0.1544,0.4134) revmatpr=dirichlet(0.6646,13.5759,1.5451,9.1131,30.5330,1.0000)</code>
===JC===
必须设置 <code>statefreqpr=fixed(equal)</code>
===Gamma (+G)===
必须根据 jModelTest 或者 ModelTest-NG 设置<code>shapepr</code>，例如<code>shapepr=exponential(0.5606)</code>
===Propinv/Invgamma (+I/+I+G)===
必须根据 jModelTest 或者 ModelTest-NG 设置<code>pinvarpr</code>，例如 <code>pinvarpr=uniform(0,0.4984)</code>
==运行==
Linux 下命令行输入 <code>mb</code> ；Windows 下点击exe文件。输入<code>exe 你的文件名</code>。Windows 下如果你的文件和exe处于同一个文件夹，那么直接输入即可，如果不是，就需要输入路径。 Linux 下如果该目录下使用 <code>mb</code> 可以直接不输入路径，只输入文件名字。

当 Average standard deviation of split frequencies 持续小于0.01时，那么，使用Ctrl+c结束，会有话语Do you really

want to stop the run (y/n)?出现，输入y即可。但是这样会缺少 log 所需要的判断输出。
==输出文件==
假设你的输入文件名为 <code>f</code>
*<code>f.ckp</code> checkpoint file，如果上次跑的不够，这次<code>mcmc append=yes</code> 就可以利用 ckp 文件继续跑
*<code>f.con.tre</code> consensus tree 合意树
*<code>f.lstat</code> marginal likelihoods 边际似然值输出
*<code>f.mcmc</code> MCMC-diagnostics test
*<code>f.mstat</code> model parameters 模型参数输出
*<code>f.parts</code> key to taxon bipartitions 样品分组关键
*<code>f.pstat</code> parameter table 参数表输出
*<code>f.run1/2.p</code> substitution model parameters 替换模型参数
*<code>f.run1/2.t</code> topology and branch lengths 带有拓扑结构和树长的树
*<code>f.trprobs</code> MCMC 过程中搜索到的树，按照 prob 排列
*<code>f.tstat</code> summary statistics for informative taxon bipartitions 样品分组信息统计
*<code>f.vstat</code>summary statistics for branch and node parameters 分支和参数统计
选择<code>f.con.tre</code>打开，进行树的美化。

== 判断 ==
tracer 查看 p1/2 总和的 ESS 是否大于 200

log 输出 <code>Chain swap information</code> 左上角必须在 0.1 到 0.7 之间，不然就调整 <code>temp</code>
       Chain swap information for run 1:
 
                    1       2       3       4 
            ----------------------------------
          1 |            0.64    0.36    0.18 
          2 |  166584            0.65    0.37 
          3 |  166540  166925            0.66 
          4 |  166842  166104  167005         
 
       Chain swap information for run 2:
 
                    1       2       3       4 
            ----------------------------------
          1 |            0.64    0.35    0.17 
          2 |  166901            0.65    0.36 
          3 |  167107  166028            0.65 
          4 |  166641  166567  166756         

== 单基因建树示例 ==

下面是 MrBayes 的进行 BI 分析的 bolck，只需要将下述 block 在 nex/nexus 文件中粘贴到数据 block 后面即可。可以调整 <code>nchains</code> <code>samplefreq</code>；<code>prset</code> 必须调整，根据'''模型调整'''。<pre>begin mrbayes;
log start append filename=any.log;
lset nucmodel=4by4 nst=6 rates=invgamma coding=all;
prset statefreqpr=dirichlet(0.3022,0.1300,0.1544,0.4134) revmatpr=dirichlet(0.6646,13.5759,1.5451,9.1131,30.5330,1.0000) shapepr=fixed(0.5606) pinvarpr=fixed(0.4984);
mcmcp savebrlen=yes ngen=2000000 nchains=4 samplefreq=2000 diagnfreq=5000 burninfrac=0.50 temp=0.05;
mcmc;
sump burninfrac=0.50;
sumt burninfrac=0.50 contype=halfcompat;
end;</pre>
==多基因建树示例==
下面是 MrBayes 的进行 BI 分析的 bolck ，只需要将下述 block 粘贴到数据 block 后面即可。 可以调整 <code>nchains</code> <code>samplefreq</code>；<code>prset</code>除了<code>applyto=(all)</code> 必须调整，根据'''模型调整'''，<code>partition</code> 和 <code>charset</code>根据自己的数据修改。<pre>begin mrbayes;
log start append filename=any.log; 
outgroup 5406.1; 
CHARSET ITS = 1-1027; 
CHARSET 16S = 1028-1391; 
partition gene=2:ITS,16S; 
set partition=gene; 
unlink revmat=(all) pinvar=(all) shape=(all) statefreq=(all);
prset applyto=(all) ratepr=variable;
lset applyto=(1) nucmodel=4by4 nst=2 rates=propinv coding=all;
prset applyto=(1) pinvarpr=fixed(0.4984);
lset applyto=(2) nucmodel=4by4 nst=1 rates=gamma; 
prset applyto=(2) statefreqpr=fixed(equal) shapepr=exponential(0.5606);
mcmcp savebrlen=yes ngen=2000000 nchains=4 samplefreq=2000 diagnfreq=5000 burninfrac=0.50 temp=0.05;
mcmc;
sump burninfrac=0.50;
sumt burninfrac=0.50 contype=halfcompat;
end;</pre>
==形态学数据建树示例==
目前 MrBayes 仅仅支持简单的 MK 模型类似与基因中的 JC 模型，形态学建树更多采取 MP 方法，详见 TNT 或者 PAUP*。不过，BI 法在形态中也是尤其优越性，可以作为一个比对结果。

与 PAUP* 相比，MrBayes 的 datatype 为 standrd 所能定义的状态要少。下面是一个形态学的实例。<pre>Begin data;
    Dimensions ntax=25 nchar=28;
    Format datatype=Standard gap=- missing=?;
    Matrix
Camaena              0000011200120222223322222220
Mastigeulota         0001000001012101100000000000
Metodontia           0001000000012101000200000000
Fruticicola          1211100000000111000210001000
Acusta               0001211200000111111210000100
Bradybaena           1121000000011101000200000000
Karaftohelix         1211111200000011000001000000
Cathaica             0001000000011110100000000000
Pliocathaica         0001011200000111100001000000
Aegista              0001011200113111000200000000
Plectotropis         0001011200113110000200000000
Pseudaspasita        0001011200111111000200000000
Platypetasus         0001000000000101100000000000
Stilpnodiscus_spa    0001100000000111100000000000
Stilpnodiscus_spb    0001100000000101100000000000
Laeocathaica         0001211200000111100000000000
Pseudobuliminus      0001000000000101100200000000
Trishoplita          0001011200100011112000000000
Euhadra              0001011200100011112000000000
Nesiohelix           0000011200100011111000010000
Aegistohadra         0000011210100111100200000010
Eueuhadra            0000011210100111100000100000
Calocochlea          1120012100000111100110000000
Pfeifferia           1120012100000111100110000001
Trichobradybaena     1121000001000101100200000000
    ;
End;</pre>其后面的 MrBayes 的 Block 如下。可以调整 <code>nchains</code> <code>samplefreq</code>；<code>coding</code> 可以更改为 <code>all</code>, <code>rates</code> 也可以根据实际更改。<pre>begin mrbayes;
log start append filename=any.log;
lset rates = equal coding = variable;
mcmcp savebrlen=yes ngen=2000000 nchains=4 samplefreq=2000 diagnfreq=5000 burninfrac=0.50 temp=0.05;
mcmc;
sump burninfrac=0.50;
sumt burninfrac=0.50 contype=halfcompat;
end;</pre>
==形态学和分子数据建树示例==
本人不建议使用二者的结合，因为随着序列的增长，形态学数据占据的比例会急剧下降，这对于形态学数据来说是一个不公的处理。

下面是数据示例。<pre>Begin data;
    Dimensions ntax=5 nchar=100;
    Format datatype=mixed(Standard:1-28,DNA:29-100) interleave=yes gap=- missing=?;
    Matrix
Camaena              0000011200120222223322222220CATTCGACCCGATGAGATTAGCAGTATCATTCGTAAACAAATTGAAGGATATGTTCCAGAAGTAAAGGTTGT
Mastigeulota         0001000001012101100000000000CATTCGACCCGATGAGATTAGCAGTATCATTCGTAAACAAATTGAAGGATATGTTCCAGAAGTAAAGGTTGT
Metodontia           0001000000012101000200000000CATTCGACCCGATGAGATTAGCAGTATCATTCGTAAACAAATTGAAGGATATGTTCCAGAAGTAAAGGTTGT
Fruticicola          1211100000000111000210001000CATTCGACCCGATGAGATTAGCAGTATCATTCGTAAACAAATTGAAGGATATGTTCCAGAAGTAAAGGTTGT
Acusta               0001211200000111111210000100CATTCGACCCGATGAGATTAGCAGTATCATTCGTAAACAAATTGAAGGATATGTTCCAGAAGTAAAGGTTGT
    ;
end;</pre>下面是 MrBayes 的 Block 。可以调整 <code>nchains</code> <code>samplefreq</code>；<code>coding</code> 可以更改为 <code>all</code>, <code>rates</code> 也可以根据实际更改，<code>prset</code>除了<code>applyto=(all)</code> 必须调整，根据'''模型调整'''，<code>partition</code> 和 <code>charset</code>根据自己的数据修改。<pre>begin mrbayes;
log start append filename=any.log;
charset morphology = 1-28;
charset gene = 29-100;
partition mixed=2:morphology,gene;
set partition=mixed;
unlink revmat=(all) pinvar=(all) shape=(all) statefreq=(all);
prset applyto=(all) ratepr=variable;
lset applyto=(1) rates=gamma; 
prset applyto=(1) shapepr=fixed(0.5606);
lset applyto=(2) nucmodel=4by4 nst=2 rates=propinv coding=all;
prset applyto=(2) pinvarpr=fixed(0.4984);
mcmcp savebrlen=yes ngen=2000000 nchains=4 samplefreq=2000 diagnfreq=5000 burninfrac=0.50 temp=0.05;
mcmc;
sump burninfrac=0.50;
sumt burninfrac=0.50 contype=halfcompat;
end;</pre>
[[分类:Bioinformatics]]