由 Fredrik Ronquist, John Huelsenbeck & Maxim Teslenko 开发的 MrBayes 是一款相当稳定好用的系统树重建软件。其基于贝叶斯推断和模型选择对系统树进行重建,命令行简单易学,教程详细,帮助文件调用方便,关键是建树结果比较稳定。在Windows、Linux和Mac系统都可运行,具体安装方法参见 INSTALL。
准备输入文件
输入文件格式一般是扩展名为 .nex
的 Nexus 文件。该格式除了包含基因序列文件之外,还可以包含一系列命令行,其格式具体如以下示例。该示例文件仅包含2个个体、4个基因。在实际运算中,可根据这一格式将序列和命令行逐次定义好(命令行也可另在MrBayes中定义)。
#NEXUS
Begin data;
Dimensions ntax=2 nchar=225;
Format datatype=DNA interleave=yes gap=- missing=?;
Matrix
Angulyagra_sp._Nam_Dinh_Angu.sp.1 CCGCTGAATTTAAGCATATCACTAAG
Bellamya_aeruginosa_Baiyang_Lake_BY1 CCGCTGAATTTAAGCATATCACTAAG
Angulyagra_sp._Nam_Dinh_Angu.sp.1 TATGGCTCGTACCAAGCAGACNGCTCGTAAATCAACCG-GAGGAAAAGC
Bellamya_aeruginosa_Baiyang_Lake_BY1 TATGGCTCGTACCAAGCAGACGGCTCGCAAATCGACCG-GAGGAAAAGC
Angulyagra_sp._Nam_Dinh_Angu.sp.1 --------------------------------------------------TAACCTGCCCGGTGA
Bellamya_aeruginosa_Baiyang_Lake_BY1 CGCCTGTTTATCAAAAACATGTCTTTTTGATTTTAGAAATATAAAAAGTCTGACCTGCCCGGTGA
Angulyagra_sp._Nam_Dinh_Angu.sp.1 ----------------------------------------------------------------GCATCCTGAGGTTTATATTTT
Bellamya_aeruginosa_Baiyang_Lake_BY1 GGTCAACAAATCATAAAGATATTGGAACTTTATATATTTTGTTTGGTGTATGGTCAGGATTAGTTGGTACTGGACTTAGTATTTT
;
End;
begin mrbayes;
[partition]
charset 28S = 1-26;
charset H3 = 27-75;
charset 16S = 76-140;
charset COI = 141-225;
partition favored = 4: 28S, H3, 16S, COI;
set partition=favored;
lset applyto=(1) rates=gamma nst=6;
lset applyto=(2,3,4) rates=invgamma nst=6;
unlink revmat=(all) pinvar=(all) shape=(all) statefreq=(all);
prset applyto=(all) ratepr=variable;
end;
对上述示例的格式说明之,即可了解 MrBayes 的初步用法。标准的 Nexus 格式的详细说明可参见 NEXUS FORMAT,而 MrBayes 有些许变动,这里仅探讨适用于 MrBayes 的 Nexus 文件格式。
-
适用于 MrBayes 的 Nexus 文件可包含核酸序列、氨基酸序列、形态学数据(也称标准数据 standard)或二项数据(也称限制数据 restriction),以及任意这四种类型数据的组合,这些数据则是由特定字符组成,如核酸序列中的核苷酸字母或蛋白质序列中的氨基酸代码字母。具体地说,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, 2, 3, 4, 5, 6, 7, 8, 9;二项数据可用字符为 0, 1。另外,一般情况下,用
-
代表序列中的 gap,用?
代表序列中的疑问位点,用N
代表序列中的缺失位点。 -
文件的前5行即对序列数据进行了定义。第一行 #NEXUS 为文件备注,第二行数据开始,固定不变。第三行开始定义:
ntax
为样本数,nchar
为样本的序列总长度。第四行定义序列类型:通常使用较多的是DNA序列数据和形态学数据的混合数据,这时datatype=mixed(DNA:1-200, standard:201-300)
,括号中定义的是不同类型序列的具体长度。如果全部是DNA序列,则datatype=DNA
,如果全部是形态学数据,则datatype=standard
。接着,这里每个样本的序列数据有4个不同的基因,并且根据基因顺序依次列出来的,因此interleave=yes
即表示每个样本的序列数据是可被分割的。然后定义gap=-
,缺失数据missing=?
。第五行固定不变。 -
按照样本个体和基因顺序,依次列出序列数据。序列数据必须先行比对,对齐后的序列才可用于 MrBayes,序列对齐可使用 MEGA 中的一系列方法。要注意,每个个体的序列要与该个体名称在同一行,用若干空格(或 Tab)隔开,但必须确保所有样本的序列起始位置是对齐的。
-
序列全部列出后,换行,添加
;
-
最后再换行,添加
end;
定义基因分区
由于示例数据是由4个基因数据组成的序列集,因此首先要定义各个基因区域。默认情况下,MrBayes 根据数据类型对序列分区,这里的示例数据则需要根据各基因的长度,手动分区。
-
第一行
begin mrbayes;
固定不变。第二行开始定义分区名称和各区长度,使用的命令是charset
,按照基因名称定义各区长度,如28S=1-26;
和H3=27-75;
,注意每一行都需要以;
结尾。 -
下一步根据基因名称定义数据的 partition ,命令行是
partition favored = 4: 28S, H3, 16S, COI;
。favored
为分区结构名称,=
连接分区结构的数量(此处是4个分区),冒号:
后则是每个分区的名称列表,名称用逗号,
分隔,最后不要忘记;
。 -
最后,命令软件依照定义的分区结构名称
favored
运行,而不是根据默认的。即定义set partition = favored;
。
这时,输入的序列文件即准备好了。为了一次性方便地设置模型参数,可以在序列数据下方,紧接着定义一些模型参数的命令行文本。
模型参数命令行
这里即对序列数据的不同分区定义不同的碱基(实为序列的位点字符)替换模型。以示例数据来说,28S 序列的替换模型为 GTR+G
,而其余三个分区的序列替换模型均为 GTR+I+G
,因此用 lset
命令设置各自分区的替换模型,具体就是 lset applyto=(1) rates=gamma nst=6;
和 lset applyto=(2,3,4) rates=invgamma nst=6;
。其中:
-
lset
定义似然模型参数,要计算似然性,就必须假设碱基替换模型,在 MrBayes 中输入命令行help lset
即可查看相应帮助文件。 -
applyto
定义分区的范围,即将模型应用到对应的分区,在括号()
中列出应用相同模型的分区。 -
nst
定义位点替换的模型(实际可以理解为 AC、AT、AG、GC、GT、TC 每对碱基的替换),用不同数字代替不同模型:1
代表所有替换率相同(如模型 JC69 或 F81);2
代表转换和颠换可以有不同的替换率(如模型 K80 或 HKY85);6
代表所有替换率都不同(如模型 GTR)。 -
rates
定义位点之间的替换率,有以下几种选择:equal
(位点的替换率无差异)、gamma
(位点的替换率呈 gamma 分布)、adgamma
(位点的替换率自相关,边缘位点替换率呈 gamma 分布,相邻位点有相关的替换率)、propinv
(一定比例位点的替换率是恒定的)、invgamma
(一定比例位点的替换率是恒定的,剩下位点的替换率呈 gamma 分布)。
接下来, unlink
命令是让每个分区序列所使用的模型都不连锁,每个分区单独应用各自的模型或参数。一般情况下,分区是根据形态或核酸等数据类型进行的,也可根据自定义的分区(如前述 set partion=<name/number>
)。软件运行时,如果应用于不同分区的参数有相同的先验,那么软件将会将这些参数连锁起来,也就是说,软件将会使用同一个参数,软件默认运行就是这样“吝啬”。但是,如果你想对各个分区使用不同的参数,比如,设置不同分区有不同的转换或颠换比率,那么就需要使用这一 unlink
命令。可在 MrBayes 中输入命令行 help unlink
查看相应帮助文件及一系列参数。在本例中,设置 unlink revmat=(all) pinvar=(all) shape=(all) statefreq=(all)
。
那么,每个分区的序列到底适用什么替换模型、什么替换率呢?即怎样选择各分区的替换模型呢?事实上,有很多软件可解决这一问题,如 jmodeltest 、BEAST 系列软件自带的 bModelTest 、以及 PartionFinder 等。个人认为 MrBayes 的模型选择使用 PartionFinder 较为方便,具体方法详见附章 1。
先验参数设定
最后, prset
命令可定义一系列系统发育模型的先验,设置先验参数分布是贝叶斯分析的必不可少的一步,先验数据代表观测数据之前的预先判断。详细帮助文件可在 MrBayes 中输入命令行 help prset
查看。本文也具体对此作了中文翻译说明,详见 附章 2。常用的先验设置主要有 revmatpr
, shapepr
, pinvarpr
, topologypr
, brlenspr
, popvarpr
, clockratepr
, clockvarpr
, ratepr
。
在本例中,applyto=(all)
表示对所有分区应用相同的参数设定。 ratepr=variable
表示设定分区的分子钟进化速率是可变的,即各分区进化速率不同,平均进化率是 1 。
运行软件
完成上述步骤后,即可运行软件,这里以 Mac 系统为例。对于已正确安装的 MrBayes ,只需在终端中输入 mb
,即可调出 MrBayes 软件。
-
导入
.nex
文件execute /path/to/finename.nex
-
其他模型和先验参数均已在
.nex
文件中设置 -
设置运行长度参数
mcmcp ngen=10000000 samplefreq=1000 diagnfreq=1000 # 可使用 help mcmcp 查看具体参数的意义 # 一般情况下,需要参考以往研究中使用的运行长度参数 # diagnfreq 默认为 5000,适用于较大的数据分析。该参数对输出的分割频率平均标准差( Average standard deviation of split frequencies)大小影响较大 # 输出的分割频率平均标准差要小于 0.05 ,才意味着能得到较可信的系统树
-
运行
mcmc
-
统计取样参数
sump
-
统计树取样
sumt
-
得到的
.tree
文件即可用 FigTree 等软件查看系统树。
附章 1 —— PatitionFinder 软件用法
PartionFinder 软件是由 Rob Lanfear 开发的,用于系统树构建之前、对序列的不同分区(如不同基因或不同类型数据)的适用替换模型进行判定和选择的软件。该软件基于 Python 开发,用法简便,结果稳定。详细操作说明参见 Manual。
下载 PartitionFinder v2.1.1 安装包 之后解压即用,当然前提是事先在系统中已安装 Python (注意,PartitionFinder v2.1.1 版本是基于 Python 2.7 编写的)。使用时主要需要编辑两个文件、并放在同一文件夹:
-
编辑
.phy
格式的序列数据文件。该格式文件与.nex
格式文件非常相似,区别是.phy
文件中的各样本的序列不能分割,并且第一行要写明样本数和序列长度。文件格式说明详见 PHYLIP 格式。就本示例数据来说,要将4个基因的序列根据样本合并成连续序列sequences_file_name.phy
,具体如下:2 225 Angulyagra_sp._Nam_Dinh_Angu.sp.1 CCGCTGAATTTAAGCATATCACTAAGTATGGCTCGTACCAAGCAGACNGCTCGTAAATCAACCG-GAGGAAAAGC--------------------------------------------------TAACCTGCCCGGTGA----------------------------------------------------------------GCATCCTGAGGTTTATATTTT Bellamya_aeruginosa_Baiyang_Lake_BY1 CCGCTGAATTTAAGCATATCACTAAGTATGGCTCGTACCAAGCAGACGGCTCGCAAATCGACCG-GAGGAAAAGCCGCCTGTTTATCAAAAACATGTCTTTTTGATTTTAGAAATATAAAAAGTCTGACCTGCCCGGTGAGGTCAACAAATCATAAAGATATTGGAACTTTATATATTTTGTTTGGTGTATGGTCAGGATTAGTTGGTACTGGACTTAGTATTTT
-
编辑
partition_finder.cfg
文件。该文件是 PartitionFinder 的运行文件,示例可参见 PartitionFinder 软件包中的exmples
文件夹中的各类型序列的示例。这里想要找出适用于 MrBayes 的模型,具体运行文件代码如下:## ALIGNMENT FILE ## alignment = sequences_file_name.phy; ## BRANCHLENGTHS: linked | unlinked ## branchlengths = unlinked; ## MODELS OF EVOLUTION: all | allx | mrbayes | beast | gamma | gammai | <list> ## models = mrbayes; # MODEL SELECCTION: AIC | AICc | BIC # model_selection = aicc; ## DATA BLOCKS: see manual for how to define ## [data_blocks] 28S = 1-26; H3 = 27-75; 16S = 76-140; COI = 141-225; ## SCHEMES, search: all | user | greedy | rcluster | rclusterf | kmeans ## [schemes] search = greedy;
-
将编辑好的
sequences_file_name.phy
序列文件和partition_finder.cfg
运行文件放到同一个文件夹中。
这样,数据运行的准备工作就完成了。打开终端(Mac)或运行(Windows),输入 python
,键入空格,将 PartitionFinder 安装包中的 PartitionFinder.py
拖入,键入空格,再将包含序列文件和运行文件的文件夹拖入,按下 Enter 键,即开始运行。
python /path/to/PartitionFinder/installpackage/PartitionFinder.py /path/to/the/folder
等待运行完成后,上述文件夹中即保存了运行结果。主要的可用运行结果在子文件夹 analysis/schemes/start_scheme.txt 文件中,文件内容包括 Nexus 文件格式、 RaxML 中使用的各分区的模型,以及 MrBayes 中使用的各分区的适用模型。下游进行 MrBayes 重建系统树时,即可将其复制到 Nexus 输入文件中,但是要注意,原封不动直接复制到 MrBayes 的 nex 文件中有可能会出错,要把 prset 的设置放在命令行的最后。
关于运行文件中的各个参数设置,简要说明如下:
-
Alignment file
即.phy
格式的序列文件。 -
Branch lengths
即分枝长度是否连锁,一般对于 MrBayes 和 BEAST 软件来说,都需要选择unlinked
。 -
Model of evolution
即各种待检验的进化模型。 PartitionFinder 软件很方便地直接将用于 MrBayes 的各种进化模型集成在一块,因此可直接设定models = mrbayes
,即检验序列适用哪一个 MrBayes 中的进化模型。 BEAST 也同理。如果想要自定义待检验的模型列表(模型之间用逗号隔开),可参见 Manual 对每个进化模型的详细说明。 -
Model selection
即模型选择的依据指标,一般情况下有 AICc 时就不要选择 AIC 了。在其操作说明中坦白指出了 AIC 指数只是由于历史原因在此保留了。 -
Data blocks
普通 DNA 序列时,该部分与 MrBayes 的 nexus 文件中的分区类似。对于密码子序列,还可单独对密码子的位置分区。 -
Schemes
如果序列数据有 12 个分区以上,就不能选择all
。一般情况下,选择greedy
,即告诉软件使用贪婪算法(greedy algorithm)寻找较好的分区体系(good partitioning scheme)。kmeans
参数被证明效果较差。其余三种参数仅适用于 raxml 的各种进化模型。 -
注意,运行文件的格式是是固定的,除了
#
号内的内容可随意更改,其余部分格式不可随意更改。 -
注意,运行文件中每一命令行都要以分号
;
结尾。
附章 2 —— 先验参数设定
先验参数设定比较繁琐复杂,其命令为 prset
,具体选项有:
-
applyto
该选项可将prset
命令应用到不同组的分区,并随之可设定各组对应的参数。如下例,第一行是设置前两个分区的参数,第二行是设置后两个分区的参数。prset applyto=(1,2) statefreqs=fixed(equal) prset applyto=(3,4) statefreqs=dirichlet(1,1,1,1)
-
tratiopr
该选项可设定转换或颠换的比率,用法如下。当选择beta
时,程序假设转换和颠换率是具有相同范围参数的独立的 gamma 分布随机变量。如果想要该 gamma 分布先验具有在1上下的相同转换/颠换率,则可以用平坦的beta
,即默认的beta(1,1)
。如果想要将更多概率聚集到 转换率在1的地方,则可用beta(x,x)
, x 就表示聚集的程度。比如,beta(20,20)
比起beta(1,1)
就有更多概率放在了转换率接近1的地方。如果认为转换率/颠换率的比值是2,那么就可设定beta(2,1)
,而beta(2,1)
就比beta(80,40)
分布得更发散。还可以用计数来设定,如有 x 个转换、 y 个颠换,那么就可设定beta(x+1, y+1)
。固定fixed
选项表示将 tratio 设定为一个固定的值。prset tratiopr = beta(<number>, <number>) prset tratiopr = fixed(<number>)
-
revmatpr
该选项设定核酸数据 GTR 模型的替换率先验,用法如下。当选择dirichlet
时,程序即认为 6 种碱基替换率是具有相同范围参数的独立 gamma 分布随机变量。括号中可设置6个数字,其顺序对应各对碱基的替换率:A<->C, A<->G, A<->T, C<->G, C<->T 以及 G<->T。默认设置是dirichlet(1,1,1,1,1,1)
,即一种“平坦”的分布。如果A<->C和A<->T的替代率分别是其他替代率的2倍和5倍,那么就可设置dirichlet(2x,x,5x,x,x,x)
,这里的 x 就代表先验集中在每个替代率的程度。当选择fixed
时表示设置替代率为某一固定值。prset revmatpr = dirichlet(<number>,<number>,...,<number>) prset revmatpr = fixed(<number>,<number>,...,<number>)
-
revratepr
该选项在nst
设置为mixed
时(参见lset
命令),可设定 GTR 模型的每个碱基对替换率的先验。用法如下。该选项可用一种修正的对称 Dirichlet 先验,将每个独立替换率联系成一个替换率的矩阵,而应用到 n 个成对碱基型的替换率有一个 alpha 值,表示该 alpha 是特定值的 n 倍。特定值越大,先验就越聚集到等比率位置,特定值的默认值是1,跟“平坦” Dirichlet 类似。prset revratepr = symdir(<number>)
-
statefreqpr
该参数设置 state (序列字符) 频率的先验。用法如下。dirichlet
可定义单个值,或有多少 state 就定义多少个值。如果定义的是单个值,则所有阶段都有相同的概率,方差则与给定的值有关。prset statefreqpr = dirichlet(<number>) prset statefreqpr = dirichlet(<number>,<number>,...,<number>) prset statefreqpr = fixed(equal) prset statefreqpr = fixed(empirical) prset statefreqpr = fixed(<number>,<number>,...,<number>)
-
shapepr
该参数设定位点替换率方差的 gamma shape 参数的先验。用法如下。( gamma 分布有两种参数,形状 shape 参数和尺度 scale 参数,形状参数代表 gamma 分布曲线的“胖瘦”,尺度参数代表 gamma 分布曲线的“高矮”)。prset shapepr = uniform(<number>,<number>) prset shapepr = exponetial(<number>) prset shapepr = fixed(<number>)
-
pinvarpr
该参数设定恒定位点比例的先验。用法如下。注意,参数的有效范围从 0 到 1 ,因此prset pinvarpr=uniform(0, 0.8)
是有效的,而prset pinvarpr=uniform(0, 10)
是无效的。默认设置是prset pinvarpr=uniform(0, 1)
。prset pinvarpr = uniform(<number>,<number>) prset pinvarpr = fixed(<number>)
-
ratecorrpr
该参数为位点替换率方差设定其自相关 gamma 分布的自相关系数的先验。用法如下。注意,参数值范围从 -1 到 1 ,默认设置为prset ratecorrpr = uniform(-1,1)
。prset ratecorrpr = uniform(<number>,<number>) prset ratecorrpr = fixed(<number>)
-
covswitchpr
该参数设定共变转换率(covarion switching rates)的先验。用法如下。共变模型(covarion model)有两个转换率:从 on 到 off 的转换率、从 off 到 on 的转换率。这些转换率被假定为独立的均匀或指数分布。还可以设定为固定的转换率,这时必须把 2 个转换率都定义。第一个值是 off -> on ,第二个值是 on -> off 。prset covswitchpr = uniform(<number>,<number>) prset covswitchpr = exponential(<number>) prset covswitchpr = fixed(<number>,<number>)
-
topologypr
该参数定义系统发育概率的先验。用法如下。当设定为默认的uniform
时,则所有可能树(possible tree)被认为有相同的先验。当设定为speciestree
时,则系统树的拓扑结构与其他(基因)树一起被折叠成物种树。当设定为constraints(<list>)
时,则可定义树的复杂先验概率(可以输入help constraint
获取更多信息),注意必须要输入一个待遵从的 constraint 列表。当设定为fixed
时,则表示固定用户指定的树,这时分枝长度仍会照常取样。prset topologypr = uniform prset topologypr = speciestree prset topologypr = constraints(<list>) prset topologypr = fixed(<treename>)
-
nodeagepr
该参数设定树内的与端节点和中间节点的年龄有关的假设,该参数只能用于 clock 树。默认模型是nodeagepr = unconstrained
,其假设所有端节点的年龄都相同而中间节点是不受约束的。另一种模型是nodeagepr = calibrated
,定义在校正设定下(详见calibrate
命令)的端节点和中间节点的概率分布先验。
示例:节点校正(node dating)。即已知某些序列 taxa 之间的节点年龄,根据这些已知的节点年龄,确定树的进化时间(dating)(用法可参见 MrBayes Manual v3.2,node-dating,p69,p73)。
参考 Manual MrBayes v3.2 中的 4.7 Node dating and total-evidence dating。
constraint humanchimp = Homo_sapiens Pan
# 第一步:定义某一个(一些)目标节点的 constraint 名称
calibrate humanchimp = uniform(5, 7)
# 第二步:给 constraint 节点赋一个校正时间的分布
prset topologypr=constraints(humanchimp)
# 第三步:强制系统树的 constraint 节点
prset nodeagepr=calibrated
# 第四步:强制节点年龄遵守校正 calibrated
-
popvarpr
在基因树——物种树中,该参数决定全部物种树中的种群大小是否恒定。默认设定是种群大小恒定popvarpr = equal
。设定物种树的种群大小是变化的popvarpr = variable
。prset popvarpr = equal prset popvarpr = variable
-
popsizepr
该参数设定coalescent
参数的种群大小成分的先验。用法如下。该参数只应用于当分枝长度先验选择 coalescent 模型的时候。注意,lset
命令中的ploidy
参数对解释该参数非常重要。prset popsizepr = uniform(<number>,<number>) prset popsizepr = exponential(<number>) prset popsizepr = fixed(<number>)
-
brlenspr
该参数设定分枝长度的分布概率的先验。用法如下。其中unconstrained
分枝的树是无根树,而clock-constraint
的树是有根树。clock
后面的参数可以定义详细的分枝长度的密度概率,当选择birthdeath
或coalescence
先验时,可能要修改这些过程的详细参数(如 birthdeath 先验中的物种形成率、灭绝率、取样概率, coalescence 先验中的种群大小、分子钟率参数)。选择clock:speciestreecoalescence
时,则基因树被折叠成物种树,该模型下,使用命令popvarpr
可设定种群大小是恒定还是变化。用fixed
分枝长度则被固定。prset brlenspr = unconstrained:uniform(<num>,<num>) prset brlenspr = unconstrained:exponential(<number>) prset brlenspr = clock:uniform prset brlenspr = clock:birthdeath prset brlenspr = clock:coalescence prset brlenspr = clock:speciestreecoalescence prset brlenspr = fixed(<treename>)
-
clockratepr
该参数设定树中与碱基替换率有关的假设先验,即每个位点在每单位时间里碱基替换的期望数目,注意,该参数仅对分子钟树有效。默认设置是固定值 0 即Fixed(1.0)
,这时实际上单位时间就等于每个位点碱基替换的期望数目。如果将对数应用在限制年龄(age constraints)上,那么默认设置就会自动变成Exponential(<x>)
,这里的<x>
可设置成比如指数的期望值是最大限制年龄的 10 倍,这会是一个比较模糊的先验,不足以解决问题。如果树中没有任何年龄校准,仍然可以使用clockratepr
来校准树。比如,如果已知某一基因序列的替代率是每百万年 0.20 ,则可以将替换率固定成 0.20 ,即prset clockratepr = fixed(0.20)
,这样树就会以百万年为单位进行校准。在不知道确切的替换率值的时候,还可以给替换率指定一个概率分布,可以选择正态分布、对数分布、指数分布和 gamma 分布。比如,如果想得出截距为 0 的关于替换率的正态分布,只允许正值,平均值为 0.20 ,标准差为 0.02 ,则可设置prset clockratepr = normal(0.20,0.02)
。自然对数分布中使用的参数是根据其平均值和标准差而来的,如prset clockratepr = lognormal(-1.61,0.10)
,表示该自然对数分布的平均 log 值是 -1.61 ,标准差是 0.10 ,在这种情况下整个对数分布的平均值将是e^(-1.61 + 0.10^2/2) = 0.20
。prset clockratepr = fixed(0.20) prset clockratepr = normal(0.20,0.02) prset clockratepr = lognormal(-1.61,0.10)
-
treeagepr
当 clock tree 的分支长度使用了brlenspr = clock:uniform
先验、并且clockratepr
不再使用默认值 1 时,就需要使用该命令设定 tree age 的分布概率先验。用法如下。 tree age 简单来讲即树中最近共同祖先的年龄。如果分子钟率(clock rate)设定为默认即固定值 1 ,则树年龄 tree age 就等于从树根部到末端的替换数量的期望值(the expected number of substitutions from the root to the tip of the tree),这时treeagepr
设置为默认的exponential(1)
一般都能运行良好。假如分子钟率设定为0.001,那么用期望均值 1 除以该分子钟率 0.001 ,即 1/0.001=1000,这时treeagepr
设置就可依据 1000 的倒数即 0.001 设置到exponential()
的括号中(MrBayes Manual v3.2, p66)。树年龄先验可以确保分枝长度的 uniform 先验模型的节点先验概率分布是正确的(proper)。默认设置是Exponential(1)
。注意,如果树的根节点已被校正,那么根校正(root calibration)设置之后,就不必再设置树年龄先验(tree age prior)treeagepr
了(用法可参见 MrBayes Manual v3.2 node-dating 和 total-evidence dating 中的用法,p73)。prset treeagepr = Gamma(<number>,<number>) prset treeagepr = Exponential(<number>) prset treeagepr = Fixed(<number>)
-
clockvarpr
该参数可设定假设的分子钟类型,注意,该参数仅对分子钟树有效。。默认是strict
,即标准分子钟模型,其假设整个树的进化率是恒定不变的。还可设置为cpp
,是一种 relaxed 分子钟模型,进化率是依据 Huelsenbeck et al. 2000 提出的复合泊松过程(Compound Poisson Process, CPP)模型。还可设置为tk02
,依据的是 Thorne & Kishino 2002 提出的 Brownian Motion 模型。还可以是另一种模型,该模型的每个分枝有一个独立的、符合 scaled gamma 分布的进化率,这样在先验中就要设定树的有效高度(the effective height of the tree)的差异,这就是一个 IGR (Independent Gamma Rate) 模型(LePage et al. 2007)。上述每个 ralaxed 分子钟模型都有额外的先验参数设置, CPP 模型有cppratepr
和cppmultdevpr
, TK02 模型是tk02varpr
, IGR 模型是igrvarpr
。在此处,bm
跟tk02
是同义的,ibr
跟igr
是同义的。 -
cppratepr
即当使用 CPP relaxed 分子钟模型时、进化率改变过程中生成的泊松过程(Poisson process)比率的先验概率分布。既可以设置为固定值,也可设置为指数分布,用法如下。比如,如果将进化速率设置为固定值 2 ,则在每个位点的平均预期替换次数等于1的长度的分枝上,平均期望可看到 2 个速率修改事件。如果将进化速率设置成指数分布exponential(0.1)
,则在期望比率为 10 (= 1/0.1) 的先验概率分布中估算进化速率。prset cppratepr = fixed(<number>) prset cppratepr = exponential(<number>)
-
cppmultdevpr
该参数可设置当使用 CPP relaxed 分子钟的进化率乘数时、对数分布的标准差。标准差以对数尺度给出。默认值是 1.0 ,因此当期望值加减一个标准差时、对应的速率乘数(rate multiplier)从 0.37 (1/e) 到 2.7 (e) 变化。乘数(multiplier)对数的期望平均值为 0 ,以确保期望平均速率是 1.0 。可使用下面用法改变默认设置,括号中的数值就是 log scale 的标准差。prset cppmultdevpr = fixed(<number>)
-
tk02varpr
该参数可定义 Thorne-Kishino (‘Brownian motion’) relaxed 分子钟模型下进化速率乘数差异的先验概率分布。特殊的是,其设定的是随着基本分子钟速率增加的方差下的参数。如果根据基本分子钟速率、有一个分枝长度对应为每位点期望替换值为 0.4 ,且tk02var
的参数值为 2.0 ,则分枝末端的速率乘数就会是一个方差为0.4*2.0
(线性的,不是log scale)的对数分布。其意义与分枝始端的速率乘数相同。也可将参数设置为固定,或设定为指数或均匀分布。此处,bmvarpr
和tk02varpr
是同义的。prset tk02varpr = fixed(<number>) prset tk02varpr = exponential(<number>) prset tk02varpr = uniform(<number>,<number>)
-
igrvarpr
该参数可设置在独立分枝进化率(IGR)relaxed 分子钟模型下分枝长度 gamma 分布的方差的先验。特别的是,特殊的是,其设定的是随着基本分子钟速率增加的方差下的参数。如果根据基本分子钟速率、有一个分枝长度对应为每位点期望替换值为 0.4 ,且igrvar
的参数值为 2.0 ,则分枝末端的速率乘数就会是一个方差为0.4*2.0
。也可将参数设置为固定,或设定为指数或均匀分布。此处,ibrvarpr
和igrvarpr
是同义的。prset igrvarpr = fixed(<number>) prset igrvarpr = exponential(<number>) prset igrvarpr = uniform(<number>,<number>)
-
ratepr
该参数可设置位点特异进化率模型或其他模型,这些模型允许不同分区以不同的速率进化。首要的是定义序列分区,比如DNA序列不同基因定义分区。定义好分区之后,就可以设定不同分区的不同速率乘数。用法如下。当设定为fixed
时,则该分区的速率乘数就设定为 1 (比如各分区平均进化速率)。当设定为variable
时,则服从限制的分区的进化率是不同的,平均进化率是 1 ,必须对至少两个分区设定不同的进化率先验,否则该参数在计算时就无效。variable
其实就是dirichlet(1,...,1)
先验。dirichlet
可设定不同分区的不同速率,并可精确控制先验的形状 shape 。dirichlet
参数中各值的顺序跟应用进去的分区的顺序一致,比如prset applyto=(1,3,4) ratepr = dirichlet(10,40,15)
就表示分区1与 10 对应,分区2与 40 对应,分区3与 15 对应。dirichlet
参数就是用来权衡各个进化速率的,就是说它根据每个分区中包含的字符数加权分区率。prset ratepr = fixed prset ratepr = variable prset ratepr = dirichlet(<number>,<number>,...,<number>)
-
speciationpr
该参数设定 net 物种形成率的先验,即 birth-death 模型中的Φ - μ
(lambda - mu)。用法如下。该参数只应用于分枝长度先验选择 birth-death 模型的时候。prset speciationpr = uniform(<number>,<number>) prset speciationpr = exponential(<number>) prset speciationpr = fixed(<number>)
-
extinctionpr
该参数设定相对灭绝率先验,即 birth-death 模型中的μ/Φ
。该参数的值范围是 0 到 1。用法如下。该参数只应用于分枝长度先验选择 birth-death 模型的时候。prset extinctionpr = beta(<number>,<number>) prset extinctionpr = fixed(<number>)
-
aamodelpr
该参数设定氨基酸数据的替换率矩阵。可以设定固定值替换率aamodelpr = fixed(<model name>)
,这里的<model name>
可以是poisson
(一种优化的 Jukes-Cantor 模型),jones
,dayhoff
,mtrev
,mtmam
,wag
,rtrev
,cprev
,vt
,blosum
,equalin
(一种优化的 Frlsenstein 1981 模型), 或gtr
。也可以通过设定aamodelpr=mixed
来将前十个模型平均化,这时马尔科夫链将根据每个模型的概率来取样。前十个模型poisson(0)
,jones(1)
,dayhoff(2)
,mtrev(3)
,mtmam(4)
,wag(5)
,rtrev(6)
,cprev(7)
,vt(8)
,blosum(9)
可以用数字编码 0-9 代替。用命令sump
即可总结 MCMC 取样并计算每个模型的后验概率估算值。 -
aarevmaptr
该参数设定氨基酸序列 GTR 模型的替换率先验。用法与revmatpr
类似:prset aarevmatpr = dirichlet(<number>,<number>,...,<number>) prset aarevmatpr = fixed(<number>,<number>,...,<number>)
-
omegapr
该参数在当核苷酸替代模型设定为密码子即lset nucmodel=codon
时,可设定其同义或非同义替换率的比例。并且,该参数只能应用在当位点的 omega 无差异的时候即lset omegavar=equal
。用法如下:prset omegapr = uniform(<number>,<number>) prset omegapr = exponential(<number>) prset omegapr = fixed(<number>)
-
Ny98omega1pr
该参数在当核苷酸替代模型设定为密码子即lset nucmodel=codon
时,可设定被选位点(sites under purifying selection)的同义或非同义替换率的比例。并且,该参数只能应用在Nielsen & Yang (1998) 模型下位点的 omega 有差异的时候即lset omegavar=ny98
。当设定为固定值fixed
的时候,括号内的值要大于0小于1。用法如下:prset Ny98omega1pr = beta(<number>,<number>) prset Ny98omega1pr = fixed(<number>)
-
Ny98omega3pr
该参数在当核苷酸替代模型设定为密码子即lset nucmodel=codon
时,可设定确切被选位点(positively selected sites)的同义或非同义替换率的比例。并且,该参数只能应用在 Nielsen & Yang (1998) 模型下位点的 omega 有差异的时候即lset omegavar=ny98
。用法如下。注意,当定义为Ny98
模型时参数值必须大于1,比如不能设定为uniform(0,10)
。prset Ny98omega3pr = uniform(<number>,<number>) prset Ny98omega3pr = exponential(<number>) prset Ny98omega3pr = fixed(<number>)
-
N3omegapr
该参数在当核苷酸替代模型设定为密码子即lset nucmodel=codon
时,可设定 M3 模型下全部3类位点的同义或非同义替换率的比例。并且,该参数只能应用在 Yang et al. (2000) 的 M3 模型下位点的 omega 有差异的时候,如设定lset omegavar=M3
。用法如下。在该指数先验下,四种替换率dN1
,dN2
,dN3
,dS
均呈独立的指数分布(无需手动定义指数分布)。dN1
、dN2
、dN3
这三种替换率的顺序是dN1
<dN2
<dN3
,它们都以同一个同义替换率dS
为尺度。prset M3omegapr = exponential prset M3omegapr = fixed(<number>,<number>,<number>)
-
codoncatfreqs
该参数在当核苷酸替代模型设定为密码子即lset nucmodel=codon
时,可设定在 purifying, neutral, 和 positive selection 下位点的频率的先验。并且,该参数只能应用在 Nielsen & Yang (1998) 模型或 Yang et al. (2000) 的 M3 模型下位点的 omega 有差异的时候,如设定lset omegavar=Ny98
或lset omegavar=M3
。用法如下。注意括号中 3 个频率的和必须是 1 。prset codoncatfreqs = dirichlet(<number>,<number>,<number>) prset codoncatfreqs = fixed(<number>,<number>,<number>)
-
symdirihyperpr
该参数设定形态学数据 the stationary frequencies of the states 的先验。形态学数据有多达 10 个 states(形态特征代码),但是这些 states 的标注有时都很武断,比如 1 在不同的形态特征中代表不同的内容,这一点跟 DNA 的碱基字母不同。软件假设这些 states 都在 dirichlet 先验下有相同的频率,设定不同的 dirichlet 先验就能用该参数实现。用法如下。当设定为fixed(infinity)
时, dirichlet 先验则是固定值,所有的字符都有相同的频率。prset Symdirihyperpr = uniform(<number>,<number>) prset Symdirihyperpr = exponential(<number>) prset Symdirihyperpr = fixed(<number>) prset Symdirihyperpr = fixed(infinity)
-
samplestart
该参数设置在分析中对物种进行采样的策略。用于 birth-death 模型树中(参见 Höhna et al 2011)。 -
sampleprob
该参数设定在分析中被采样的物种的比例。用于 birth-death 模型树中(参见 Yang & Rannala 1997)。