R/ReadMrBayes.R

In TreeTools: Create, Modify and Analyse Phylogenetic Trees

#' @importFrom ape read.nexus
.ReadMrBayesTrees <- function(treeFile, burninFrac) {
  trees <- read.nexus(treeFile)
  trees[-seq_len(burninFrac * length(trees))]
}

#' Read posterior tree sample produced by MrBayes
#'
#' Read posterior trees from ['MrBayes'](
#' https://nbisweden.github.io/MrBayes/) output files, discarding burn-in
#' generations.
#'
#' `ReadMrBayesTrees()` samples trees from the posterior distributions
#' computed using the Bayesian inference software ['MrBayes'](
#' https://nbisweden.github.io/MrBayes/) 
#' 
#' @param filepath character string specifying path to `.nex` input file used
#' to initialize the MrBayes analysis, 
#' relative to the R working directory (visible with `getwd()`).
#' @param n Integer specifying number of trees to sample from posterior.
#' @param burninFrac Fraction of trees to discard from each run as burn-in.
#' If `NULL` (the default), this will be read from the last `mcmc` or `mcmcp`
#' command in `filepath`.
#'
#' @return `ReadMrBayesTrees()` returns a 'multiPhylo' object containing 
#' `n` trees sampled evenly from all runs generated by analysis of `filepath`,
#' or `NULL` if no trees are found.
#'
#' @examples
#' \dontrun{ # Download will take a few seconds
#'   url <- 
#'   "https://raw.githubusercontent.com/ms609/hyoliths/master/MrBayes/hyo.nex"
#'   trees <- ReadMrBayesTrees(url, n = 40)
#'   plot(Consensus(trees, p = 0.5))
#' }
#' @family tree import functions
#' @importFrom RCurl url.exists
#' @template MRS
#' @export
ReadMrBayesTrees <- function(filepath, n = NULL, burninFrac = NULL) {
  if (is.null(burninFrac)) {
    lines <- readLines(filepath, warn = FALSE)
    burninPattern <- ".*\\bmcmcp?\\b.*burninfr?a?c?\\s*=\\s*([\\d\\.]+)\\b.*"
    burninFrac <- rev(as.numeric(
      gsub(burninPattern, "\\1", lines[grep(burninPattern, lines, TRUE, TRUE)],
           perl = TRUE, ignore.case = TRUE)
    ))[1]
  }
  runFmt <- paste0(filepath, ".run%s.t")
  run1 <- sprintf(runFmt, 1)
  treeFiles <- if (file.exists(run1)) {
    list.files(dirname(filepath),
               paste0(basename(filepath), "\\.run\\d+\\.t"),
               full.names = TRUE)
  } else if (url.exists(run1)) {
    i <- 2
    repeat{
      if (!url.exists(sprintf(runFmt, i))) break
      i <- i + 1
    }
    sprintf(runFmt, seq_len(i - 1))
  } else {
    stop("Cannot find ", run1)
  }
  trees <- tryCatch(do.call(c, lapply(treeFiles, .ReadMrBayesTrees,
                                      burninFrac = burninFrac)),
                    error = function(e) {
                      warning(e, "Could not read trees from ", filepath)
                      NULL
                    })
  if (is.null(trees)) {
    return(NULL)
  }
  if (!is.null(n)) {
    trees <- trees[seq(1, length(trees), length.out = n)]
  }
}

#' @rdname ReadMrBayesTrees
#' @export
ReadMrBayes <- ReadMrBayesTrees

#' @rdname ReadMrBayesTrees
#' @export
MrBayesTrees <- ReadMrBayesTrees