R/seq-gen_wrapper.R

In wadedismukes/treeduckenTools: treeduckenTools

#'
#' Simulate sequences using Seq-Gen
#'
#' @description
#' This function provides a more friendly interface around the seqgen wrapper
#' used in phyclust. Parameters are provided as numeric values and sent to
#' seqgen and the alignment is returned as a matrix.
#'
#' @importFrom phyclust seqgen

run_seqgen_HKY <- function(tree, equilibrium_probabilities, tstv_ratio, seq_length,
    num_replicates, tree_scale = 1, branch_scale = 1, seed = NULL) {
    equil_prob <- paste(equilibrium_probabilities, collapse = ",")
    opts_list <- c(paste0("-m", "HKY"), paste0("-l", seq_length), paste0("-t",
    tstv_ratio), paste0("-f", equil_prob), paste0("-n", num_replicates))
    seqgen(opts = paste(opts_list, collapse = " "), newick.tree = tree)
}

run_seqgen_JC <- function(tree, tree_scale = 1, branch_scale = 1,
                          seq_length, num_replicates, seed = NULL) {
    opts_list <- c(paste0("-m", "HKY"), paste0("-l", seq_length),
                   paste0("-n", num_replicates))
    seqgen(opts = paste(opts_list, collapse = " "), newick.tree = tree)
}


run_seqgen_GTR <- function(tree, equilibrium_prob, rate_matrix, seq_length,
                               prop_invariant = NULL, gamma_rate_categories = NULL,
                               scale_tree = NULL, scale_branches = NULL,
                               replicates = NULL, rng_seed = NULL) {
    opts <- "-mGTR"
    equil_prob_str <- paste(equilibrium_prob, collapse = ",")
    rate_mat_str <- paste(rate_matrix, collapse = ",")
    opts <- paste(opts, paste0("-l", seq_length))
    opts <- paste(opts, paste0("-f", equil_prob_str))
    opts <- paste(opts, paste0("-r", rate_mat_str))
    if (!is.null(prop_invariant)) {
        opts <- paste(opts, paste0("-i", prop_invariant))
    }
    if (!is.null(gamma_rate_categories)) {
        opts <- paste(opts, paste0("-g", gamma_rate_categories))
    }
    if (!is.null(scale_tree)) {
        opts <- paste(opts, paste0("-d", scale_tree))
    }
    if (!is.null(scale_branches)) {
        opts <- paste(opts, paste0("-s", scale_branches))
    }
    seqgen(opts = opts, newick.tree = tree)
}