MrBayes

MrBayes #

MrBayes是贝叶斯推断法(Bayesian Inference)构建系统发育的经典软件, 下载地址点击此处。 

前期数据的准备与处理

  • 先对数据进行align。可以使用的软件有MAFFT,ClustalW。
  • 对序列保守区的选择,可以手动,也可以不删除。可以使用的软件有GBlocks。 使用DAMBE的Xia et al.的方法对序列进行饱和度检验。
  • 将fas等格式文件转换至nexus文件。可以手动,但是建议使用 SeqCombGo 这个软件转换,其它如DAMBE或者MEGA等写得过于繁琐。
  • 格式如下最好,名称必须是连续的,一般推荐用_来连接,"-" “:” “.“等不要出现在名称中,并且不要出现中文的各种符号,不然会报错。
  • 使用 ModelTest-NG 基于BIC数值对模型进行选择。BIC数值越低越适合。

参数解释 #

lset #

设置模型

lset applyto=(1) nucmodel=4by4 nst=2 rates=propinv coding=all;

nucmodel 指定 DNA 模型。4by4 标准情况 Doublet 核糖体 DNA 配对茎区 Codon 密码子编码区

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

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

coding character 采样方法。all 全部;variable 只对变化的特征进行采样。

prset #

prset applyto=(all) statefreqpr=dirichlet(1,1,1,1) ratepr=variable;

revmatpr 6个 GTR 概率矩阵的模型替换率。例如如果是 GTR,modeltes-ngSubst. Rates:可以获得 revmatpr=dirichlet(0.6646,13.5759,1.5451,9.1131,30.5330,1.0000)

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

shapepr GAMMA 速率变化分布的 shape parameter 。如果分布为 gammamodeltes-ngGamma shape:可以获得 shapepr=fixed(0.5606)

pinvarpr 不变位点的比例。modeltes-ngInv. sites prop: 可以获得 prset pinvarpr=fixed(0.4984)

ratepr 设置让分区以不同的 rates 进化

topologypr对于树的拓扑结构

brelenspr 对于树的枝长

applyto 设置适用范围,比如 applyto=(2,3,4,5)

MCMC #

MCMC 相关

mcmcp savebrlen=yes ngen=2000000 nchains=4 samplefreq=2000 diagnfreq=5000 burninfrac=0.50 append=yes;
mcmc;

mcmc 运行 MCMC

mcmcp MCMC 设置

ngen 运行 MCMC 代数

nchains 是 MCMCMC 的设置,4 表示 1 code chain + 3 heated chain

samplefreq 每多少代采样一次,当ngen大的时候,samplefreq变小

diagnfreq 每多少代 diagnostics test 一次

printfreq 打印在屏幕上的频率

burninfrac 丢弃 cold chain 的前百分之多少。

append 分析过早停止可以续上 ,如果append=yes

savebrlen 报错枝长

sum parameters #

sump burninfrac=0.50;

sump summarize 所有的 parameters 用 mcmc 一样的 burnin

sum trees #

sumt burninfrac=0.50 contype=halfcompat;

sumt summarize 所有的 trees 用 mcmc 一样的 burnin

Contype consensus 树类型 halfcompat majority rule consensus tee allcompat 所有兼容的进入这个树

unlink revmat=(all) pinvar=(all) shape=(all) statefreq=(all);

各个分区不要将这些关联在一起

模型调整 #

下述均是 prset 下的设置

GTR #

GTR 中 prset 下的 statefreqprrevmatpr 是必须要根据 jModelTest 或者 ModelTest-NG 设置的,例如statefreqpr=dirichlet(0.3022,0.1300,0.1544,0.4134) revmatpr=dirichlet(0.6646,13.5759,1.5451,9.1131,30.5330,1.0000)

JC #

必须设置 statefreqpr=fixed(equal)

Gamma (+G) #

必须根据 jModelTest 或者 ModelTest-NG 设置shapepr,例如shapepr=fixed(0.5606)

Propinv/Invgamma (+I/+I+G) #

必须根据 jModelTest 或者 ModelTest-NG 设置pinvarpr,例如 pinvarpr=fixed(0.4984)

运行 #

Linux 下命令行输入 mb ;Windows 下点击exe文件。输入exe 你的文件名。Windows 下如果你的文件和exe处于同一个文件夹,那么直接输入即可,如果不是,就需要输入路径。 Linux 下如果该目录下使用 mb 可以直接不输入路径,只输入文件名字。

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

want to stop the run (y/n)?出现,输入y即可。

输出文件 #

假设你的输入文件名为 f

  • f.ckp checkpoint file,如果上次跑的不够,这次mcmc append=yes 就可以利用 ckp 文件继续跑
  • f.con.tre consensus tree 合意树
  • f.lstat marginal likelihoods 边际似然值输出
  • f.mcmc MCMC-diagnostics test
  • f.mstat model parameters 模型参数输出
  • f.parts key to taxon bipartitions 样品分组关键
  • f.pstat parameter table 参数表输出
  • f.run1/2.p substitution model parameters 替换模型参数
  • f.run1/2.t topology and branch lengths 带有拓扑结构和树长的树
  • f.trprobs MCMC 过程中搜索到的树,按照 prob 排列
  • f.tstat summary statistics for informative taxon bipartitions 样品分组信息统计
  • f.vstatsummary statistics for branch and node parameters 分支和参数统计

选择f.con.tre打开,进行树的美化。

单基因建树示例 #

下面是 MrBayes 的进行 BI 分析的 bolck,只需要将下述 block 在 nex/nexus 文件中粘贴到数据 block 后面即可。可以调整 nchains samplefreqprset 必须调整,根据模型调整

begin mrbayes;
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;
mcmc;
sump burninfrac=0.50;
sumt burninfrac=0.50 contype=halfcompat;
end;

多基因建树示例 #

下面是 MrBayes 的进行 BI 分析的 bolck ,只需要将下述 block 粘贴到数据 block 后面即可。 可以调整 nchains samplefreqprset除了applyto=(all) 必须调整,根据模型调整partitioncharset根据自己的数据修改。

begin mrbayes;
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=fixed(0.5606);
mcmcp savebrlen=yes ngen=2000000 nchains=4 samplefreq=2000 diagnfreq=5000 burninfrac=0.50;
mcmc;
sump burninfrac=0.50;
sumt burninfrac=0.50 contype=halfcompat;
end;

形态学数据建树示例 #

目前 MrBayes 仅仅支持简单的 MK 模型类似与基因中的 JC 模型,形态学建树更多采取 MP 方法,详见 TNT 或者 PAUP*。不过,BI 法在形态中也是尤其优越性,可以作为一个比对结果。

与 PAUP* 相比,MrBayes 的 datatype 为 standrd 所能定义的状态要少。下面是一个形态学的实例。

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;

其后面的 MrBayes 的 Block 如下。可以调整 nchains samplefreqcoding 可以更改为 all, rates 也可以根据实际更改。

begin mrbayes;
lset rates = equal coding = variable;
mcmcp savebrlen=yes ngen=2000000 nchains=4 samplefreq=2000 diagnfreq=5000 burninfrac=0.50;
mcmc;
sump burninfrac=0.50;
sumt burninfrac=0.50 contype=halfcompat;
end;

形态学和分子数据建树示例 #

本人不建议使用二者的结合,因为随着序列的增长,形态学数据占据的比例会急剧下降,这对于形态学数据来说是一个不公的处理。

下面是数据示例。

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;

下面是 MrBayes 的 Block 。可以调整 nchains samplefreqcoding 可以更改为 all, rates 也可以根据实际更改,prset除了applyto=(all) 必须调整,根据模型调整partitioncharset根据自己的数据修改。

begin mrbayes;
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;
mcmc;
sump burninfrac=0.50;
sumt burninfrac=0.50 contype=halfcompat;
end;