BAMM的全称是Bayesian Analysis of Macroevolutionary Mixtures，使用的方法就是贝叶斯，运行后首先就要看看运行是否收敛
> library(BAMMtools)
> setBAMMpriors(tree)
> mcmcout <- read.csv("mcmc_out.txt",header = T) #读取运行数据
> plot(mcmcout$logLik~mcmcout$generation)# 。i输出运行数据查看是否收敛
> burnstart <- floor(0.2 * nrow(mcmcout))
> postburn <- mcmcout[burnstart:nrow(mcmcout), ]
> library(coda) #用coda包来检查有效事件数，数值在200以上算是比较可靠
> effectiveSize(postburn$N_shifts)
var1
467.1729
> effectiveSize(postburn$logLik)
var1
263.742
> tree <- read.tree("FigTree_2023070802.convert.topo_gesneri.tree") #读取时间树
> edata <- getEventData(phy = tree,eventdata ='event_data.txt',burnin = 0.1, type = "diversification") #读取时间树和分化速率数据
> plot(edata, lwd=2.5, labels= T, cex = .4) #对读取的数据进行绘图，lwd为枝的宽度，labels为是否显示枝顶的label, cex为label的字的大小
结果如下（画出的树形图，跟原树形图排列好像是相反的）
cc82ff7f613322bac86d1786dc9f863.png
提取在置信区间95%以内的结果集
css <- credibleShiftSet(edata, expectedNumberOfShifts = 1,set.limit = 0.95)
css$number.distinct
# [1] 9
##位于95%置信区间内的预测配置（这里姑且将速率转换位置的系统发育拓扑翻译为‘配置’）结果数，这个数字会随着运行代数的增加而增加，当运行代数从这里的十万提升到一千万时，95%置信区间内的配置数为457，多次平行的运行BAMM可能得到不太相同但相近的结果数。
summary(css)
#通过这个函数可以查看更多关于配置集的
plot(css)
#绘制数据集的所有配置，这里就不展示结果了，但是可以粗略观察到，有两个速率变换位置在大多数配置中是比较固定的。
获取后验概率最高的配置
best <- getBestShiftConfiguration(edata, expectedNumberOfShifts=1)
plot.bammdata(best, lwd = 2)
addBAMMshifts(best, cex=2.5)
结果类似下图 （不知道为啥点的位置跑偏了）
1689929645(1).png
为了更直观的来接不同配置，这里将后验和相关配置绘制出来。
dsc <- distinctShiftConfigurations(edata, expectedNumberOfShifts=1, threshold=5)
plot(dsc,edata,rank = 1,legend = F)
#最佳配置的速率转换概率
plot(dsc,edata,rank = 2,legend = F)
用 plotRateThroughTime 绘制分化速率曲线
plotRateThroughTime(edata,ratetype = "'auto", xlim = c(100,0), ylim=c(0,0.2)) # auto就是speciation ，如果曲线尾部出现异常，可以使用 end.time = 0.1参数去掉最尾部的一截
plotRateThroughTime(edata,ratetype = "netdiv", xlim = c(100,0), ylim=c(0,0.2))
plotRateThroughTime(edata,ratetype = "extinction", xlim = c(100,0), ylim=c(0,0.2)) #xlim设置x轴范围，ylim设置Y轴范围
获取特定分支的速率图
先读取两个能代表该分支的物种，分别为最顶端和最基部的样品
species1 <- "GN007"
species2 <- "GN045"#读取两个样品
tipnode1 <- which(tree$tip.label == species1)
tipnode2 <- which(tree$tip.label == species2)#将两个样品赋予到节点
clade1 <- getMRCA(tree, tip = c(tipnode1, tipnode2))#将节点赋予到分支
plotRateThroughTime(edata, ratetype = "diversification", node = clade1, nodetype="include", xlim = c(100,0), ylim=c(0,0.2)) #nodetyoe="include" 表示只显示该分支的速率，exclude表示显示除了该分支以为其他分支的速率
可以将speciation、netdiv、和extinction 3个曲线合并到一张图中
par(mfrow=c(3,1)) #使用par函数将图形分成3行1列的分布模式，函数中第一个数字表示行数，第二个数字表示列数
plotRateThroughTime(edata,ratetype = "'auto", xlim = c(100,0), ylim=c(0,0.2)) # auto为speciation
plotRateThroughTime(edata,ratetype = "netdiv", xlim = c(100,0), ylim=c(0,0.2))
plotRateThroughTime(edata,ratetype = "extinction", xlim = c(100,0), ylim=c(0,0.2)) #xlim设置x轴范围，ylim设置Y轴范围
结果如下图
1690006637(1).png
网友评论