祖先分布区重建：BiogeoBEARS使用示例
收录于文集
共38篇
BiogeoBEARS是Matzke开发的功能强大的祖先分布区重建R包，介绍参考http://phylo.wikidot.com/biogeobears
示例脚本都可在这个R包的exdata文件夹中a_scripts文件夹中获取
个人能力有限，对脚本的理解或有错误，仅供参考
对应文件可在exdata文件夹中搜索，也可参考http://phylo.wikidot.com/biogeobears#links_to_files
安装
代码块
R
自动换行
复制代码
install.packages("rexpokit")
install.packages("cladoRcpp")
library(devtools)
devtools::install_github(repo="nmatzke/BioGeoBEARS")复制成功
加载及数据了解
加载包
代码块
R
自动换行
复制代码
library(optimx)
library(GenSA)
library(FD)
library(snow)
library(parallel)
library(rexpokit)
library(cladoRcpp)
library(BioGeoBEARS)
复制成功
设置工作路径，输入物种系统发生树信息及地理分布信息
树文件需要是带枝长信息的超度量树（时间树），推荐使用BEAST或MCMCtree等分子钟估算得到的树，分析的时候先用R包ape的 keep.label() 保留内类群tip，去除外类群。需要注意的是，BioGeoBEARS的参数默认时间单位是百万年（My），用MCMCtree生成的树不能直接使用，需要先要对枝长的单位进行处理。
也可以不是超度量树，但需要化石标定，不推荐。
地理信息文件内容应符合如下格式
代码块
R
自动换行
复制代码
3 4 (a b c d)
sp1 1000
sp2 1001
sp3 1111
#第一行为tip数 地理区域数 （各个区域的名称）
#后面每一行对应一个tip的分布，1表示有，0无，sp1就是只分布在a，sp2分布在a和d复制成功
代码块
R
自动换行
复制代码
#示例数据放在这个文件夹中，开头提到的脚本也在这里，可以打开路径查看
extdata_dir = np(system.file("extdata", package="BioGeoBEARS"))
extdata_dir
list.files(extdata_dir)
#树文件输入，需要带枝长信息
trfn = np(paste(addslash(extdata_dir), "Psychotria_5.2.newick", sep=""))
moref(trfn)
#绘制树图
pdffn = "tree.pdf"
pdf(file=pdffn, width=9, height=12)
tr = read.tree(trfn)#这一条很重要，后面需要用到tr这个对象
plot(tr)
title("Example Psychotria phylogeny from Ree & Smith (2008)")
axisPhylo()
dev.off()
cmdstr = paste0("open ", pdffn)
system(cmdstr)
#输入地理分布信息
geogfn = np(paste(addslash(extdata_dir), "Psychotria_geog.data", sep=""))
moref(geogfn)
tipranges = getranges_from_LagrangePHYLIP(lgdata_fn=geogfn)#重要
tipranges
max(rowSums(dfnums_to_numeric(tipranges@df)))
max_range_size = 4 #重要
#max_range_size可以设置为区域总数至树上采样物种分布区域最大数之间的任意数字，现在大部分的研究取后者复制成功
当命令有误，会产生一个错误的pdf文件，打开便会显示如下图，这时候需要删除这个pdf文件，或更换另一个名字运行（正确的命令）。先在命令框运行输入一遍 dev.off() 确保pdf文件处于非打开状态再删除。后面有类似的pdf生成命令同理。
区域数多，状态数就多，状态的数量超过500左右就会难以运行。用以下命令检测状态数。
代码块
R
自动换行
复制代码
numstates_from_numareas(numareas=4, maxareas=4, include_null_range=TRUE)
#>[1] 16
numstates_from_numareas(numareas=4, maxareas=4, include_null_range=FALSE)
#>[1] 15
#状态数过多
numstates_from_numareas(numareas=10, maxareas=10, include_null_range=TRUE)
#>[1] 1024复制成功
运行基础模型
DEC
代码块
R
自动换行
复制代码
#新建项目
BioGeoBEARS_run_object = define_BioGeoBEARS_run()
#输入树和分布信息
BioGeoBEARS_run_object$trfn = trfn
BioGeoBEARS_run_object$geogfn = geogfn
#设置分布区域数最大值
BioGeoBEARS_run_object$max_range_size = max_range_size
#设置最小枝长长度，小于这个枝长的末端会被认作直接祖先（没有物种形成事件），默认的适用枝长单位为百万年。
BioGeoBEARS_run_object$min_branchlength = 0.000001
BioGeoBEARS_run_object$include_null_range = TRUE
## set to FALSE for e.g. DEC* model, DEC*+J, etc. DEC*是对DEC的一个修正，能更好的推断物种在局部范围内灭绝
#提速、多核处理
BioGeoBEARS_run_object$on_NaN_error = -1e50
BioGeoBEARS_run_object$speedup = TRUE
BioGeoBEARS_run_object$use_optimx = TRUE ##FALSE用optim()；如果输入GenSA，使用Generalized Simulated Annealing，适用于高次元模型（参数多于5），但有时无法优化简单的模型
BioGeoBEARS_run_object$num_cores_to_use = 1 ##多线程运算，需要library(parallel) 和/或 library(snow)，数据量大的话一般改为2就行了。
BioGeoBEARS_run_object$force_sparse = FALSE #是否使用稀疏矩阵的项，稀疏矩阵对小规模数据的提速不明显，而且会使结果不精确
#加载dispersal multiplier matrix等文本到模型对象中，并会检测的一些错误。
BioGeoBEARS_run_object = readfiles_BioGeoBEARS_run(BioGeoBEARS_run_object)
#获取良好的默认状态设置
BioGeoBEARS_run_object$return_condlikes_table = TRUE
BioGeoBEARS_run_object$calc_TTL_loglike_from_condlikes_table = TRUE
BioGeoBEARS_run_object$calc_ancprobs = TRUE
#检查输入
check_BioGeoBEARS_run(BioGeoBEARS_run_object)
#运行并储存结果
resfn = "Psychotria_DEC_M0_unconstrained_v1.Rdata"
res = bears_optim_run(BioGeoBEARS_run_object)
save(res, file=resfn)
resDEC = res
#需要重新载入结果时
#load(resfn)
#resDEC = res
复制成功
DEC+j
代码块
R
自动换行
复制代码
BioGeoBEARS_run_object = define_BioGeoBEARS_run()
BioGeoBEARS_run_object$trfn = trfn
BioGeoBEARS_run_object$geogfn = geogfn
BioGeoBEARS_run_object$max_range_size = max_range_size
BioGeoBEARS_run_object$min_branchlength = 0.000001
BioGeoBEARS_run_object$include_null_range = TRUE
BioGeoBEARS_run_object$on_NaN_error = -1e50
BioGeoBEARS_run_object$speedup = TRUE
BioGeoBEARS_run_object$use_optimx = TRUE
BioGeoBEARS_run_object$num_cores_to_use = 1
BioGeoBEARS_run_object$force_sparse = FALSE
BioGeoBEARS_run_object = readfiles_BioGeoBEARS_run(BioGeoBEARS_run_object)
BioGeoBEARS_run_object$return_condlikes_table = TRUE
BioGeoBEARS_run_object$calc_TTL_loglike_from_condlikes_table = TRUE
BioGeoBEARS_run_object$calc_ancprobs = TRUE
#设定DEC+j模型
#从参数嵌套模型中获取ML参数值（需要建立在DEC模型运算的结果之上）
dstart = resDEC$outputs@params_table["d","est"]
estart = resDEC$outputs@params_table["e","est"]
jstart = 0.0001
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["d","init"] = dstart
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["d","est"] = dstart
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["e","init"] = estart
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["e","est"] = estart
#加入j作为自由参数
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","type"] = "free"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","init"] = jstart
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","est"] = jstart
check_BioGeoBEARS_run(BioGeoBEARS_run_object)
resfn = "Psychotria_DEC+J_M0_unconstrained_v1.Rdata"
res = bears_optim_run(BioGeoBEARS_run_object)
save(res, file=resfn)
resDECj = res
#load(resfn)
#resDECj = res
复制成功
DIVALIKE
代码块
R
自动换行
复制代码
BioGeoBEARS_run_object = define_BioGeoBEARS_run()
BioGeoBEARS_run_object$trfn = trfn
BioGeoBEARS_run_object$geogfn = geogfn
BioGeoBEARS_run_object$max_range_size = max_range_size
BioGeoBEARS_run_object$min_branchlength = 0.000001
BioGeoBEARS_run_object$include_null_range = TRUE
BioGeoBEARS_run_object$on_NaN_error = -1e50
BioGeoBEARS_run_object$speedup = TRUE
BioGeoBEARS_run_object$use_optimx = TRUE
BioGeoBEARS_run_object$num_cores_to_use = 1
BioGeoBEARS_run_object$force_sparse = FALSE
BioGeoBEARS_run_object = readfiles_BioGeoBEARS_run(BioGeoBEARS_run_object)
BioGeoBEARS_run_object$return_condlikes_table = TRUE
BioGeoBEARS_run_object$calc_TTL_loglike_from_condlikes_table = TRUE
BioGeoBEARS_run_object$calc_ancprobs = TRUE
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["s","type"] = "fixed"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["s","init"] = 0.0
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["s","est"] = 0.0
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["ysv","type"] = "2-j"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["ys","type"] = "ysv*1/2"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["y","type"] = "ysv*1/2"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["v","type"] = "ysv*1/2"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01v","type"] = "fixed"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01v","init"] = 0.5
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01v","est"] = 0.5
# BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","type"] = "free"
# BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","init"] = 0.01
# BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","est"] = 0.01
check_BioGeoBEARS_run(BioGeoBEARS_run_object)
resfn = "Psychotria_DIVALIKE_M0_unconstrained_v1.Rdata"
res = bears_optim_run(BioGeoBEARS_run_object)
save(res, file=resfn)
resDIVALIKE = res
#load(resfn)
#resDIVALIKE = res复制成功
DIVALIKE+j
代码块
R
自动换行
复制代码
BioGeoBEARS_run_object = define_BioGeoBEARS_run()
BioGeoBEARS_run_object$trfn = trfn
BioGeoBEARS_run_object$geogfn = geogfn
BioGeoBEARS_run_object$max_range_size = max_range_size
BioGeoBEARS_run_object$min_branchlength = 0.000001
BioGeoBEARS_run_object$include_null_range = TRUE
BioGeoBEARS_run_object$on_NaN_error = -1e50
BioGeoBEARS_run_object$speedup = TRUE
BioGeoBEARS_run_object$use_optimx = TRUE
BioGeoBEARS_run_object$num_cores_to_use = 1
BioGeoBEARS_run_object$force_sparse = FALSE
BioGeoBEARS_run_object = readfiles_BioGeoBEARS_run(BioGeoBEARS_run_object)
BioGeoBEARS_run_object$return_condlikes_table = TRUE
BioGeoBEARS_run_object$calc_TTL_loglike_from_condlikes_table = TRUE
BioGeoBEARS_run_object$calc_ancprobs = TRUE
dstart = resDIVALIKE$outputs@params_table["d","est"]
estart = resDIVALIKE$outputs@params_table["e","est"]
jstart = 0.0001
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["d","init"] = dstart
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["d","est"] = dstart
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["e","init"] = estart
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["e","est"] = estart
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["s","type"] = "fixed"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["s","init"] = 0.0
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["s","est"] = 0.0
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["ysv","type"] = "2-j"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["ys","type"] = "ysv*1/2"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["y","type"] = "ysv*1/2"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["v","type"] = "ysv*1/2"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01v","type"] = "fixed"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01v","init"] = 0.5
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01v","est"] = 0.5
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","type"] = "free"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","init"] = jstart
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","est"] = jstart
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","min"] = 0.00001
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","max"] = 1.99999
check_BioGeoBEARS_run(BioGeoBEARS_run_object)
resfn = "Psychotria_DIVALIKE+J_M0_unconstrained_v1.Rdata"
res = bears_optim_run(BioGeoBEARS_run_object)
save(res, file=resfn)
#load(resfn)
#resDIVALIKEj = res复制成功
BAYAREALIKE
代码块
R
自动换行
复制代码
BioGeoBEARS_run_object = define_BioGeoBEARS_run()
BioGeoBEARS_run_object$trfn = trfn
BioGeoBEARS_run_object$geogfn = geogfn
BioGeoBEARS_run_object$max_range_size = max_range_size
BioGeoBEARS_run_object$min_branchlength = 0.000001
BioGeoBEARS_run_object$include_null_range = TRUE
BioGeoBEARS_run_object$on_NaN_error = -1e50
BioGeoBEARS_run_object$speedup = TRUE
BioGeoBEARS_run_object$use_optimx = TRUE
BioGeoBEARS_run_object$num_cores_to_use = 1
BioGeoBEARS_run_object$force_sparse = FALSE
BioGeoBEARS_run_object = readfiles_BioGeoBEARS_run(BioGeoBEARS_run_object)
BioGeoBEARS_run_object$return_condlikes_table = TRUE
BioGeoBEARS_run_object$calc_TTL_loglike_from_condlikes_table = TRUE
BioGeoBEARS_run_object$calc_ancprobs = TRUE
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["s","type"] = "fixed"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["s","init"] = 0.0
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["s","est"] = 0.0
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["v","type"] = "fixed"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["v","init"] = 0.0
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["v","est"] = 0.0
# BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","type"] = "free"
# BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","init"] = 0.01
# BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","est"] = 0.01
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["ysv","type"] = "1-j"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["ys","type"] = "ysv*1/1"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["y","type"] = "1-j"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01y","type"] = "fixed"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01y","init"] = 0.9999
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01y","est"] = 0.9999
check_BioGeoBEARS_run(BioGeoBEARS_run_object)
resfn = "Psychotria_BAYAREALIKE_M0_unconstrained_v1.Rdata"
res = bears_optim_run(BioGeoBEARS_run_object)
save(res, file=resfn)
resBAYAREALIKE = res
#load(resfn)
#resBAYAREALIKE = res复制成功
BAYAREALIKE+j
代码块
R
自动换行
复制代码
BioGeoBEARS_run_object = define_BioGeoBEARS_run()
BioGeoBEARS_run_object$trfn = trfn
BioGeoBEARS_run_object$geogfn = geogfn
BioGeoBEARS_run_object$max_range_size = max_range_size
BioGeoBEARS_run_object$min_branchlength = 0.000001
BioGeoBEARS_run_object$include_null_range = TRUE
BioGeoBEARS_run_object$on_NaN_error = -1e50
BioGeoBEARS_run_object$speedup = TRUE
BioGeoBEARS_run_object$use_optimx = “GenSA”
BioGeoBEARS_run_object$num_cores_to_use = 1
BioGeoBEARS_run_object$force_sparse = FALSE
BioGeoBEARS_run_object = readfiles_BioGeoBEARS_run(BioGeoBEARS_run_object)
BioGeoBEARS_run_object$return_condlikes_table = TRUE
BioGeoBEARS_run_object$calc_TTL_loglike_from_condlikes_table = TRUE
BioGeoBEARS_run_object$calc_ancprobs = TRUE
#从参数嵌套模型中获取ML参数值
dstart = resBAYAREALIKE$outputs@params_table["d","est"]
estart = resBAYAREALIKE$outputs@params_table["e","est"]
jstart = 0.0001
#输入d，e参数起始值
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["d","init"] = dstart
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["d","est"] = dstart
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["e","init"] = estart
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["e","est"] = estart
#无子集共域
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["s","type"] = "fixed"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["s","init"] = 0.0
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["s","est"] = 0.0
#无隔离成种
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["v","type"] = "fixed"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["v","init"] = 0.0
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["v","est"] = 0.0
#允许创始事件
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","type"] = "free"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","init"] = jstart
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","est"] = jstart
#BAYAREALIKE+J模型下，j的上限为1
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","max"] = 0.99999
#调整参数关联性
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["ysv","type"] = "1-j"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["ys","type"] = "ysv*1/1"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["y","type"] = "1-j"
#允许广域的同域物种形成，范围成种。子代范围总与祖先一致。
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01y","type"] = "fixed"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01y","init"] = 0.9999
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01y","est"] = 0.9999
#微调参数上下界，防止崩溃
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["d","min"] = 0.0000001
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["d","max"] = 4.9999999
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["e","min"] = 0.0000001
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["e","max"] = 4.9999999
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","min"] = 0.00001
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","max"] = 0.99999
check_BioGeoBEARS_run(BioGeoBEARS_run_object)
resfn = "Psychotria_BAYAREALIKE+J_M0_unconstrained_v1.Rdata"
res = bears_optim_run(BioGeoBEARS_run_object)
save(res, file=resfn)
#load(resfn)
#resBAYAREALIKEj = res复制成功
绘制结果图
代码块
R
自动换行
复制代码
#仅列举了DEC和DEC+j结果的绘制
pdffn = "Psychotria_M0_unconstrained_v1.pdf"
pdf(pdffn, width=6, height=6)
analysis_titletxt ="BioGeoBEARS DEC on Psychotria M0_unconstrained"
results_object = resDEC
scriptdir = np(system.file("extdata/a_scripts", package="BioGeoBEARS")) #根据树调整边幅
res2 = plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="text", label.offset=0.45, tipcex=0.7, statecex=0.7, splitcex=0.6, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges)
plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="pie", label.offset=0.45, tipcex=0.7, statecex=0.7, splitcex=0.6, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges)
analysis_titletxt ="BioGeoBEARS DEC+J on Psychotria M0_unconstrained"
results_object = resDECj
scriptdir = np(system.file("extdata/a_scripts", package="BioGeoBEARS"))
res1 = plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="text", label.offset=0.45, tipcex=0.7, statecex=0.7, splitcex=0.6, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges)
plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="pie", label.offset=0.45, tipcex=0.7, statecex=0.7, splitcex=0.6, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges)
dev.off()
cmdstr = paste("open ", pdffn, sep="")
system(cmdstr)复制成功
test model
代码块
R
自动换行
复制代码
restable = NULL
teststable = NULL
#建空表储存比较结果
LnL_2 = get_LnL_from_BioGeoBEARS_results_object(resDEC)
LnL_1 = get_LnL_from_BioGeoBEARS_results_object(resDECj)
numparams1 = 3
numparams2 = 2
stats = AICstats_2models(LnL_1, LnL_2, numparams1, numparams2)
stats
res2 = extract_params_from_BioGeoBEARS_results_object(results_object=resDEC, returnwhat="table", addl_params=c("j"), paramsstr_digits=4)
res1 = extract_params_from_BioGeoBEARS_results_object(results_object=resDECj, returnwhat="table", addl_params=c("j"), paramsstr_digits=4)
rbind(res2, res1)
tmp_tests = conditional_format_table(stats)
restable = rbind(restable, res2, res1)
teststable = rbind(teststable, tmp_tests)
LnL_2 = get_LnL_from_BioGeoBEARS_results_object(resDIVALIKE)
LnL_1 = get_LnL_from_BioGeoBEARS_results_object(resDIVALIKEj)
numparams1 = 3
numparams2 = 2
stats = AICstats_2models(LnL_1, LnL_2, numparams1, numparams2)
stats
res2 = extract_params_from_BioGeoBEARS_results_object(results_object=resDIVALIKE, returnwhat="table", addl_params=c("j"), paramsstr_digits=4)
res1 = extract_params_from_BioGeoBEARS_results_object(results_object=resDIVALIKEj, returnwhat="table", addl_params=c("j"), paramsstr_digits=4)
rbind(res2, res1)
conditional_format_table(stats)
tmp_tests = conditional_format_table(stats)
restable = rbind(restable, res2, res1)
teststable = rbind(teststable, tmp_tests)
LnL_2 = get_LnL_from_BioGeoBEARS_results_object(resBAYAREALIKE)
LnL_1 = get_LnL_from_BioGeoBEARS_results_object(resBAYAREALIKEj)
numparams1 = 3
numparams2 = 2
stats = AICstats_2models(LnL_1, LnL_2, numparams1, numparams2)
stats
res2 = extract_params_from_BioGeoBEARS_results_object(results_object=resBAYAREALIKE, returnwhat="table", addl_params=c("j"), paramsstr_digits=4)
res1 = extract_params_from_BioGeoBEARS_results_object(results_object=resBAYAREALIKEj, returnwhat="table", addl_params=c("j"), paramsstr_digits=4)
rbind(res2, res1)
conditional_format_table(stats)
tmp_tests = conditional_format_table(stats)
restable = rbind(restable, res2, res1)
teststable = rbind(teststable, tmp_tests)
#合并表格
teststable$alt = c("DEC+J", "DIVALIKE+J", "BAYAREALIKE+J")
teststable$null = c("DEC", "DIVALIKE", "BAYAREALIKE")
row.names(restable) = c("DEC", "DEC+J", "DIVALIKE", "DIVALIKE+J", "BAYAREALIKE", "BAYAREALIKE+J")
restable = put_jcol_after_ecol(restable)
restable
teststable
#储存
save(restable, file="restable_v1.Rdata")
load(file="restable_v1.Rdata")
save(teststable, file="teststable_v1.Rdata")
load(file="teststable_v1.Rdata")
#write.table(restable, file="restable.txt", quote=FALSE, sep="\t")
#write.table(unlist_df(teststable), file="teststable.txt", quote=FALSE, sep="\t")
restable2 = restable
AICtable = calc_AIC_column(LnL_vals=restable$LnL, nparam_vals=restable$numparams)
restable = cbind(restable, AICtable)
restable_AIC_rellike = AkaikeWeights_on_summary_table(restable=restable, colname_to_use="AIC")
restable_AIC_rellike = put_jcol_after_ecol(restable_AIC_rellike)
restable_AIC_rellike
samplesize = length(tr$tip.label)
AICtable = calc_AICc_column(LnL_vals=restable$LnL, nparam_vals=restable$numparams, samplesize=samplesize)
restable2 = cbind(restable2, AICtable)
restable_AICc_rellike = AkaikeWeights_on_summary_table(restable=restable2, colname_to_use="AICc")
restable_AICc_rellike = put_jcol_after_ecol(restable_AICc_rellike)
restable_AICc_rellike
复制成功
时间分层分析（time-stratified）以DEC+j为例
有的时候，在一段时间内，一些区域不会有物种分布，可以使用这一约束条件，示例是夏威夷群岛的研究数据，有的岛屿在某一时间节点后才出现，在此时间节点前不可能有物种分布。同样的约束思路可以套用到其他情况中去，来引入大的地理事件带来的影响。
代码块
R
自动换行
复制代码
BioGeoBEARS_run_object = define_BioGeoBEARS_run()
BioGeoBEARS_run_object$trfn = trfn
BioGeoBEARS_run_object$geogfn = geogfn
BioGeoBEARS_run_object$max_range_size = max_range_size
BioGeoBEARS_run_object$min_branchlength = 0.000001
BioGeoBEARS_run_object$include_null_range = TRUE
#####导入分层约束条件
BioGeoBEARS_run_object$timesfn = "timeperiods.txt"
BioGeoBEARS_run_object$areas_allowed_fn = "areas_allowed.txt"
BioGeoBEARS_run_object$on_NaN_error = -1e50
BioGeoBEARS_run_object$speedup = TRUE
BioGeoBEARS_run_object$use_optimx = TRUE
BioGeoBEARS_run_object$num_cores_to_use = 1
BioGeoBEARS_run_object$force_sparse = FALSE
BioGeoBEARS_run_object = readfiles_BioGeoBEARS_run(BioGeoBEARS_run_object)
#####分层分析多了下面这一条命令将树按照时间段分割
BioGeoBEARS_run_object = section_the_tree(inputs=BioGeoBEARS_run_object, make_master_table=TRUE, plot_pieces=FALSE)
BioGeoBEARS_run_object$return_condlikes_table = TRUE
BioGeoBEARS_run_object$calc_TTL_loglike_from_condlikes_table = TRUE
BioGeoBEARS_run_object$calc_ancprobs = TRUE
dstart = resDEC$outputs@params_table["d","est"]
estart = resDEC$outputs@params_table["e","est"]
jstart = 0.0001
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["d","init"] = dstart
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["d","est"] = dstart
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["e","init"] = estart
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["e","est"] = estart
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","type"] = "free"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","init"] = jstart
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","est"] = jstart
check_BioGeoBEARS_run(BioGeoBEARS_run_object)
resfn = "Psychotria_DEC+J_M0_time_stratified.Rdata"
res = bears_optim_run(BioGeoBEARS_run_object)
save(res, file=resfn)
resDECj_time_str= res
#load(resfn)
#resDECj_time_str= res复制成功
25.12.19 加笔
最近看到一篇使用了时间分层的文章，分享一下，有兴趣的可以了解一下
Dupin, J., Hong-Wa, C., Gaudeul, M., & Besnard, G. (2024). Phylogenetics and biogeography of the olive family (Oleaceae). Annals of Botany, 134(4), 577-592. https://doi.org/10.1093/aob/mcae100
BSM（biogeographical stochastic mapping）
用于估计生物地理事件的类型与数量
代码块
R
自动换行
复制代码
#load（file="Psychotria_DEC+J_M0_unconstrained_v1.Rdata")
#resDECj=res
model_name = "DEC+J"
res = resDECj
###时间分层约束则改成下面的
#model_name = "DEC+J_time_str"
#res=resDECj_time_str
clado_events_tables = NULL
ana_events_tables = NULL
lnum = 0
BSM_inputs_fn = "BSM_inputs_file.Rdata"
stochastic_mapping_inputs_list = get_inputs_for_stochastic_mapping(res=res)
save(stochastic_mapping_inputs_list, file=BSM_inputs_fn)
#load(BSM_inputs_fn)
##这部分时间分层约束的似乎无法正常运行
names(stochastic_mapping_inputs_list)
stochastic_mapping_inputs_list$phy2
stochastic_mapping_inputs_list$COO_weights_columnar
stochastic_mapping_inputs_list$unconstr
set.seed(seed=as.numeric(Sys.time()))
BSM_output = runBSM(res, stochastic_mapping_inputs_list=stochastic_mapping_inputs_list, maxnum_maps_to_try=100, nummaps_goal=50, maxtries_per_branch=40000, save_after_every_try=TRUE, savedir=getwd(), seedval=12345, wait_before_save=0.01)
#BSM_output = runBSM(res, stochastic_mapping_inputs_list=stochastic_mapping_inputs_list, maxnum_maps_to_try=100, nummaps_goal=50, maxtries_per_branch=40000, save_after_every_try=TRUE, savedir=getwd(), seedval=12345, wait_before_save=0.01, master_nodenum_toPrint=0)
##分层分析的脚本范例多了“master_nodenum_toPrint=0”
RES_clado_events_tables = BSM_output$RES_clado_events_tables
RES_ana_events_tables = BSM_output$RES_ana_events_tables
#load(file="RES_clado_events_tables.Rdata")
#load(file="RES_ana_events_tables.Rdata")
#BSM_output = NULL
#BSM_output$RES_clado_events_tables = RES_clado_events_tables
#BSM_output$RES_ana_events_tables = RES_ana_events_tables
#提取输出
clado_events_tables = BSM_output$RES_clado_events_tables
ana_events_tables = BSM_output$RES_ana_events_tables
include_null_range = TRUE
areanames = names(tipranges@df)
areas = areanames
max_range_size = 4
#library（cladoRcpp）
states_list_0based = rcpp_areas_list_to_states_list(areas=areas, maxareas=max_range_size, include_null_range=include_null_range)
colors_list_for_states = get_colors_for_states_list_0based(areanames=areanames, states_list_0based=states_list_0based, max_range_size=max_range_size, plot_null_range=TRUE)
colors_list_for_states[colors_list_for_states == "#FFFFFF"] = "#FFF5EE"
scriptdir = np(system.file("extdata/a_scripts", package="BioGeoBEARS"))
stratified = FALSE
#stratified = TURE
##分层分析注意这里的替换，要不然会报错
clado_events_table = clado_events_tables[[1]]
ana_events_table = ana_events_tables[[1]]
pdffn = paste0(model_name, "_single_stochastic_map_n1.pdf")
pdf(file=pdffn, width=6, height=6)
master_table_cladogenetic_events = clado_events_tables[[1]]
resmod = stochastic_map_states_into_res(res=res, master_table_cladogenetic_events=master_table_cladogenetic_events, stratified=stratified)
plot_BioGeoBEARS_results(results_object=resmod, analysis_titletxt="Stochastic map", addl_params=list("j"), label.offset=0.5, plotwhat="text", cornercoords_loc=scriptdir, root.edge=TRUE, colors_list_for_states=colors_list_for_states, skiptree=FALSE, show.tip.label=TRUE)
paint_stochastic_map_branches(res=resmod, master_table_cladogenetic_events=master_table_cladogenetic_events, colors_list_for_states=colors_list_for_states, lwd=5, lty=par("lty"), root.edge=TRUE, stratified=stratified)
plot_BioGeoBEARS_results(results_object=resmod, analysis_titletxt="Stochastic map", addl_params=list("j"), plotwhat="text", cornercoords_loc=scriptdir, root.edge=TRUE, colors_list_for_states=colors_list_for_states, skiptree=TRUE, show.tip.label=TRUE)
dev.off()
cmdstr = paste("open ", pdffn, sep="")
system(cmdstr)
include_null_range = include_null_range
areanames = areanames
areas = areanames
max_range_size = max_range_size
states_list_0based = rcpp_areas_list_to_states_list(areas=areas, maxareas=max_range_size, include_null_range=include_null_range)
colors_list_for_states = get_colors_for_states_list_0based(areanames=areanames, states_list_0based=states_list_0based, max_range_size=max_range_size, plot_null_range=TRUE)
colors_list_for_states[colors_list_for_states == "#FFFFFF"] = "#FFF5EE"
scriptdir = np(system.file("extdata/a_scripts", package="BioGeoBEARS"))
stratified = stratified
pdffn = paste0(model_name, "_", length(clado_events_tables), "BSMs_v1.pdf")
pdf(file=pdffn, width=6, height=6)
#50个结果图全部绘制
nummaps_goal = 50
for (i in 1:nummaps_goal)
{
clado_events_table = clado_events_tables[[i]]
analysis_titletxt = paste0(model_name, " - Stochastic Map #", i, "/", nummaps_goal)
plot_BSM(results_object=res, clado_events_table=clado_events_table, stratified=stratified, analysis_titletxt=analysis_titletxt, addl_params=list("j"), label.offset=0.5, plotwhat="text", cornercoords_loc=scriptdir, root.edge=TRUE, colors_list_for_states=colors_list_for_states, show.tip.label=TRUE, include_null_range=include_null_range)
}
dev.off()
cmdstr = paste("open ", pdffn, sep="")
system(cmdstr)
#统计结果的输出保存
areanames = names(tipranges@df)
actual_names = areanames
actual_names
dmat_times = get_dmat_times_from_res(res=res, numstates=NULL)
dmat_times
clado_events_tables = BSM_output$RES_clado_events_tables
ana_events_tables = BSM_output$RES_ana_events_tables
BSMs_w_sourceAreas = simulate_source_areas_ana_clado(res, clado_events_tables, ana_events_tables, areanames)
clado_events_tables = BSMs_w_sourceAreas$clado_events_tables
ana_events_tables = BSMs_w_sourceAreas$ana_events_tables
counts_list = count_ana_clado_events(clado_events_tables, ana_events_tables, areanames, actual_names)
summary_counts_BSMs = counts_list$summary_counts_BSMs
print(conditional_format_table(summary_counts_BSMs))
hist_event_counts(counts_list, pdffn=paste0(model_name, "_histograms_of_event_counts.pdf"))
tmpnames = names(counts_list)
cat("\n\nWriting tables* of counts to tab-delimited text files:\n(* = Tables have dimension=2 (rows and columns). Cubes (dimension 3) and lists (dimension 1) will not be printed to text files.) \n\n")
for (i in 1:length(tmpnames))
{
cmdtxt = paste0("item = counts_list$", tmpnames[i])
eval(parse(text=cmdtxt))
if (length(dim(item)) != 2)
{
next()
}
outfn = paste0(tmpnames[i], ".txt")
if (length(item) == 0)
{
cat(outfn, " -- NOT written, *NO* events recorded of this type", sep="")
cat("\n")
} else {
cat(outfn)
cat("\n")
write.table(conditional_format_table(item), file=outfn, quote=FALSE, sep="\t", col.names=TRUE, row.names=TRUE)
} # END if (length(item) == 0)
} # END for (i in 1:length(tmpnames))
cat("...done.\n")
#比对ML与BSM的祖先状态和范围概率的均值（没太懂）
library(MultinomialCI)
check_ML_vs_BSM(res, clado_events_tables, model_name, tr=NULL, plot_each_node=FALSE, linreg_plot=TRUE, MultinomialCI=TRUE)
#####后面的感觉作用不大就直接复制过来了
#######################################################
# Convert BioGeoBEARS BSM output to phytools
#######################################################
library(phytools)
Q = matrix(c(-3,1,1,1,1,-3,1,1,1,1,-3,1,1,1,1,-3),4,4)
rownames(Q) = letters[1:4]
colnames(Q) = letters[1:4]
Q
set.seed(seed=54321)
simdata = sim.history(tree=tr, Q=Q, nsim=1)
tipdata = simdata$states
tipdata
set.seed(seed=54321)
tr2 = make.simmap(tree=tr, x=tipdata, model="ER")
tr2
tr2$mapped.edge
set.seed(seed=54321)
tr2 = make.simmap(tree=tr, x=sort(tipdata), model="ER")
tr2
tr2$mapped.edge
# Same!
# Plot to PDF
pdffn = "phytools_simmap.pdf"
pdf(file=pdffn, width=6, height=6)
plotSimmap(tr2,lwd=3)
dev.off()
cmdstr = paste0("open ", pdffn)
system(cmdstr)
#######################################################
# Get the ranges_list (useful for colors)
#######################################################
returned_mats = get_Qmat_COOmat_from_BioGeoBEARS_run_object(BioGeoBEARS_run_object=resDEC$inputs)
ranges_list = returned_mats$ranges_list
#######################################################
# Convert a time-stratified BSM to a phytools BSM
#######################################################
res = resDEC
clado_events_table = clado_events_tables[[1]]
ana_events_table = ana_events_tables[[1]]
tr_wSimmap = BSM_to_phytools_SM(res=resDEC, clado_events_table=clado_events_tables[[6]], ana_events_table=ana_events_tables[[6]])
summary(tr_wSimmap)
print(tr_wSimmap)
countSimmap(tr_wSimmap)
# Plot to PDF
pdffn = "BSM_in_phytools_simmap_format.pdf"
pdf(file=pdffn, width=6, height=6)
cols = setNames(colors_list_for_states, ranges_list)
plotSimmap(tr_wSimmap, lwd=3, colors=cols)
dev.off()
cmdstr = paste0("open ", pdffn)
system(cmdstr)
#######################################################
# Convert a list of time-stratified BSMs to a phytools list of BSMs
#######################################################
simmaps_list = BSMs_to_phytools_SMs(res=resDEC, clado_events_tables=clado_events_tables, ana_events_tables=ana_events_tables)
summary(simmaps_list)
print(simmaps_list)
countSimmap(simmaps_list) # Gets big fast, for large states lists
# Plot to PDF
pdffn = "50BSMs_in_phytools_simmap_format.pdf"
pdf(file=pdffn, width=6, height=6)
cols = setNames(colors_list_for_states, ranges_list)
plotSimmap(simmaps_list, colors=cols, lwd=3, hold=FALSE, add=FALSE, plot=TRUE)
dev.off()
cmdstr = paste0("open ", pdffn)
system(cmdstr)
复制成功
2022.11.27加笔
如何指定分布状态的颜色
颜色的设定可以参考BSM的示例代码里的绘图部分，设定colors_list_for_states，然后在绘图命令plot_BioGeoBEARS_results里添加就可以了
下面进行一个演示
代码块
R
自动换行
复制代码
analysis_titletxt ="change_color_before"
results_object = resDEC
scriptdir = np(system.file("extdata/a_scripts", package="BioGeoBEARS"))
plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="pie", label.offset=0.45, tipcex=0.7, statecex=0.7, splitcex=0.6, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges)
复制成功
结果是这样的原设定高饱和亮眼图
下面设置一下状态跟颜色的对应
代码块
R
自动换行
复制代码
include_null_range = TRUE
areanames = names(tipranges@df)
areas = areanames
max_range_size = 4
#library（cladoRcpp）
#排列组合列举所有状态
states_list_0based = rcpp_areas_list_to_states_list(areas=areas, maxareas=max_range_size, include_null_range=include_null_range)
#给所有状态先匹配个颜色
colors_list_for_states = get_colors_for_states_list_0based(areanames=areanames, states_list_0based=states_list_0based, max_range_size=max_range_size, plot_null_range=TRUE)
#更换颜色
colors_list_for_states[colors_list_for_states == "#FFFFFF"] = "#FFF5EE"复制成功
下面截图示例一下
先看一下都有哪些状态
然后看一下自动匹配的颜色列表
抽取一个幸运颜色编号替换所有颜色编号
截一部分是个意思
看一下替换后的颜色列表
全部都一样了呢
绘图
代码块
R
自动换行
复制代码
plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="pie", label.offset=0.45, tipcex=0.7, statecex=0.7, splitcex=0.6, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges,colors_list_for_states=colors_list_for_states)
#要记得把设置好的颜色列表加到参数设定括号里，这里将colors_list_for_states=colors_list_for_states放在了最后复制成功
麻了，#806600这个色号好丑
换色很麻烦，主要是不知道色号，得去翻表，上网查，直接导PS里换颜色方便多了……
25.2.27 加笔
增加了一点注释。
下面介绍一下限制区域间扩散的模型，参考：
Vasconcellos, M. M., Colli, G. R., & Cannatella, D. C. (2020). Paleotemperatures and recurrent habitat shifts drive diversification of treefrogs across distinct biodiversity hotspots in sub‐Amazonian South America. Journal of Biogeography, 48(2), 305-320. https://doi.org/10.1111/jbi.13997
先要准备一个限制区域间扩散的扩散系数 dispersal_multipliers.txt
代码块
R
自动换行
复制代码
a b c d
1 1 1 0.1
1 1 1 1
1 1 1 1
0.1 1 1 1
#这个矩阵设置了扩散权重系数，例子里的意思是a和d区域间的直接扩散的概率权重下调到了其他扩散路径的十分之一（设置为0时禁止a和d之间的直接扩散）复制成功
以DEC模型为例
代码块
R
自动换行
复制代码
#运行DEC*模型
BioGeoBEARS_run_object = define_BioGeoBEARS_run()
BioGeoBEARS_run_object$trfn = trfn
BioGeoBEARS_run_object$geogfn = geogfn
BioGeoBEARS_run_object$max_range_size = max_range_size
BioGeoBEARS_run_object$min_branchlength = 0.000001
BioGeoBEARS_run_object$include_null_range = FALSE
#设置为FALSE则为修正后的非空模型“* model”
#约束区域间扩散
BioGeoBEARS_run_object$dispersal_multipliers_fn = "dispersal_multipliers.txt"
BioGeoBEARS_run_object$on_NaN_error = -1e50
BioGeoBEARS_run_object$speedup = TRUE
BioGeoBEARS_run_object$use_optimx = "GenSA"
BioGeoBEARS_run_object$num_cores_to_use = 1
BioGeoBEARS_run_object$force_sparse = FALSE
BioGeoBEARS_run_object = readfiles_BioGeoBEARS_run(BioGeoBEARS_run_object)
BioGeoBEARS_run_object$return_condlikes_table = TRUE
BioGeoBEARS_run_object$calc_TTL_loglike_from_condlikes_table = TRUE
BioGeoBEARS_run_object$calc_ancprobs = TRUE
check_BioGeoBEARS_run(BioGeoBEARS_run_object)
#设置结果文件名称
resfn = "DEC.Rdata"
#运行
res = bears_optim_run(BioGeoBEARS_run_object)
#储存结果
save(res, file=resfn)
#环境中暂存DEC*模型结果
resDEC = res
复制成功
设置有一定的用处，但是有可能会过拟合或与 +j 模型的理念冲突，谨慎使用。
原文的代码中还有去除指定分布区组合的部分，这一功能在RASP里可以方便的设置，个人十分推荐这样做，可以剔除不合理的分布区域组合（如不连续的间断分布形式）。
代码块
R
自动换行
复制代码
########################
# Modify the states_list of allowed ranges manually for 3 max_range_size (edits Mariana Vasconcellos)
BioGeoBEARS_run_object$states_list
areas = names(tipranges@df)
states_list = rcpp_areas_list_to_states_list(areas=areas, maxareas=max_range_size, include_null_range=TRUE)
states_list
#这里会从0开始对应区域状态，对应的时候需要注意
TF = rep(TRUE, 42) ## Total number of ranges for 6 areas and 3 max_range_size = 42
# Removing 10 areas in your model (non-neighboring areas in the ancestral range)
TF[10] = FALSE
TF[12] = FALSE
TF[14] = FALSE
TF[19] = FALSE
TF[21] = FALSE
TF[24] = FALSE
TF[29] = FALSE
TF[31] = FALSE
TF[37] = FALSE
TF[40] = FALSE
states_list = states_list[TF]
#######################
BioGeoBEARS_run_object$states_list = states_list
复制成功
利用BSM数据推断物种形成方式的概率可视化还是挺有意思的，这里顺便搬一下
代码块
R
自动换行
复制代码
library("tidyverse")
library("dplyr")
library("ggplot2")
####BSM运行了100个平行
####下面这一部分中的34和35需要根据自己的情况进行修改
myData <- dplyr::as_tibble(matrix(NA, ncol = 34, nrow = 100)) #34改成树中tip的数量-1
for (i in 1:100) {
for (j in 1:34) { #34改成树中tip的数量-1
colnames(myData)[j] <- paste0("node", j+35)#35改成tip数
myData[i, j] <- clado_events_tables[[i]]$clado_event_type[j+35]#35改成tip数
}
}
myTidyData <- pivot_longer(myData, everything(), names_to = "node", values_to = "rangeshift")
# make a two-way table to use as source ofr piechart data on tree
myTable <- table(myTidyData$node, myTidyData$rangeshift)
myTable <- as.matrix(myTable)
# calculate proportions from raw numbers
myTable <- myTable/100
tr = read.tree("./input_files/hyptree_R.tre")
# truncate tiplabels to fit page
tr$tip.label <- str_sub(tr$tip.label, 1, 20)
# Set colors for piecharts and subtitle of plot
if (model_name == "DECJ") {
plot_subtitle <- "red=jump, white=subset, gray=sympatry, green=vicariance"
myColors <- c("red", "white", "gray", "green")
} else { # DEC
plot_subtitle <- "white=subset, gray=sympatry, green=vicariance"
myColors <- c("white", "gray", "green")#这里需要注意，这个颜色映射向量直接根据列数来对应的，如果所有节点的subset都为0，则在表格中不会存在这一列，这样会导致颜色对应错位，饼图上sympatry就变成了白色
}
# Open a PDF
pdffn = paste0(outputdir, model_name,
"_Cladogenetic_Events.pdf")
pdf(file=pdffn, width=8.5, height=11)
plot(tr, sub = plot_subtitle)
title(paste0(model_name, ": Prob of types of cladogenetic events"))
axisPhylo()
# Add pie diagrams with proportion of type of cladogenetic event
nodelabels(pie = myTable, piecol = myColors, cex = 0.9)
dev.off()
cmdstr = paste("open ", pdffn, sep="")
system(cmdstr)
### END ###复制成功
cv18805296
分享至
投诉或建议
