library(tidyverse)
library(fs)
library(ape)
library(phangorn)
library(glue)
set.seed(123)
MrBayes postprocessing
After MrBayes
is finished, a little postprocessing is needed. I will do this in R
, because it has useful libraries for this kind of work.
The goal is
- to throw away the first half of each mcmc,
- remove dummy taxa,
- sample 100 trees at random from the remaining posterior distributions for each family
- compute a maximum clade credibility tree for each family.
preamble
loading libraries and setting a random seed for reproducability
workflow
Full disclosure: most of the code was written by ChatGPT, with some natural language input from me.
Here are some helper function, whose purpose should be clear from their names.
# Helper function to extract latter half of trees
<- function(tree_list) {
get_latter_half ::map(tree_list, ~ .x[round(length(.x) / 2):length(.x)])
purrr
}
# Helper function to remove dummy taxa if it exists
<- function(tree) {
remove_dummy if ("dummy" %in% tree$tip.label) {
drop.tip(tree, "dummy")
else {
}
tree
}
}
# Check file completion
<- function(file) {
is_file_complete if (file_exists(file)) {
<- read_lines(file) %>% tail(1)
last_line return(str_detect(last_line, "end;"))
} return(FALSE)
}
Now the posterior tree changes are extracted, cleaned and sampled.
# Get the list of family names
<- dir_ls("mrbayes/output") %>%
families ::map_chr(fs::path_file) %>%
purrrstr_remove("\\..*") %>%
unique() %>%
::keep(str_detect, "^0")
purrr
# Check if all files in the family have completed
<- purrr::map_lgl(families, function(fm) {
families_completed <- glue::glue("mrbayes/output/{fm}.run{1:2}.t")
files all(purrr::map_lgl(files, is_file_complete))
})
names(families_completed) <- families
# Extract and process trees
<- purrr::map(families, function(fm) {
fm_trees print(fm)
<- glue::glue("mrbayes/output/{fm}.run{1:2}.t")
files
if (families_completed[fm]) {
<- purrr::map(files, read.nexus)
trees else {
} <- purrr::map(files, function(file) {
trees_as_string <- readLines(file, warn = FALSE)
tree_string c(tree_string, "end;")
})<- purrr::map(trees_as_string, ~ read.nexus(textConnection(.x)))
trees
}
<- do.call(c, get_latter_half(trees))
latter_half_trees if (length(latter_half_trees) < 100) {
stop(glue("Not enough trees in {fm}. Expected at least 100 but found {length(latter_half_trees)}."))
}::map(sample(latter_half_trees, 100), remove_dummy)
purrr
})
names(fm_trees) <- families
writing the posterior samples to disk…
# Create directories and write trees
= "../data/posterior_trees"
posterior_dir
dir_create(posterior_dir)
::walk(names(fm_trees), function(fm) {
purrr<- glue("{posterior_dir}/{fm}_posterior.tree")
file_name write.tree(fm_trees[[fm]], file = file_name)
})
computing the maximum clade credibility trees and writing them to disk…
= "../data/mcc_trees"
mcc_dir
dir_create(mcc_dir)
::walk(families, function(fm) {
purrr<- glue("{mcc_dir}/{fm}_mcc.tre")
file_name <- maxCladeCred(fm_trees[[fm]])
tr write.tree(tr, file = file_name)
})