主题:翻译一个用贝叶斯方法进行系统发育分析的软件的使用说明 -- 空格
构建系统发育树是分子进化研究中的相当重要的一个研究手段。构建树的方法,过去主要有距离法,简约法,最大似然法这三大类。进入新世纪以来,一种叫做贝叶斯分析的方法横空出世。很快得到了广泛的应用。用于系统发育的贝氏分析的软件,也冒出来很多。笔者在学习分子进化的过程中,被若干大牛推荐了一款据称比较可靠的贝氏分析的软件。名叫MrBayes。因为最近翻译翻上了瘾。所以决定把这个软件的说明也翻译出来。希望为想从事这个领域工作的盆友们提供一点便利。
按说,这里应该先介绍一下贝叶斯公式是个什么东西。但是,很遗憾地说,我不太懂这个公式的细节。照大里说,我们知道一个事件的概率是多少,那么就大概能知道,这个事情发生或者不发生的可能性有多大。但是贝大师剑走偏锋,他非要在这个事情已经按一定概率发生了之后,反过来修正我们之前假定的那个概率。个人私下以为,这个过程就是俗称“马后炮”或者“找后帐”的便是。能把这样一个过程发展成一套理论进而形成一个学派。还是顶着若干年正统统计学的口水来实现的。这位贝大师也真可算是人杰。
再说一遍,统计学我不懂,上面关于贝氏分析的理解是我自己的半瓶水,河里如果有懂行的大牛看到我放的厥词,千万务必要技痒难熬义愤填膺地站出来下笔千言指若悬河以正一下视听。免得我的流毒伤害到很多不名真相的群众。
至于我的翻译,好在主要涉及的是软件的使用,关于其理论背景并不很多。希望不会有很多理论上的错误。我尽力去做。具体格式和前面翻译的讲适应进化的统计检验的综述文章一样。
另外,仍然很抱歉地说,这篇文章同样可能不太适合入门级的读者。至少要了解分子生物学,进化生物学的基本知识才会比较好懂。上一次发那篇综述就有朋友提出这一点。不过我实在想不出有什么办法能把自己还要很费力才能看明白的东西写得能更浅显一点。我只能尽力翻译,如果有朋友有问题请提出,我尽量解释。解释不了的就请河里的大小水牛们指点。总以大家都能有所收获为佳。
这个话题或许会很小众,但是我希望对有兴趣的人会是非常有营养的一栋楼。
软件的网址在这里。下载或者查阅都很方便。另外,网上还有一些使用这个软件的教程貌似是繁体字版的。貌似也和原软件的使用手册颇有渊源。仅向这位疑似吴仲义老师学生的兄弟致敬,不管蓝绿。
另外,原手册中有些内容是需要字体和格式才能清楚理解的。在帖子里这些格式通通看不出来了。我只能先把译好的东西贴出来。然后自己整理一个doc格式的版本。回头全文译完后,如果哪位朋友需要,请给我短信留个信箱,我给您发一份全本无码的。
下面是手册的目录部分:
目录 ………………………………………………………………………………………………………………………… 2
1 介绍 ………………………………………………………………………………………………………………… 4
1.1 本手册的一些规范 ……………………………………………………………………………… 4
1.2 获取和安装 ……………………………………………………………………………………………… 4
1.3 开始 …………………………………………………………………………………………………………… 6
1.4 改变MrBayes程序窗口的大小 ……………………………………………………… 6
1.5 获取帮助 ……………………………………………………………………………………………………… 7
1.6 报告和修正程序bug ……………………………………………………………………………… 7
1.7 许可和担保 ……………………………………………………………………………………………… 8
2 教程:一个简单的分析 ……………………………………………………………………………… 8
2.1 快速启动版本 …………………………………………………………………………………………… 9
2.2 向MrBayes程序导入数据 …………………………………………………………………… 9
2.3 指定一个模型 …………………………………………………………………………………………… 11
2.4 设定先验 …………………………………………………………………………………………………… 13
2.5 检查模型 ……………………………………………………………………………………………………… 15
2.6 设置分析参数 …………………………………………………………………………………………… 16
2.7 运行分析 …………………………………………………………………………………………………… 19
2.8 何时终止分析 …………………………………………………………………………………………… 21
2.9 取代模型参数的摘要抽样 …………………………………………………………………… 22
2.10 树和枝长的摘要抽样 ………………………………………………………………………… 24
3.一个分区的数据集的分析 ………………………………………………………………………… 27
3.1.向MrBayes程序导入混合的数据 …………………………………………………… 27
3.2.给数据分区 …………………………………………………………………………………………… 28
3.3.给分区指定模型 …………………………………………………………………………………… 29
3.4.运行分析 ………………………………………………………………………………………………… 30
4.Bayes 3中使用的模型 …………………………………………………………………………… 31
4.1.核苷酸模型 …………………………………………………………………………………………… 31
4.1.1. 简单核苷酸模型 …………………………………………………………………………… 31
4.1.2. 对偶模型 …………………………………………………………………………………………… 33
4.1.3. 密码子模型 ……………………………………………………………………………………… 35
4.2.氨基酸模型 …………………………………………………………………………………………… 36
4.2.1. 固定速率模型 ………………………………………………………………………………… 37
4.2.2. 固定速率模型的评估 …………………………………………………………………… 37
4.2.3. 速率可变模型 ……………………………………………………………………………… 37
4.3.限制位点(二进制)模型 ………………………………………………………………… 39
4.4.标准离散(表型)模型 ……………………………………………………………………… 40
4.5.简约模型 …………………………………………………………………………………………………… 42
4.6.位点间速率差异模型 …………………………………………………………………………… 43
4.6.1. gamma分布的速率模型 ……………………………………………………………… 43
4.6.2. 自动互关联的{autocorrelated}gamma模型 ……………… 43
4.6.3. 不变位点的比例 …………………………………………………………………………… 44
4.6.4. 分区速率模型(位点特异的) ……………………………………………… 45
4.6.5. 推断位点速率 ……………………………………………………………………………… 45
4.7.树间速率差异:Covarion 模型 ………………………………………………… 46
4.8.拓扑和枝长模型 ………………………………………………………………………………… 47
4.8.1. 非约束和约束的拓扑 ………………………………………………………………… 47
4.8.2. 非时钟树(标准树) ……………………………………………………………… 48
4.8.3. 严格时钟树 ………………………………………………………………………………… 48
4.8.4. 宽松时钟树 ………………………………………………………………………………… 49
4.9.分区的模型 ………………………………………………………………………………………… 49
4.10.祖先状态重建 ………………………………………………………………………………… 50
5.常见问题 ………………………………………………………………………………………………… 51
6.程序第二版和第三版的区别 …………………………………………………………… 57
7.高级话题 ………………………………………………………………………………………………… 59
7.1.编译MrBayes …………………………………………………………………………………… 59
7.1.1. 用GNU的make来编译 ……………………………………………………………… 59
7.1.2. 用Code Warrior 或 Visual Studio编译 ……………… 61
7.2.编译和运行并行版的MrBayes …………………………………………………… 61
7.2.1.苹果机上的并行版 ……………………………………………………………………… 61
7.2.2. linux集群上的MPI版 ……………………………………………………………… 62
7.3.源代码 …………………………………………………………………………………………………… 62
8.致谢 ……………………………………………………………………………………………………………… 63
9.参考文献 …………………………………………………………………………………………………… 63
附录:MrBayes中用到的进化模型和建议{proposals}
内容已整理,此楼已合并。
MrBayes 3是一个将Bayes推断应用于系统发育{phylogeny}分析的程序。本程序使用命令行界面,可以在多种系统平台上运行,包括苹果和linux的集群{cluster}。注意程序使用的计算机应该足够快而且有相当数量的内存(内存的需要量依赖于数据矩阵的大小,程序可能最多需要到几百M的内存)程序是按照最快运算速度优化的而不是按照最小内存量来优化的。
本手册告诉您如何使用程序。在这一部分只是介绍一下本程序的概貌,之后我们会领您看一些简单的分析(本手册的第二部分)。之后在第三部分我们会进入一个复杂些的分析,向您展示程序的更多的性能{capabilities}。第四部分我们向您介绍一下程序所使用的模型。之后的第五部分会回答一些常见问题。第六部分讨论一些第二版和第三版的差异。最后的第七部分,我们会给出一些更细节的建议,主要是关于如何编译和运行并行版程序的问题。第七部分同样包含一些开发者最关心的问题,包括如何对程序代码作微调以及如何加入MrBayes项目并作出贡献。最后有一个附录部分,给出了一些图表,图示化地展示了程序中用到的模型和建议的机制{proposal mechanisms}的概貌。至于MrBayes程序更细节的命令和选项的信息,请从网站上下载命令参考或查看命令行的帮助(见第1.4节获取帮助)。
所有的命令行帮助的信息都可以在使用在线程序时获得。
本手册假定读者熟悉Bayes系统发育分析的基本概念。如果您是这个领域的新手,我们推荐您看看最近的几篇综述,包括 Holder 和 Lewis (2003), Lewis (2001) 和 Huelsenbeck 等人 (2001, 2002)的。同样有价值的是一些早期的介绍Bayes系统发育方法的一些文章(Li 1996; Mau, 1996; Rannala and Yang, 1996; Mau and Newton, 1997; Rannala and Yang, 1997; Larget and Simon, 1999; Mau, Newton and Larget, 1999; Newton, Mau and Larget, 1999)。基本的MCMC技术在Metropolis et al. (1953) 和Hastings (1970)中有提到。而MrBayes 程序中用到的 Metropolis-coupled MCMC 方法在Geyer (1991)的文章中有提到。
1.1 本手册的一些规范
读者在屏幕上和输入文件中看到的内容用斜体字体表示。需要用户输入的字体用加粗的斜体字体表示。<手册原文使用的是plain typewriter字体和加粗的plain typewriter字体来表示>
1.2 获取和安装
MrBayes程序可以在MrBayes网站上免费得到(http://mrbayes.net)。如果有人给来您一个MrBayes 3程序,我们强烈在建议您从这个网站上下载最新的版本。该网站上同样给出了关于MrBayes用户邮件组和如何参加开发项目的信息。
MrBayes是一个简单的使用命令行界面的程序,因而在不同的操作平台上会有相同的表现,无论苹果{Macintosh}系统,瘟到死{windows}系统,Unix系统。各平台有彼此独立的下载包。苹果的包有两个版本,名为MrBayes3.1的标准版本(program icon one copy of Reverend Bayes’s portrait)和名为MrBayes3.1p的用于苹果机集群上的并行计算版本(program icon four portraits of Reverend Bayes)。关于MrBayes苹果机并行版的更多信息请参考POOCH的安装信息,具体可参考本手册后面第七部分。而瘟到死版只包含串行程序,解压缩后可以直接运行。
如果您需要在unix或犁牛渴死{linux}上运行MrBayes程序,您需要从源代码编译。在这种情况下,您需要用下面的命令解压缩mrbayes3.1_src.tar.gz文件即可:
gunzip MrBayes_3.1_src.tar.gz
然后
tar –xf MrBayes_3.1_src.tar
注意gunzip命令解开压缩包。而 tar -xf 命令从解压后的.tar 文件中释放打包的内容(注意 .gz 后缀在解压缩后就被去掉了)。然后一步的工作是编译程序。解压后的包中包含一个名为Makefile的文件,可以生成串行版的编译程序。只要输入make就可以编译这些内容。典型的编译时的输入类似下面的样子:
ronquistg5:mrbayes>ls
mrbayes3.1_src.tar.gz
ronquistg5:mrbayes>gunzip MrBayes_3.1_src.tar.gz
ronquistg5:mrbayes>ls
mrbayes3.1_src.tar
ronquistg5:mrbayes>tar -xf MrBayes_3.1_src.tar
ronquistg5:mrbayes>make
gcc -DUNIX_VERSION -O3 -Wall -Wno-uninitialized -c -o mb.o mb.c
gcc -DUNIX_VERSION -O3 -Wall -Wno-uninitialized -c -o mcmc.o mcmc.c
gcc -DUNIX_VERSION -O3 -Wall -Wno-uninitialized -c -o bayes.o bayes.c
gcc -DUNIX_VERSION -O3 -Wall -Wno-uninitialized -c -o command.o command.c
gcc -DUNIX_VERSION -O3 -Wall -Wno-uninitialized -c -o mbmath.o mbmath.c
gcc -DUNIX_VERSION -O3 -Wall -Wno-uninitialized -c -o model.o model.c
gcc -DUNIX_VERSION -O3 -Wall -Wno-uninitialized -c -o plot.o plot.c
gcc -DUNIX_VERSION -O3 -Wall -Wno-uninitialized -c -o sump.o sump.c
gcc -DUNIX_VERSION -O3 -Wall -Wno-uninitialized -c -o sumt.o sumt.c
gcc -DUNIX_VERSION -O3 -Wall -Wno-uninitialized -lm mb.o bayes.o command.o
mbmath.o mcmc.o model.o plot.o sump.o sumt.o -o mb
ronquistg5:mrbayes>
编译程序一般会在mcmc.c文件编译时停留若干分钟。这是正常现象。这是最大的一个源代码文件,优化它们需要比较长一点的时间。
我们缺省使用的C编译器是gcc。在绝大多数系统上都缺省安装了这个编译器。如果您的机器上没有这个编译器,或者您需要生成MPI版本,您需要修改编译信息中的内容。具体参见本手册第七部分。可执行的串行程序名为”mb”。若要运行它。只需在编译的子目录里键入
./mb
即可。其中“./”前缀是告知unix系统你要运行的是本地本子目录中的该程序的一个拷贝。如果你需要经常运行mb程序,你可能需要把这个程序加到你的系统路径{path}中去。请查阅您的unix手册或咨询您身边的unix系统专家来解决此问题。
所有三个MrBayes程序都附带样例数据文件。他们被设计来展示程序所能运行的各种分析。您可以以这些例子作为模板来设计您自己的分析。其中有两个文件primates.nex和cynmix.nex在后文第二和第三部分中被使用。
1.3 开始
双击MrBayes应用图标就可以开始程序,或者键入./mb后回车。您将能看到如下的信息。
MrBayes v3.1
(Bayesian Analysis of Phylogeny)
by
Fredrik Ronquist and John P. Huelsenbeck
School of Computational Science
Florida State University
Section of Ecology, Behavior and Evolution
Division of Biological Sciences
University of California, San Diego
Distributed under the GNU General Public License
Type "help" or "help <command>" for information
on the commands that are available.
MrBayes >
注意作者的次序是每次程序启动时随机变化的。所以如果发现这一次启动程序后的作者次序和上一次启动程序不一样时请保持情绪稳定,不必大惊小怪。注意底部的MrBayes >提示符。出现这个提示符说明程序已经作好工作准备。
1.4 改变MrBayes程序窗口的大小
有的MrBayes的命令会生成很多行的输出结果。您可能会需要把窗口尺寸扩大以方便地读您的结果。在苹果和Unix机器上,您只要拖动窗口的边缘就可以了。但是在瘟到死系统上(包括2000,XP,NT血色),您就不能这么简单地操作。您需要用右键单击MrBayes窗口顶部的蓝色标题栏,然后在出现的菜单中选择“属性”。在跳出来的新窗口中,确认“布局”{Layout}项被选定,然后将窗口尺寸和窗口缓冲大小按您的需要修改。
1.5 获取帮助
在mrbayes>提示符下键入help 可以看到一个MrBayes中可用的命令的列表。绝大多数命令允许你为不同的参数设定值(选项)。如果你键入 help <command> ,其中的尖括号表示某个特定的命令。那你就可以得到这个特定命令的帮助,包括了所有可选的参数的介绍。对于大多数命令,在介绍的结尾能看到所有参数目前的设定值。例如,试试执行help lset 或 help mcmc命令。则lset的设置会被显示如下:
Default model settings:
Parameter Options Current Setting
----------------------------------------------------------------------------------------------
Nucmodel 4by4/Doublet/Codon 4by4
Nst 1/2/6 1
Code Universal/Vertmt/Mycoplasma/
Yeast/Ciliates/Metmt Universal
Ploidy Haploid/Diploid Diploid
Rates Equal/Gamma/Propinv/Invgamma/Adgamma Equal
Ngammacat <number> 4
Nbetacat <number> 5
Omegavar Equal/Ny98/M3 Equal
Covarion No/Yes No
Coding All/Variable/Noabsencesites/
Nopresencesites All
Parsmodel No/Yes No
------------------------------------------------------------------------------------------------
注意MrBayes3中支持命令和选项的缩写,所以在很多情况下可以键入命令开头的少数几个字母来代替命令的全称。
在命令行界面下,有一个全部命令和选项的列表。这个列表可以在程序的网站上下载得到(http://www.mrbayes.net)。你同样可以得到一个ASCII码格式的版本方式为在命令行提示符后输入manual命令。更多的帮助有Jeff Bates 在mrbayes网站上做好的一系列html格式的超链接。最后,你可以通过mrbayes的用户邮件组与其他mrbayes程序的使用者保持联系(提交信息到www.mrbayes.net )。
1.6 报告和修正程序bug
如果您在使用中发现了bug。您可以在SourceForge上的报告错误功能来通知我们。我们会对您的行为深表感谢。提交错误报告的功能在网站上也有介绍(www.mrbayes.net)当您提交一个错误报告时,请确认您上传了一个包含数据集和产生错误的命令列表的数据文件。
如果问题出现在运行MCMC分析的过程中(在执行”mcmc”命令之后),请确认您的错误报告中包含了mcmc命令行中的固定种子{fixed seed}和交换种子{swapseed},最好还包括一个小的数据集。这对我们分析您报告的错误会非常有帮助。如果这个问题被确定是MrBayes源代码的问题。SourceForge的追踪器程序会给您发一封email通知您。然而,需要注意的是,确认修正了源码错误的新的可执行程序会在晚一点的时间内发布。
高级用户可能会倾向于自行确定源码中的错误。请到本手册第七部分查询如何参与源码维护,算法改进或功能扩展等工作。
1.7 许可和担保
MrBayes是一个自由软件。您可以遵循自由软件基金会的GNU通用公共授权{General Public License}来修改或重新发布它。无论第二版或以后的版本都是如此。
本程序希望是有用处的,但是并不保证这一点!甚至没有暗示此软件可以对某类问题具有特别的适用性!请参考GNU通用公共授权来获得更多的细节(http://www.gnu.org/copyleft/gpl.html)。
抱歉拖了很久才回来更新。河里有墙只是小部分原因。一个主要原因是这个工作量比我想像的大,也难。只能说尽力去做。因为这部分太长了。所以我会把开头第一段放在这里,后面各小节独立出帖。(也免得没有新帖出先这个楼被沉下去,呵呵。。。)
---废话分隔线---
这一部分的内容基于primates.nex数据文件。它将引导读者进行一个基本的系统发育Bayes的MCMC分析,并解释程序最重要的若干特性。为了照顾部分性急的读者,本部分的教程将分为两个部分,第一部分是快速型的,可以立即开始一个分析并得到结果。第二部分是详细型的,里面会介绍这个分析的过程中的更多细节。
占楼
占楼
崭露
-
希望不久的将来我可以把这里也填满。那样我就可以去翻译PAML了。呵呵。。。
怎么都是占楼啊?
期待,虽然估计自己看不懂。
典型的使用MrBayes程序进行系统发育分析的Bayes分析有四个步骤:
首先,读入nexus数据文件
其次,设定进化模型
然后,进行分析
最后,总结结果
下面会依次介绍每个步骤在执行中的细节内容。
1.在MrBayes> 提示符后输入execute primates.nex 命令。这将把数据导入程序中。数据文件(primates.nex)必须和程序在同一个子目录中。如果遇到了什么问题,请看2.2节有没有什么可能的解决方案。如果你运行的是自己的数据文件,请注意里面不要有可能改变程序行为的MrBayes的命令,删除这些命令,或者用中扩弧把它们注释掉。这样你就可以得到和本教程一样的输出。
2.在MrBayes> 提示符后输入lset nst=6 rates=invgamma 命令。此命令将设定进化模型为GTR模型,同时位点间速率的差异将遵循gamma分布,并有一定比例的不变位点。如果你的序列是蛋白质序列,或者或者你想使用非缺省的设定,请参考本手册其他部分,特别是3到5部分及附录。
3.在MrBayes> 提示符后输入mcmc ngen=10000 samplefreq=10 命令。这一命令将使你得到后验概率分布中的1000次抽样。如果数据集更大,则可能需要更长的时间或降低抽样的频率(目前的缺省频率是每百次抽一个样)。你可以在屏幕输出的最后一栏中看到一个预测的程序执行时间。
如果切断频率的标准变异{standard deviation of split frequencies}在100000代之后其数值小于 0.01,对程序询问是否继续的问题”Continue the analysis?(yes/no)”回复no。这样就可以终止程序。或者,继续追加世代数,直到频率值小于0.01。
4.在MrBayes> 提示符后输入sump burnin=250命令以显示分析使用的参数值(或者使用一个相当于你的样本四分之一大小的数值)。程序将给出一个样本的取代模型参数的摘要的列表。包括了每个参数的均值,模式,95%置信区间。请确认所有参数的潜在的标尺缩减参数(PSRF){potential scale reduction factor}非常接近1.0。如果不是,就需要运行更长时间的分析。
输入sumt burnin=250(或者其他随便什么接近你样本大小25%达到数值)以分析树形。程序会给出一个包含每枝的枝长和后验概率值的进化树。此树会同样保存到文件,以用其他树形软件打开并查看,例如TreeView,MacClade或Mesquite等等。
分析过程就是这么回事,不会再复杂了。但是,熟悉了这个套路之后,读者八成会对程序执行时后台发生的事情有兴趣。本节的其他部分为您解释了每个步骤的细节内容,并介绍了在分析过程中所有的隐含假定以及MrBayes程序的设置等等你可能想知道的东西。
向mrBayes丞相导入数据需要Nexus格式的文件。这种格式可以使用四种类型的数据,核苷酸序列,氨基酸序列,表型数据(标准数据)限制性位点数据。或者可以是这四种数据的混合。在本教程中使用的Nexus文件(primates.nex)包含12种灵长类的线粒体DNA序列。
Nexus文件是一种简单的ASCII文本格式的文件。它以#NEXUS作为首行的开头。其余部分被分割为几个不同的块。primates.nex文件看起来是这个样子的:
#NEXUS
begin data;
dimensions ntax=12 nchar=898;
format datatype=dna interleave=no gap=-;
matrix
Saimiri_sciureus AAGCTTCATAGGAGC ... ACTATCCCTAAGCTT
Tarsius_syrichta AAGCTTCACCGGCGC ... ATTATGCCTAAGCTT
Lemur_catta AAGCTTCACCGGCGC ... ACTATCTATTAGCTT
Macaca_fuscata AAGCTTCACCGGCGC ... CCTAACGCTAAGCTT
M_mulatta AAGCTTCACCGGCGC ... CCTAACACTAAGCTT
M_fascicularis AAGCTTTACAGGTGC ... CCTAACACTAAGCTT
M_sylvanus AAGCTTTTCCGGCGC ... CCTAACATTAAGCTT
Homo_sapiens AAGCTTTTCTGGCGC ... GCTCTCCCTAAGCTT
Pan AAGCTTCTCCGGCGC ... GCTCTCCCTAAGCTT
Gorilla AAGCTTCTCCGGTGC ... ACTCTCCCTAAGCTT
Pongo AAGCTTCACCGGCGC ... ACTCTCACTAAGCTT
Hylobates AAGTTTCATTGGAGC ... ACTCTCCCTAAGCTT
;
end;
文件中只有一个数据块。数据块用begin data; 标识块起始,跟着是维度标识{dimensions statement},格式标识{format statement}和矩阵标识{matrix statement}。dimensions statement 必须包含指定分类单元数目的ntax参数,以及指定每条联配好的序列的联配字符数目的nchar参数。format statement 必须指定数据类型,例如 datatype=DNA(或者RNA或Protein或standard或Restriction)。format statement 中还必须有以下内容
gap=-(或任意字符以表示联配中的空位)
missing=?(或任意其他字符以表示文件中的丢失数据)
interleave=yes,如果数据矩阵是被分成多行的
以及match=.(或者任意字符以表示联配中匹配到的字符)。
跟着format statement的是matrix statement,通常以独立行中的matrix这个单词开始,然后跟着的是联配的序列。每行以序列名称开头,以至少一个空格将名称与序列分开。数据块以end; 标识结束。注意这四个statments都要以分号结尾。这个分号必须有而不能省略<习惯用C或JAVA的朋友们有福了>。另外要注意的是,虽然文件中联配的序列会有许多行,矩阵的描述仍然是一个有规律的标识。所以它也是要以分号结束的。我们的例子文件包含的序列是没有分成多行的{non-interleaved}块格式。不分行的情况是缺省的,不过你可以在格式标识中增加interleave=no来使之允许多行。
Nexus 数据文件可以从诸如MacClade或 Mesquite这样的程序生成。但是,需要注意的是,MrBayes程序并不支持全部的Nexus 标准,所以读者可能必须先编辑一下您的数据文件以适应MrBayes程序的要求。尤其要说的是,MrBayes程序对每种数据类型使用的是固定的标志{symbols},而不支持对数据格式使用自定义的标志,即不支持在 “format”命令中使用“symbols” 和 “equate”选项。所支持的标志对DNA数据是{A, C, G, T, R, Y, M, K, S, W, H, B, V, D, N},对RNA数据是{A, C, G, U, R, Y, M, K, S, W, H, B, V, D, N},对蛋白质数据是{A, R, N, D, C, Q, E, G, H, I, L, K, M, F, P, S, T, W, Y, V, X},对限制性数据是{0,1}二进制的,对标准表型数据是{0, 1, 2, 3, 4, 5, 6, 5, 7, 8, 9}。另外,对DNA和RNA的标准单字母模糊表示之外,不确定的核苷酸还可以用Nexus parenthesis 或 curly braces notation来表示。例如,一个分类单元多态性的状态2和3可以被编码为(23),(2,3),{23},或 {2,3};而如果是氨基酸A或F则可以表示为(AF),(A,F),{AF}或{A,F}。和其他大多数系统发育的程序类似,MrBayes能有效地处理多态性和不确定数据。所以无论使用圆括号或花括号都可以。如果你的数据矩阵中有其他MrBayes不支持的标志,你需要在执行MrBayes之前替换掉这些符号。你同样需要删除格式行中的Equate和Symbols标识。与标准Nexus格式不同,MrBayes支持在数据块中爆可混合数据类型,请参看本手册第三部分以了解详细内容。
向MrBayes程序导入数据需要在Mrbayes>提示符后输入execute <filename>命令,其中<filename>是输入文件名。例如,为了执行我们的例子文件。要输入的命令是:
execute primates.nex,或者比较省事的方法是输入 exe primates.nex命令。注意输入文件名必须和MrBayes程序在同一个子目录下(或者你必须给出到文件名的路径),并且输入文件名不能有空格。如果每一个步骤都正常进行,MrBayes将会回馈说他已经读入了Nexus文件里DATA block中的数据。其屏幕显示如下:
Executing file "primates.nex"
DOS line termination
Longest line length = 915
Parsing file
Expecting NEXUS formatted file
Reading data block
Allocated matrix
Matrix has 12 taxa and 898 characters
Data is Dna
Data matrix is not interleaved
Gaps coded as -
Setting default partition (does not divide up characters).
Taxon 1 -> Tarsius_syrichta
Taxon 2 -> Lemur_catta
Taxon 3 -> Homo_sapiens
Taxon 4 -> Pan
Taxon 5 -> Gorilla
Taxon 6 -> Pongo
Taxon 7 -> Hylobates
Taxon 8 -> Macaca_fuscata
Taxon 9 -> M_mulatta
Taxon 10 -> M_fascicularis
Taxon 11 -> M_sylvanus
Taxon 12 -> Saimiri_sciureus
Setting output file names to "primates.nex.<run<i>.p/run<i>.t>"
Successfully read matrix
Exiting data block
Reached end of file
所有的命令都在MrBayes>提示符后面输入。有两个命令是必须的,lset和prset。因为这哥俩的作用是指定分析使用的是哪个进化模型。showmodel命令可以用来检查模型设定是否正确,通常会使用这个命令以策万全。另,lset命令指定模型的结构,prset命令按给定模型的参数定义先验概率的分布。接下来,我们会为线粒体序列的进化指定GTR+I+_模型(即通用时间可逆{General Time Reversible}模型,带可变位点比例,位点间速率呈gamma分布)并检查所有相关的先验。若您对分子进化的离散模型不甚了了,我们建议您看一些相关的介绍性文章,例如Felsenstein的(2004)这篇文章。
一般说来,使用help lset 命令来入门是一个好办法。忽略开头的帮助信息而关注于输出末尾的表格。该表格显示了当前的设置,一般说来会是下面这个样子的:
Model settings for partition 1:
Parameter Options Current Setting
---------------------------------------------------
Nucmodel 4by4/Doublet/Codon 4by4
Nst 1/2/6 1
Code Universal/Vertmt/Mycoplasma/
Yeast/Ciliates/Metmt Universal
Ploidy Haploid/Diploid Diploid
Rates Equal/Gamma/Propinv/
Invgamma/Adgamma Equal
Ngammacat <number> 4
Nbetacat <number> 5
Omegavar Equal/Ny98/M3 Equal
Covarion No/Yes No
Coding All/Variable/Noabsencesites/
Nopresencesites All
Parsmodel No/Yes No
------------------------------------------------
首先注意这个表的题目是Model settings for partition 1。在缺省的情况下,MrBayes为你数据块中的每类数据建立一个分区{partition}。如果你的数据只有一种类型。那么就会之有一个分区。改变分区的方法见本手册第三部分。
名为Nucmodel的设置项允许你设置DNA模型的各类参数。其中,Doublet选项针对的是核糖体DNA匹配成茎的区域;Codon选项是按照其密码子来分析DNA序列。在例子中,我们将用标准的核苷酸取代模型,其中缺省的 4by4 选项就已经是最合适的了。所以我们会保留Nucmodel的缺省设置。
取代模型的通用设置是由Nst选项的设置来确定的。在缺省的情况下,所有的取代具有相同的速率(Nst=1)。这对应的是F81模型。(或者,可以相当于设定prset参数强制稳定状态频率彼此相等的JC模型,详见下文)。我们所需要的是GTR模型(Nst=6)而不是F81模型。所以我们键入并执行lset nst=6 命令。MrBayes会声明已经修改了模型的设定。
Code设置只有当Nucmodel设定为codon时才有效。而Ploidy在目前的例子里同样用不到。需要修改的是Rates的设置,从缺省的Equal(意即取代在不同位点间无差别)改为Invgamma(位点间速率差异服从gamma分布同时保留一部分非变异位点)。实现此改变所要执行的命令是lset rates=invgamma 。同样地,MrBayes会声明模型的参数已经被重新设定过了。
上面两步可以合并为如下一行的执行方式:lset nst=6 rates=invgamma
我们保留了Ngammacat的设置为缺省的4。通常情况下,4个速率的域{categories}已经足够了。增加域的数目有可能会增加似然计算的精确度。然而,计算所需要的时间随设定的域数目的增加而线性增加,而且多数情况下对结果的影响是负的。<这段话前后逻辑有点不太明白>
其他的设置里,只有Covarion 和 Parsmodel 两个设置是与单取代模型有关系的。不过,因为我们的数据既不需要俭约{parsimony}模型也不需要协方差{covariotide}模型。所以这两个的设置就不用改了。如果你键入help lset 命令来检查目前的设置是否正确,所看到的应该是下面的表的样子:
Model settings for partition 1:
Parameter Options Current Setting
------------------------------------------------------------------
Nucmodel 4by4/Doublet Codon 4by4
Nst 1/2/6 6
Code Universal/Vertmt/Mycoplasma/
Yeast/Ciliates/Metmt Universal
Ploidy Haploid/Diploid Diploid
Rates Equal/Gamma/Propinv/Invgamma/Adgamma Invgamma
Ngammacat <number> 4
Nbetacat <number> 5
Omegavar Equal/Ny98/M3 Equal
Covarion No/Yes No
Coding All/Variable/Noabsencesites/
Nopresencesites All
Parsmodel No/Yes No
------------------------------------------------------------------
现在我们要为我们的模型设定先验{priors}。模型中有6种类型的参数:拓扑,支长,4种核苷酸的稳定频率,6种不同的核苷酸取代的速率,非稳定位点的比例,以及速率差异所服从的gamma分布的形状参数。在MrBayes中,缺省的先验一般都能很好地工作。所以如果你没耐性的话可以不用修改这些参数。所以忽略这一节直接进入下一节“检查模型”也是可以的。然而,学习这写先验并理解这些缺省的设置及可能的其他选项是一个很好的体验。键入help prset 可以得到你的模型的缺省设置的列表。帮助信息最后的表格是这样的:
Model settings for partition 1:
Parameter Options Current Setting
------------------------------------------------------------------
Tratiopr
Beta/Fixed
Beta(1.0,1.0)
Revmatpr
Dirichlet/Fixed
Dirichlet(1.0,1.0,1.0,1.0,1.0,1.0)
Aamodelpr
Fixed/Mixed
Fixed(Poisson)
Aarevmatpr
Dirichlet/Fixed
Dirichlet(1.0,1.0,...)
Omegapr
Dirichlet/Fixed
Dirichlet(1.0,1.0)
Ny98omega1pr
Beta/Fixed
Beta(1.0,1.0)
Ny98omega3pr
Uniform/Exponential/Fixed
Exponential(1.0)
M3omegapr
Exponential/Fixed
Exponential
Codoncatfreqs
Dirichlet/Fixed
Dirichlet(1.0,1.0,1.0)
Statefreqpr
Dirichlet/Fixed
Dirichlet(1.0,1.0,1.0,1.0)
Ratepr
Fixed/Variable=Dirichlet
Fixed
Shapepr
Uniform/Exponential/Fixed
Uniform(0.0,50.0)
Ratecorrpr
Uniform/Fixed
Uniform(-1.0,1.0)
Pinvarpr
Uniform/Fixed
Uniform(0.0,1.0)
Covswitchpr
Uniform/Exponential/Fixed
Uniform(0.0,100.0)
Symmetricbetapr
Uniform/Exponential/Fixed
Fixed(Infinity)
Topologypr
Uniform/Constraints
Uniform
Brlenspr
Unconstrained/Clock
Unconstrained:Exp(10.0)
Speciationpr
Uniform/Exponential/Fixed
Uniform(0.0,10.0)
Extinctionpr
Uniform/Exponential/Fixed
Uniform(0.0,10.0)
Sampleprob
<number>
1
Thetapr
Uniform/Exponential/Fixed
Uniform(0.0,10.0)
Growthpr
Uniform/Exponential/
Fixed/Normal
Fixed(0.0)
------------------------------------------------------------------
我们需要关注的是如下选项:
Revmatpr 设定GTR速率矩阵中的6种取代速率
Statefreqpr 设定GTR速率矩阵中的稳定核苷酸频率
Shapepr 设定 速率变异的gamma分布的形状参数
Pinvarpr 设定非编译位点的比例
Topologypr 设定树拓扑
Brlenspr 设定支长
对于Revmatpr和 Statefreqpr 模型来说,缺省的先验概率密度是一个扁平的狄利克雷{a flat Dirichlet},即所有的值都是1.0。当我们需要从不知道任何先验知识的数据中估计这些参数时。这个设定还算是合适的。虽然固定速率和核苷酸频率是可能的,但是通常情况下并不推荐做这样的假设。不过,有的情况下,固定核苷酸频率相等也是有必要的,例如JC或SYM模型。实现这样的修改的命令如下:
prset statefreqpr=fixed(equal)
如果我们想设定一个比缺省的扁平狄利克雷更强调相等的核苷酸频率的先验,可以用这个命令:
prset statefreqpr = Dirichlet(10,10,10,10)
或者如果要更强调相等的核苷酸频率,可以用这个命令
prset statefreqpr=Dirichlet(100,100,100,100)
狄利克雷分布的数字的总数决定了分布在多大程度上聚集起来,而数字间的平衡决定了每种核苷酸的期望比例(按ACGT的顺序)。通常,在狄立克雷分布的参数和观察值之间会有一个连接{connection}。你可以考虑一个(150,100,90,140) 这样的分布,如果在某组序列中有150个A,100个C,90个G以及140个T,而这组序列与你所要分析的序列组是彼此独立的,但是又与你所分析的序列确定地有某种联系。那么就可以用这样的分布作为先验。
在我们的分析中,我们会比较谨慎地处理这个问题。我们保留了缺省的设定。如果你想按上面说的方法修改这个先验。可以使用下面的命令:
prset statefreqpr=Dirichlet(1,1,1,1)
或者用这样比较省事的写法:
prsst = Dir(1,1,1,1)
类似地,我们也可以设定取代速率的先验为缺省的狄利克雷分布值(1,1,1,1,1,1)。
ShapeprI 参数设定了速率变异所服从的gamma分布的alpha形状参数的先验值。我们同样会保留这个值的缺省值,使之为一个alpha值跨度很大的统一的分布。非变异位点的比例的先验值由Pinvarpr 参数设定。其缺省值为0到1之间的统一分布。因为我们没有假定任何关于非变异位点比例的先验的知识,所以这个设定是合适的。
对于拓扑来说,缺省的Uniform设置Topologypr 参数为在所有的区间{distinct}上有相等的概率,充分解出拓扑{fully resolved topologies}。<这两句话貌似是拓扑学中的专业说法,不是很懂>。可选的是约束树上的某些节点使之总是在场{The alternative is to constrain some nodes in the tree to always be present}。但是我们在这里不会用这个方法。
Brlenspr 参数可以设成时钟约束的{clock-constrained}或者非约束的。对于没有分子钟的树,其支长的先验可以是指数的或者uniform的。缺省的参数值为10的指数先验适用于绝大多数分析。其期望值为1/10 = 0.1但是允许更广泛的支长的数值(从0到无穷)。因为短支的似然值比长支的似然值变异更快,所以支长的指数先验比uniform先验更为接近信息统一{an exponential prior on branch lengths is closer to being uninformative than a uniform prior}。<这里的uninformative很难有好的翻译>
一直很想用之
但始终静不下来学
有空多请教你