2008-1-5 05:14
bioin
MrBayes命令块的含义
MrBayes命令块的含义.谢谢!最近在用 MrBaye:pai.1构建进化树,下面是我用到的命令块.请高手指教一下下列命令的含义.另外,还想请问一下,如何设置命令块中ngen和burnin的数值?谢谢了!
begin mrbayes;
log start filename=拟南芥.log replace;
1set rates=gamma;
prset aamodelpr=mixed;
set autoclose=yes;
mcmcp ngen=10000 printfreq=1000 samplefreq=100\
nchains=4
savebrlens=yes filename=拟南芥;
mcmc;
sumt filename=拟南芥.t burnin=20 contype=halfcompat;
plot filename=拟南芥.p;
log stop;
end;
================================================================================
haha 看问题也是学习，我改天看看这个软件 ，
================================================================================
在这先谢谢mawencai2001了!
================================================================================
begin mrbayes; 开始mrbayes执行
log start filename=拟南芥.log replace; 记录日志文件（这一步可以略去）
1set rates=gamma; 设定likelihood参数，碱基替换率使用gamma参数
prset aamodelpr=mixed; 设定模型优先值
set autoclose=yes; 自动关闭
mcmcp ngen=10000 printfreq=1000 samplefreq=100\ 设定隐马氏模型，循环次数为10000，后面两个是取样和显示结果的参数，默认值就是这个，可以略去
nchains=4 进行4次隐马氏链计算（默认值是2次，一般默认即可）
savebrlens=yes filename=拟南芥; 设定需要保存的树文件
mcmc; 开始计算
sumt filename=拟南芥.t burnin=20 contype=halfcompat; 计算一致树，前2000(20×100)次结果被舍去。
plot filename=拟南芥.p; 计算距离矩阵？
log stop;
end;
MrBayes的手册是入门指引。详细命令文档请在MrBayes界面下使用manual命令来查看。或者用help+命令名来参考。附件是我保存下来的文档。
几点意见：
既然你用到了gamma值，为什么没有给出设定的数值？建议你先用modeltest做一下，选取适当距离模型和gamma值。
ngen太低，一般都是1000000。因为mcmc是渐进式逼近的，重复次数太低达不到最优结果。
同样，sumt的burnin值也太低。根据我的经验和文献报道，至少前5000个generation都是要burnin的，也就是burnin至少要设置到50以上。
最后提醒注意命令格式，每条命令一定要用分号结束。
[url=http://www.dxy.cn/bbs/user/download/9180732/commref_mb3.1.2.txt][color=#0000ff]commref_mb3.1.2.txt[/color][/url]
(175.73k)
================================================================================
谢谢Growlywolf的帮助.
我还有3个问题想继续请教您.
1.模型:我手头分析的是蛋白质序列,因为对模型的选择不是很清楚,于是就根据文献报道,采用了mixed-amino acid model.以下是我运行日志中的Model settings:
Datatype = Protein
Aamodel = Mixture of models with fixed rate matrices
Covarion = No
# States = 20
State frequencies come from the mixture of models
Rates = Gamma
Gamma shape parameter is uniformly dist-
ributed on the interval (0.00,200.00).
Gamma distribution is approximated using 4 categories.
请问在此模型下,还需要具体设置Gamma数值吗?
此外,您建议我用modeltest做一下选取适当距离模型和gamma值,能具体说明一下吗?
2.关于ngen和burnin的设置.
文献报道中的ngen确实都是一两百万,只是我机子的内寸是256M的,分析38个序列运行14小时,也就是产生60000generation.所以我想适当地选择ngen数值.
在运行了上述的命令块后,能得到两个后缀为P的文件(见附件38planta.run1.p和38planta.run2.p),我想请教的是,如果P文件中显示的log likelihood进入收敛状态(即得到稳定值),是否就可以认为命名块中设置的ngen是合适的?我查的资料是,burnin一般选择为ngen的0.1%.不知道我的看法是否是正确的.
3.Clade credibility values:是否大于0.5,就可以认为此节点的可信度可信吗?
再次感谢您的关注.
[url=http://www.dxy.cn/bbs/user/download/9183030/38planta.run1.p][color=#0000ff]38planta.run1.p[/color][/url]
(24.89k)
================================================================================
还有个附件.
[url=http://www.dxy.cn/bbs/user/download/9183035/38planta.run2.p][color=#0000ff]38planta.run2.p[/color][/url]
(24.89k)
================================================================================
[b]markfan wrote:[/b]
1.模型:我手头分析的是蛋白质序列,因为对模型的选择不是很清楚,于是就根据文献报道,采用了mixed-amino acid model.以下是我运行日志中的Model settings:
Datatype = Protein
Aamodel = Mixture of models with fixed rate matrices
Covarion = No
# States = 20
State frequencies come from the mixture of models
Rates = Gamma
Gamma shape parameter is uniformly dist-
ributed on the interval (0.00,200.00).
Gamma distribution is approximated using 4 categories.
请问在此模型下,还需要具体设置Gamma数值吗?
此外,您建议我用modeltest做一下选取适当距离模型和gamma值,能具体说明一下吗?
如果位点替换率不均一，rate选择了gamma，那么就应该设置gamma值。当然你也可以不选这个。
modeltest是一个用来选择核苷酸替代模型的软件，具体你可以google一下。
你用的是氨基酸序列，我看到的一篇文献中是使用tree-puzzle来计算gamma值。
[b]markfan wrote:[/b]
2.关于ngen和burnin的设置.
文献报道中的ngen确实都是一两百万,只是我机子的内寸是256M的,分析38个序列运行14小时,也就是产生60000generation.所以我想适当地选择ngen数值.
在运行了上述的命令块后,能得到两个后缀为P的文件(见附件38planta.run1.p和38planta.run2.p),我想请教的是,如果P文件中显示的log likelihood进入收敛状态(即得到稳定值),是否就可以认为命名块中设置的ngen是合适的?我查的资料是,burnin一般选择为ngen的0.1%.不知道我的看法是否是正确的.
理论上进入收敛状态后是可以的，但是考虑到mcmc的随机性还是比较强，因此我仍然建议ngen不能太少。这个计算速度主要跟CPU有关，你可以考虑借一台CPU强悍一些的电脑来算，^_^。我曾经用AMD Barton 2500+和512M内存算57个序列（序列长度在130 nt左右）也就花了14个小时左右。
burnin的意义是剔除最开始那些处于mcmc低位的随机数据。所谓的0.1%，是对于ngen比较大的情况。burnin后面的数值的是sample数，所以如果ngen是1000000，前10000个不要，而你默认的取样数是每100次取一个杨，那么burnin就是10000/100=100。这个请查阅软件文档。
[b]markfan wrote:[/b]
3.Clade credibility values:是否大于0.5,就可以认为此节点的可信度可信吗?
一般对于MrBayes结果，文献报道是大于90%比较可信，也有取95%的。
================================================================================
非常感谢Growlywolf的热情指导.
看来我还是需要去弄台强悍的机机,不然真是会崩溃的.
^&^smile.gif[/img]
^&^smile.gif[/img]
================================================================================
您有MrBaye:pai.1软件吗?我现在需要这个软件,您能给我吗?!谢谢![email]bingjianla2008@hotmail.com[/email]
================================================================================
呵呵，Growlywolf 解释已很详细，我就不罗嗦了
对于bingjianla2008的需求：MrBayes是免费软件，在其主页可以下载