365bet网址搜索器

好好先生 —— MrBayes 操作说明

发布时间: 2025-09-19 06:31:14 作者: admin 阅读量: 9007 评论数: 258

由 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=)。软件运行时,如果应用于不同分区的参数有相同的先验,那么软件将会将这些参数连锁起来,也就是说,软件将会使用同一个参数,软件默认运行就是这样“吝啬”。但是,如果你想对各个分区使用不同的参数,比如,设置不同分区有不同的转换或颠换比率,那么就需要使用这一 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 | ##

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(, )

prset tratiopr = fixed()

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(,,...,)

prset revmatpr = fixed(,,...,)

revratepr 该选项在 nst 设置为 mixed 时(参见 lset 命令),可设定 GTR 模型的每个碱基对替换率的先验。用法如下。该选项可用一种修正的对称 Dirichlet 先验,将每个独立替换率联系成一个替换率的矩阵,而应用到 n 个成对碱基型的替换率有一个 alpha 值,表示该 alpha 是特定值的 n 倍。特定值越大,先验就越聚集到等比率位置,特定值的默认值是1,跟“平坦” Dirichlet 类似。

prset revratepr = symdir()

statefreqpr 该参数设置 state (序列字符) 频率的先验。用法如下。 dirichlet 可定义单个值,或有多少 state 就定义多少个值。如果定义的是单个值,则所有阶段都有相同的概率,方差则与给定的值有关。

prset statefreqpr = dirichlet()

prset statefreqpr = dirichlet(,,...,)

prset statefreqpr = fixed(equal)

prset statefreqpr = fixed(empirical)

prset statefreqpr = fixed(,,...,)

shapepr 该参数设定位点替换率方差的 gamma shape 参数的先验。用法如下。( gamma 分布有两种参数,形状 shape 参数和尺度 scale 参数,形状参数代表 gamma 分布曲线的“胖瘦”,尺度参数代表 gamma 分布曲线的“高矮”)。

prset shapepr = uniform(,)

prset shapepr = exponetial()

prset shapepr = fixed()

pinvarpr 该参数设定恒定位点比例的先验。用法如下。注意,参数的有效范围从 0 到 1 ,因此 prset pinvarpr=uniform(0, 0.8) 是有效的,而 prset pinvarpr=uniform(0, 10) 是无效的。默认设置是 prset pinvarpr=uniform(0, 1) 。

prset pinvarpr = uniform(,)

prset pinvarpr = fixed()

ratecorrpr 该参数为位点替换率方差设定其自相关 gamma 分布的自相关系数的先验。用法如下。注意,参数值范围从 -1 到 1 ,默认设置为 prset ratecorrpr = uniform(-1,1) 。

prset ratecorrpr = uniform(,)

prset ratecorrpr = fixed()

covswitchpr 该参数设定共变转换率(covarion switching rates)的先验。用法如下。共变模型(covarion model)有两个转换率:从 on 到 off 的转换率、从 off 到 on 的转换率。这些转换率被假定为独立的均匀或指数分布。还可以设定为固定的转换率,这时必须把 2 个转换率都定义。第一个值是 off -> on ,第二个值是 on -> off 。

prset covswitchpr = uniform(,)

prset covswitchpr = exponential()

prset covswitchpr = fixed(,)

topologypr 该参数定义系统发育概率的先验。用法如下。当设定为默认的 uniform 时,则所有可能树(possible tree)被认为有相同的先验。当设定为 speciestree 时,则系统树的拓扑结构与其他(基因)树一起被折叠成物种树。当设定为 constraints() 时,则可定义树的复杂先验概率(可以输入 help constraint 获取更多信息),注意必须要输入一个待遵从的 constraint 列表。当设定为 fixed 时,则表示固定用户指定的树,这时分枝长度仍会照常取样。

prset topologypr = uniform

prset topologypr = speciestree

prset topologypr = constraints()

prset topologypr = fixed()

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(,)

prset popsizepr = exponential()

prset popsizepr = fixed()

brlenspr 该参数设定分枝长度的分布概率的先验。用法如下。其中 unconstrained 分枝的树是无根树,而 clock-constraint 的树是有根树。clock 后面的参数可以定义详细的分枝长度的密度概率,当选择 birthdeath 或 coalescence 先验时,可能要修改这些过程的详细参数(如 birthdeath 先验中的物种形成率、灭绝率、取样概率, coalescence 先验中的种群大小、分子钟率参数)。选择 clock:speciestreecoalescence 时,则基因树被折叠成物种树,该模型下,使用命令 popvarpr 可设定种群大小是恒定还是变化。用 fixed 分枝长度则被固定。

prset brlenspr = unconstrained:uniform(,)

prset brlenspr = unconstrained:exponential()

prset brlenspr = clock:uniform

prset brlenspr = clock:birthdeath

prset brlenspr = clock:coalescence

prset brlenspr = clock:speciestreecoalescence

prset brlenspr = fixed()

clockratepr 该参数设定树中与碱基替换率有关的假设先验,即每个位点在每单位时间里碱基替换的期望数目,注意,该参数仅对分子钟树有效。默认设置是固定值 0 即 Fixed(1.0) ,这时实际上单位时间就等于每个位点碱基替换的期望数目。如果将对数应用在限制年龄(age constraints)上,那么默认设置就会自动变成 Exponential() ,这里的 可设置成比如指数的期望值是最大限制年龄的 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(,)

prset treeagepr = Exponential()

prset treeagepr = Fixed()

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()

prset cppratepr = exponential()

cppmultdevpr 该参数可设置当使用 CPP relaxed 分子钟的进化率乘数时、对数分布的标准差。标准差以对数尺度给出。默认值是 1.0 ,因此当期望值加减一个标准差时、对应的速率乘数(rate multiplier)从 0.37 (1/e) 到 2.7 (e) 变化。乘数(multiplier)对数的期望平均值为 0 ,以确保期望平均速率是 1.0 。可使用下面用法改变默认设置,括号中的数值就是 log scale 的标准差。

prset cppmultdevpr = fixed()

tk02varpr 该参数可定义 Thorne-Kishino (‘Brownian motion’) relaxed 分子钟模型下进化速率乘数差异的先验概率分布。特殊的是,其设定的是随着基本分子钟速率增加的方差下的参数。如果根据基本分子钟速率、有一个分枝长度对应为每位点期望替换值为 0.4 ,且 tk02var 的参数值为 2.0 ,则分枝末端的速率乘数就会是一个方差为 0.4*2.0 (线性的,不是log scale)的对数分布。其意义与分枝始端的速率乘数相同。也可将参数设置为固定,或设定为指数或均匀分布。此处, bmvarpr 和 tk02varpr 是同义的。

prset tk02varpr = fixed()

prset tk02varpr = exponential()

prset tk02varpr = uniform(,)

igrvarpr 该参数可设置在独立分枝进化率(IGR)relaxed 分子钟模型下分枝长度 gamma 分布的方差的先验。特别的是,特殊的是,其设定的是随着基本分子钟速率增加的方差下的参数。如果根据基本分子钟速率、有一个分枝长度对应为每位点期望替换值为 0.4 ,且 igrvar 的参数值为 2.0 ,则分枝末端的速率乘数就会是一个方差为 0.4*2.0 。也可将参数设置为固定,或设定为指数或均匀分布。此处, ibrvarpr 和 igrvarpr 是同义的。

prset igrvarpr = fixed()

prset igrvarpr = exponential()

prset igrvarpr = uniform(,)

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(,,...,)

speciationpr 该参数设定 net 物种形成率的先验,即 birth-death 模型中的 Φ - μ (lambda - mu)。用法如下。该参数只应用于分枝长度先验选择 birth-death 模型的时候。

prset speciationpr = uniform(,)

prset speciationpr = exponential()

prset speciationpr = fixed()

extinctionpr 该参数设定相对灭绝率先验,即 birth-death 模型中的 μ/Φ。该参数的值范围是 0 到 1。用法如下。该参数只应用于分枝长度先验选择 birth-death 模型的时候。

prset extinctionpr = beta(,)

prset extinctionpr = fixed()

aamodelpr 该参数设定氨基酸数据的替换率矩阵。可以设定固定值替换率 aamodelpr = fixed() ,这里的 可以是 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(,,...,)

prset aarevmatpr = fixed(,,...,)

omegapr 该参数在当核苷酸替代模型设定为密码子即 lset nucmodel=codon 时,可设定其同义或非同义替换率的比例。并且,该参数只能应用在当位点的 omega 无差异的时候即 lset omegavar=equal 。用法如下:

prset omegapr = uniform(,)

prset omegapr = exponential()

prset omegapr = fixed()

Ny98omega1pr 该参数在当核苷酸替代模型设定为密码子即 lset nucmodel=codon 时,可设定被选位点(sites under purifying selection)的同义或非同义替换率的比例。并且,该参数只能应用在Nielsen & Yang (1998) 模型下位点的 omega 有差异的时候即 lset omegavar=ny98 。当设定为固定值 fixed 的时候,括号内的值要大于0小于1。用法如下:

prset Ny98omega1pr = beta(,)

prset Ny98omega1pr = fixed()

Ny98omega3pr 该参数在当核苷酸替代模型设定为密码子即 lset nucmodel=codon 时,可设定确切被选位点(positively selected sites)的同义或非同义替换率的比例。并且,该参数只能应用在 Nielsen & Yang (1998) 模型下位点的 omega 有差异的时候即 lset omegavar=ny98 。用法如下。注意,当定义为 Ny98 模型时参数值必须大于1,比如不能设定为 uniform(0,10)。

prset Ny98omega3pr = uniform(,)

prset Ny98omega3pr = exponential()

prset Ny98omega3pr = fixed()

N3omegapr 该参数在当核苷酸替代模型设定为密码子即 lset nucmodel=codon 时,可设定 M3 模型下全部3类位点的同义或非同义替换率的比例。并且,该参数只能应用在 Yang et al. (2000) 的 M3 模型下位点的 omega 有差异的时候,如设定 lset omegavar=M3 。用法如下。在该指数先验下,四种替换率 dN1, dN2, dN3, dS 均呈独立的指数分布(无需手动定义指数分布)。dN1、 dN2、 dN3 这三种替换率的顺序是 dN1< dN2

prset M3omegapr = exponential

prset M3omegapr = fixed(,,)

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(,,)

prset codoncatfreqs = fixed(,,)

symdirihyperpr 该参数设定形态学数据 the stationary frequencies of the states 的先验。形态学数据有多达 10 个 states(形态特征代码),但是这些 states 的标注有时都很武断,比如 1 在不同的形态特征中代表不同的内容,这一点跟 DNA 的碱基字母不同。软件假设这些 states 都在 dirichlet 先验下有相同的频率,设定不同的 dirichlet 先验就能用该参数实现。用法如下。当设定为 fixed(infinity) 时, dirichlet 先验则是固定值,所有的字符都有相同的频率。

prset Symdirihyperpr = uniform(,)

prset Symdirihyperpr = exponential()

prset Symdirihyperpr = fixed()

prset Symdirihyperpr = fixed(infinity)

samplestart 该参数设置在分析中对物种进行采样的策略。用于 birth-death 模型树中(参见 Höhna et al 2011)。

sampleprob 该参数设定在分析中被采样的物种的比例。用于 birth-death 模型树中(参见 Yang & Rannala 1997)。

相关文章