yzy2020的个人博客分享 http://blog.sciencenet.cn/u/yzy2020 技术体现的是术,可以通过多次练习掌握,不要迷恋术,idea是道,需要通过文献加强训练。总之,孰能生巧!

博文

[转载]MrBayes 操作说明

已有 1411 次阅读 2023-5-20 15:55 |个人分类:linux学习|系统分类:科研笔记|文章来源:转载

https://bin-ye.com/post/2019/10/19/%E5%A5%BD%E5%A5%BD%E5%85%88%E7%94%9F-mrbayes-%E6%93%8D%E4%BD%9C%E8%AF%B4%E6%98%8E/

好好先生 —— 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 后面的参数可以定义详细的分枝长度的密度概率,当选择 birthdeathcoalescence 先验时,可能要修改这些过程的详细参数(如 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 模型有 cpprateprcppmultdevpr , TK02 模型是 tk02varpr , IGR 模型是 igrvarpr 。在此处, bmtk02 是同义的,ibrigr 是同义的。

  • 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)的对数分布。其意义与分枝始端的速率乘数相同。也可将参数设置为固定,或设定为指数或均匀分布。此处, bmvarprtk02varpr 是同义的。

    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 。也可将参数设置为固定,或设定为指数或均匀分布。此处, ibrvarprigrvarpr 是同义的。

    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 均呈独立的指数分布(无需手动定义指数分布)。dN1dN2dN3 这三种替换率的顺序是 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=Ny98lset 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://bin-ye.com/post/2019/10/19/%E5%A5%BD%E5%A5%BD%E5%85%88%E7%94%9F-mrbayes-%E6%93%8D%E4%BD%9C%E8%AF%B4%E6%98%8E/




https://blog.sciencenet.cn/blog-3434047-1388735.html

上一篇:[转载]MrBayes 操作说明
下一篇:conda create -n py27 python=2.7.10 #安装“package”并创建特定环境
收藏 IP: 221.11.67.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-3-28 20:12

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部