博文
[转载]MrBayes 操作说明
|||
好好先生 —— MrBayes 操作说明
Tags: summary
由 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）。
https://blog.sciencenet.cn/blog-3434047-1388734.html
上一篇：[转载]使用 ProtTest 来选择最优氨基酸替代模型
下一篇：[转载]MrBayes 操作说明
