使用BAMM估算分化速率

BAMM(Bayesian Analysis of Macroevolutionary Mixtures)是一种根据系统发生树对物种形成、灭绝和性状进化的复杂动力学建模的程序。

官网：http://bamm-project.org/
BAMM下载页面：http://bamm-project.org/download.html
说明书链接：http://bamm-project.org/quickstart.html#incomplete-taxon-sampling

win版本下载后只包含了.exe和.dll，需要放置在同一目录下使用。为了学习需要示例数据，下载源文件或仅下载文件列表里的example。https://github.com/macroevolution/bamm
首先将以下文件放在同一目录下：
1.分子钟估算后带时间标记的系统发生树
2.control文件（参数设置）
3.BAMM程序
在使用前用R的ape包检查树是否符合分析条件

library(ape)
v <- read.tree("myPhylogeny.tre")
is.ultrametric(v)
is.binary.tree(v)
# Now to check min branch length:
min(v$edge.length) 

树的枝长以百万年为单位

然后是control文件，类似于一个脚本，里面是对参数的设置，在example文件夹里有diversification和 trait，在diversification里找divcontrol.txt文件，打开可以看到里面的参数设定，每行都有注释解说。默认的设定对大多是百万年级别的树有用，根据情况需要调整来获取理想结果。先复制一份到BAMM所在目录下，这里也放出来看看，主要需要修改的是前面几个。

# BAMM configuration file for speciation/extinction analysis
# ==========================================================
#
# Format
# ------
#
#     - Each option is specified as: option_name = option_value
#     - Comments start with # and go to the end of the line
#     - True is specified with "1" and False with "0"


################################################################################
# GENERAL SETUP AND DATA INPUT
################################################################################

#设置模式
modeltype = speciationextinction
# Specify "speciationextinction" or "trait" analysis

#系统发生树文件名称（树文件要跟BAMM运行程序放在一个文件路径里）
treefile = whaletree.tre
# File name of the phylogenetic tree to be analyzed

#输出文件的名称
runInfoFilename = run_info.txt
# File name to output general information about this run

sampleFromPriorOnly = 0
# Whether to perform analysis sampling from prior only (no likelihoods computed)

#后略 

其中，先验概率模块可以用R的BAMMtools包自动生成

library(BAMMtools)
tr <- read.tree("treePathway")
setBAMMpriors(tr) 

运行后会生成一个叫做 myPriors.txt 的文本问件在工作目录下，内容如下图【输入getwd（）查询R目前的工作路径】

1689839234(1).png

将这些数据替换掉控制文件中相对应的数据

需要注意的是，这个生成的数据是一套对多个数据有效的设定根据树的时间枝长尺度进行缩放的结果，并不一定适用于所有数据。作者对于先验概率对数据的影响有较大篇幅的解说，看了一遍啥都没看懂，数学能力强的可以进一步了解。

将myPriors.txt 里的内容复制，替换control文件里的 # PRIORS 部分，检查一下后最好重命名一下control文件，如命名为 ‘ testcontrol ’。

这里用BioGeobears包里的Psychotria_5.2.newick树文件从头试一遍。

树文件放在BAMM所在目录下，将工作目录转过来比较方便。

打开Rstudio，输入

library(BAMMtools)
wd <- "C://BAMM所在路径//BAMM"
setwd（wd）
tr <- read.tree("Psychotria_5.2.newick")
setBAMMpriors(tr) 

得到myPriors.txt 内容如下，复制替换control文件，开头的树文件名称更改，保存关闭后重命名为 testcontrol.txt

（2023.3.3补：我看到有的文章先验设置了多个，1、2、3、4、5、10，都跑了一遍，结果一致。在不太好解释的情况下这样设置可能是比较好的选择，就怕不同的先验算出不同后验结果）

###############################################

# Prior block chosen by BAMMtools::setBAMMpriors

expectedNumberOfShifts = 1.0

lambdaInitPrior = 0.461956997597441

lambdaShiftPrior = 0.221402412787889

muInitPrior = 0.461956997597441

#### End Prior block
###################### 

打开win系统的命令提示符（不是那个bamm.exe, 是从电脑开始菜单里找的CMD）

图片.png

从cmd去调用bamm.exe，不是直接点开用的，直接点只会闪一下指令框

cd "C://BAMM所在路径//BAMM"
bamm -c testcontrol.txt 

再回到Rstudio里实行可视化

evdata <- getEventData(tr, eventdata = "event_data.txt",burnin = 0.1)
plot(eventdata,method='polar',lwd=3,pal='RdYlGn') 

incompsampling模式

如果是采样不充分的系统树，可以使用incompsampling模式，说明书链接为：http://bamm-project.org/advanced.html#incompsampling

对于较少不完整的系统发育：不完整的分类单元抽样可能会对系统发育树的物种形成和灭绝分析产生偏差。如果您有一个不完全采样的树（例如，您的目标分类群的总现存物种丰富度的某些分数 < 1），您可以轻松指定此采样分数以在物种缺失的假设下生成对物种形成和灭绝的无偏估计随机从树上。在您的控制文件中，您可以通过设置参数来指定已采样物种的百分比globalSamplingFraction。具体来说，您应该在控制文件中进行以下设置：

useGlobalSamplingProbability =  1
globalSamplingFraction = XYZ

其中 XYZ 是进化枝中包含在系统发育中的总分类群的百分比（例如，如果你的进化枝有 100 个物种而你的树有 73 个，你应该设置）。XYZ = 0.73

但是如果你的物种不是随机抽样的呢？BAMM 允许您将多个级别的此类非随机性合并到您的分析中。考虑以下示例。假设您希望对由两个属组成的进化枝执行多样化分析：fu和bar。在 fu分支中，您对某些物种 A 和 B 进行了采样，但另外 8 个物种仍未采样。在第二个属栏中，您对物种 C、D、E 和 F 进行了采样，并且仅缺少一个物种。您可以得出结论，您在fu中有 0.2 (2 / 10) 个物种，在bar中有 0.8 个物种（4 / 5）。此外，由于这些是您的焦点树中仅有的两个属，您知道系统发育的骨干具有 1.0 的采样概率（因为进化枝中没有包含 fu 和 bar 的其他属）。我们可以如下在这棵树上“绘制”相关的采样分数，分别使用红色、蓝色和黑色作为 fu 、 bar和backbone 采样分数。

图片.png

我们可以直接解释 BAMM 中这种类型的不完整抽样。我们需要做两件事。首先，我们告诉 BAMM 我们不会使用单一的全局抽样概率，而是指定进化枝或物种特定的抽样概率。在控制文件中，您设置. 然后，您必须指定一个文件名，该文件名使用参数给出相关的采样分数。这个文件的格式有点奇特。在文件的第一行，您给出主干采样概率（此处示例中为 1.0）。useGlobalSamplingProbability = 0sampleProbsFilename

现在，您为系统发育中的每个样本物种添加一条单独的线（本例中为 A、B、C、D、E、F）。每行包含以下内容，中间用tab符合隔开（不能用空格）：

speciesName     cladeName       samplingFraction

对于物种 A，相关行将如下所示：

A       fu      0.2

所有物种的完整文件如下所示：

1.0
A       fu      0.2
B       fu      0.2
C       bar     0.8
D       bar     0.8
E       bar     0.8
F       bar     0.8

对于许多数据集，这是使用不完整采样信息的一种简单方法。假设您正在为所有雀形目鸟类的多样化动态建模。如果您拥有鸟类的所有科，但这些科中的抽样不完整，则可以直接合并此抽样，前提是您可以合理地假设单系科。

如果主干抽样分数不太可能等于 1.0，则指定主干抽样分数可能会带来一些挑战。这是一个不等于 1 的示例。继续上面的例子，假设在我们的焦点进化枝中还有其他几个总计 7 个物种没有被采样。前提是我们（合理地）确信该属不属于我们的属fu和bar，这意味着这些谱系必须从系统发育的骨干（黑色）部分分支出来。一种可能的修正是假设 2 个采样的骨干谱系是总共 9 个中的 2 个（包括来自其他属的 7 个未采样的分类群），或 0.22。另一种可能性是将主干分数设置为等于整个进化枝的总采样分数。本例中，我们从fu和bar中一共采样了 6 个物种, 共计 15 种。另有7种未从其他属中取样。您可以将骨干分数设置为 6 / (15 + 7) = 0.27。这些都不是完美的解决方案，但我们认为 - 一般来说 - 如果采样分类单元的百分比在您的树中变化很大，则最好使用进化枝特定的采样分数