Wrapper to run IQ-TREE

iq_tree.R

#' Run IQ-TREE
#' 
#' For details, see http://www.iqtree.org/doc/
#'
#' @param alignment DNA alignment to use for phylogenetic analysis. Must be matrix
#' (i.e., aligned sequences) of class DNAbin 
#' @param wd Path to working directory. The alignment and IQ-TREE intermediate files
#' and results will be written here.
#' @param bb Optional; number of ultrafast bootstrap replicates to run.
#' @param nt Optional; number of cores to use. Set to "AUTO" to determine automatically.
#' @param alrt Optional; number of SH-aLRT tests to run.
#' @param m Optional; specify model. If no model is given, ModelTest will be run
#' to identify the best model for the data.
#' @param redo Logical; should the analysis be redone from scratch if output from
#' previous runs is present?
#' @param echo Logical; should STDERR be written to the screen?
#' @param ... Other arguments not used by this function but used by
#' drake for tracking.
#' 
#'
#' @return Phylogenetic tree (list of class "phylo")
#' 
#' @examples
#' data(woodmouse)
#' # Rapid boot-strap tree with 1000 replicates on best-fitting model
#' tree <- iq_tree(woodmouse, tempdir(), bb = 1000, echo = TRUE)
#' plot(tree)
#' # Check the optimum number of cores to use for GTR+I+G model
#' iq_tree(tempdir(), woodmouse, m = "GTR+I+G", nt = "AUTO", echo = TRUE, redo = TRUE)
iq_tree <- function (alignment, wd, 
                     bb = NULL, nt = NULL, alrt = NULL, m = NULL, redo = FALSE, 
                     echo = FALSE, ...) {
  
  assertthat::assert_that(inherits(alignment, "DNAbin"),
                          msg = "alignment must be of class 'DNAbin'")
  assertthat::assert_that(is.matrix(alignment),
                          msg = "alignment must be a matrix (not a list of unaligned sequences)")
  
  assertthat::assert_that(assertthat::is.dir(wd))
  
  assertthat::assert_that(is.logical(echo))
  
  assertthat::assert_that(is.logical(redo))
  
  if(!is.null(bb))
    assertthat::assert_that(assertthat::is.number(bb))
  
  if(!is.null(alrt))
    assertthat::assert_that(assertthat::is.number(alrt))
  
  if(!is.null(nt))
    assertthat::assert_that(assertthat::is.number(nt) | assertthat::is.string(nt))
  
  if(!is.null(m))
    assertthat::assert_that(assertthat::is.string(m))
  
  wd <- fs::path_norm(wd)
  
  # check that iqtree is installed and on the PATH
  tryCatch({
    processx::run("iqtree", "-h", echo = FALSE)
  }, warning = function(w) {
    stop("iqtree not installed and on path")
  }, error = function(e) {
    stop("iqtree not installed and on path")
  }, finally = {
    TRUE
  })
  
  # Write alignment to working directory in phylip format
  aln_path <- fs::path(wd, deparse(substitute(alignment))) %>%
    fs::path_ext_set("phy")
  
  phangorn::write.phyDat(alignment, aln_path, format = "phylip")
  
  # Set up arguments
  iqtree_arguments <- c(
    "-s", fs::path_file(aln_path),
    if(!is.null(bb)) "-bb", 
    bb,
    if(!is.null(alrt)) "-alrt", 
    alrt,
    if(!is.null(nt)) "-nt", 
    nt,
    if(!is.null(m)) "-m", 
    m,
    if(isTRUE(redo)) "-redo"
  )
  
  # Run iqtree command
  processx::run(
    "iqtree", 
    iqtree_arguments, wd = wd, echo = echo)
  
  # Read in resulting consensus tree
  tree_path <- fs::path(wd, deparse(substitute(alignment))) %>%
    fs::path_ext_set("phy.contree")
  
  ape::read.tree(tree_path)
  
}