Skip to content

Instantly share code, notes, and snippets.

@cmt2
Created April 7, 2022 21:09
Show Gist options
  • Save cmt2/f8f875461ae993e35f51f2b41877913c to your computer and use it in GitHub Desktop.
Save cmt2/f8f875461ae993e35f51f2b41877913c to your computer and use it in GitHub Desktop.
load RevBayes stochastic mapping files into R to plot with phytools, etc.
readStochasticMaps <- function(map_paths,
burnin = 0.1,
combine_runs = TRUE) {
# read in map logs and remove burnin
maps <- vector("list", length = length(map_paths))
for (i in 1:length(map_paths)) {
message(paste0("Reading in map log file ", i, "\n",
"This may take a while."))
maps[[i]] <- read.table(file = map_paths[i],
sep = "\t",
header = T,
row.names = 1)
end_burnin <- round(burnin*nrow(maps[[i]]), digits = 0)
maps[[i]] <- maps[[i]][c(end_burnin : nrow(maps[[i]])),]
}
# combine log files
if (combine_runs == T) {
maps <- list(do.call(rbind, maps))
}
# get maps from output
maps_only <- vector("list", length = length(maps))
for (i in 1:length(maps)) {
maps_only[[i]] <- vector("list", length = nrow(maps[[i]]))
message(paste0("Converting trace ",
i,
" to a list of SIMMAP objects",
"\n",
"This may take a while."))
for (j in 1:nrow(maps[[i]])){
maps_only[[i]][[j]] <-
phytools::read.simmap(text = maps[[i]][j,"simmap"])
}
}
# reverse branch order
# revbayes outputs the maps in reverse order along branches
for (i in 1:length(maps_only)) {
for (j in 1:length(maps_only[[i]])) {
for (k in 1:length(maps_only[[i]][[j]]$maps))
maps_only[[i]][[j]]$maps[[k]] <- rev(maps_only[[i]][[j]]$maps[[k]])
}
}
return(maps_only)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment