MalacoKnowledge

MrBayes

Malacology（留言 | 贡献）2023年1月1日 (日) 03:34的版本 （add temp related）
(差异) ←上一版本 | 最后版本 (差异) | 下一版本→ (差异)

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-ng 的 Subst. Rates:可以获得 revmatpr=dirichlet(0.6646,13.5759,1.5451,9.1131,30.5330,1.0000)

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

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

pinvarpr 不变位点的比例。modeltes-ng 的Inv. 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 temp=0.05 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 的前百分之多少。

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

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

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

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

log

log start append filename=any.log;

start 开始记录 append 及时存在也在后面添加

模型调整

下述均是 prset 下的设置

GTR

GTR 中 prset 下的 statefreqpr 和 revmatpr 是必须要根据 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=exponential(0.5606)

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

必须根据 jModelTest 或者 ModelTest-NG 设置pinvarpr，例如 pinvarpr=uniform(0,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即可。但是这样会缺少 log 所需要的判断输出。

输出文件

假设你的输入文件名为 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打开，进行树的美化。

判断

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

log 输出 Chain swap information 左上角必须在 0.1 到 0.7 之间，不然就调整 temp 

      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 后面即可。可以调整 nchains samplefreq；prset 必须调整，根据模型调整。

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;

多基因建树示例

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

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;

形态学数据建树示例

目前 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 samplefreq；coding 可以更改为 all, rates 也可以根据实际更改。

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;

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

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

下面是数据示例。

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 samplefreq；coding 可以更改为 all, rates 也可以根据实际更改，prset除了applyto=(all) 必须调整，根据模型调整，partition 和 charset根据自己的数据修改。

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;
检索自“https://knwl.malacology.net/index.php?title=MrBayes&oldid=91”