diff --git a/.Rbuildignore b/.Rbuildignore index 321037d..cf17d11 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -1,3 +1,5 @@ +^Meta$ +^doc$ ^appveyor\.yml$ ^.*\.Rproj$ ^\.Rproj\.user$ diff --git a/.gitignore b/.gitignore index cf1d36a..e7122ba 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,5 @@ +Meta +doc .Rproj.user .Rhistory .RData diff --git a/R/aaa-classes.R b/R/aaa-classes.R index 9e83a3e..4d09f3a 100644 --- a/R/aaa-classes.R +++ b/R/aaa-classes.R @@ -228,8 +228,6 @@ ref_genome$lock() #' \describe{ #' \item{`mu()`}{Calculates the average overall mutation rate at equilibrium.} #' \item{`q()`}{Calculates the mutation rate for each nucleotide.} -#' \item{`to_ptr()`}{Converts information in this object to a C++ pointer. -#' You shouldn't need to use this. Ever.} #' } #' #' @return An object of class \code{mevo}. @@ -330,30 +328,6 @@ mevo <- R6::R6Class( # Mutation rates by nucleotides: q <- rowSums(self$Q) + indel return(q) - }, - - - # -------* - # Convert to a XPtr<[Chunk]MutationSampler> object - # -------* - to_ptr = function() { - - stopifnot(is.numeric(self$chunk_size) & !is.na(self$chunk_size)) - - if (self$chunk_size <= 0) { - sampler_ptr <- make_mutation_sampler_base(self$Q, - self$pi_tcag, - self$insertion_rates, - self$deletion_rates) - } else { - sampler_ptr <- make_mutation_sampler_chunk_base(self$Q, - self$pi_tcag, - self$insertion_rates, - self$deletion_rates, - self$chunk_size) - } - - return(sampler_ptr) } ), diff --git a/R/create_variants.R b/R/create_variants.R index 29f822d..48b8f6a 100644 --- a/R/create_variants.R +++ b/R/create_variants.R @@ -1,4 +1,33 @@ +#' Convert to a XPtr<[Chunk]MutationSampler> object +#' +#' @noRd +#' +mevo_obj_to_ptr <- function(mevo_obj) { + + if (!single_integer(mevo_obj$chunk_size, 0)) { + err_msg("create_variants", "mevo_obj", "a \"mevo\" object with a `chunk_size`", + "field that's a single integer >= 0") + } + + if (mevo_obj$chunk_size <= 0) { + sampler_ptr <- make_mutation_sampler_base(mevo_obj$Q, + mevo_obj$pi_tcag, + mevo_obj$insertion_rates, + mevo_obj$deletion_rates) + } else { + sampler_ptr <- make_mutation_sampler_chunk_base(mevo_obj$Q, + mevo_obj$pi_tcag, + mevo_obj$insertion_rates, + mevo_obj$deletion_rates, + mevo_obj$chunk_size) + } + + return(sampler_ptr) +} + + + # doc start ---- #' Create variants from a reference genome. #' @@ -70,7 +99,9 @@ #' @param method_info Object containing information used for the given method. #' See "Method arguments" section for which arguments are used for each method. #' @param mevo_obj A `mevo` object that stores molecular-evolution information. +#' This argument is needed for all methods except `"vcf"`. #' See \code{\link{make_mevo}} for more information. +#' Defaults to `NULL`. #' @param n_cores Number of cores to use for parallel processing. #' This argument is ignored if OpenMP is not enabled. #' Cores are spread across sequences, so it @@ -93,7 +124,7 @@ create_variants <- function(reference, method, method_info, - mevo_obj, + mevo_obj = NULL, n_cores = 1, show_progress = FALSE) { @@ -107,21 +138,17 @@ create_variants <- function(reference, # ---------* if (!inherits(reference, "ref_genome")) { - stop("\nCreating variants can only be done to a ref_genome object.", - call. = FALSE) + err_msg("create_variants", "reference", "a \"ref_genome\" object") } ref_genome_ptr <- reference$genome if (!inherits(ref_genome_ptr, "externalptr")) { - stop("\nYou're attempting to create variants using a \"ref_genome\" object with ", - "a `genome` field that is not of class \"externalptr\". ", - "Restart by reading a FASTA file or by simulating a genome. ", - "And do NOT change the `genome` field manually.", - call. = FALSE) + err_msg("create_variants", "mevo_obj", "a \"mevo\" object with a `genome`", + "field of class \"externalptr\".", + "Restart by reading a FASTA file or by simulating a genome,", + "and do NOT change the `genome` field manually") } if (!single_integer(n_cores, .min = 1)) { - stop("\nThe `n_cores` argument supplied to `create_variants` is not a single", - "whole number greater than or equal to 1.", - call. = FALSE) + err_msg("create_variants", "n_cores", "a single integer >= 1") } @@ -134,18 +161,11 @@ create_variants <- function(reference, # -------+ # Check mevo_obj argument # -------+ - mevo_obj_err <- FALSE - if (missing(mevo_obj)) { - mevo_obj_err <- TRUE - } else if (!inherits(mevo_obj, "mevo")) { - mevo_obj_err <- TRUE - } - if (mevo_obj_err) { - stop("\nIf you want to use a method other than \"vcf\" in ", - "`create_variants`, you must provide the `mevo_obj` argument that is ", - "of class \"mevo\". ", - "You should use the `make_mevo` function to create this object.", - call. = FALSE) + if (is.null(mevo_obj) || !inherits(mevo_obj, "mevo")) { + err_msg("create_variants", "mevo_obj", "a \"mevo\" object (if you", + "want to use a method other than \"vcf\").", + "You should use the `make_mevo` function to create this object") + } # -------+ @@ -160,7 +180,7 @@ create_variants <- function(reference, # -------+ # Make sampler_base_ptr # -------+ - sampler_base_ptr <- mevo_obj$to_ptr() + sampler_base_ptr <- mevo_obj_to_ptr(mevo_obj) # -------+ # Make Gamma matrices (for mutation-rate variability among sites): diff --git a/R/mevo.R b/R/mevo.R index 316756a..9030dfd 100644 --- a/R/mevo.R +++ b/R/mevo.R @@ -1,3 +1,47 @@ +#' Check validity of input parameters for substitions and indels. +#' +#' +#' @noRd +#' +validate_sub_pars <- function(par_list) { + + # Defaults err_msg fxn doesn't quite work + err_msg_ <- "\nThe parameter `%s` used for substititions should be %s." + + for (x in c("lambda", "alpha", "beta", "alpha_1", "alpha_2", "kappa")) { + z <- par_list[[x]] + if (!is.null(z) && !single_number(z, 0)) { + stop(sprintf(err_msg_, x, "a single number >= 0"), call. = FALSE) + } + } + if (!is.null(par_list[["pi_tcag"]]) && + (!is_type(par_list[["pi_tcag"]], "numeric", 4) || + any(par_list[["pi_tcag"]] < 0) || + all(par_list[["pi_tcag"]] == 0))) { + stop(sprintf(err_msg_, "pi_tcag", + paste("numeric vector of length 4, where", + "no values should be < 0 and", + "at least one value should be > 0")), call. = FALSE) + } + if (!is.null(par_list[["abcdef"]]) && + (!is_type(par_list[["abcdef"]], "numeric", 6) || + any(par_list[["abcdef"]] < 0) || + all(par_list[["abcdef"]] == 0))) { + stop(sprintf(err_msg_, "abcdef", + paste("numeric vector of length 6, where", + "no values should be < 0 and", + "at least one value should be > 0")), call. = FALSE) + } + if (!is.null(par_list[["Q"]]) && (!is_type(par_list[["Q"]], "matrix") || + !identical(dim(par_list[["Q"]]), c(4L, 4L)) || + any(par_list[["Q"]] < 0))) { + stop(sprintf(err_msg_, "Q", "a 4x4 matrix, where no values are < 0"), + call. = FALSE) + } + + invisible(NULL) +} + #' Make substitution rate matrix (Q) and vector of equilibrium frequencies (pis). #' @@ -8,57 +52,61 @@ #' substitutions <- function(sub) { - err_msg <- paste("\nNot all required parameters provided for the given", + err_msg_ <- paste("\nNot all required parameters provided for the given", "substitution model (\"%s\").", "See `vignette(\"sub-model\")` for each model's required", "parameter(s).") if (sub$model == "TN93") { if (any(! c("pi_tcag", "alpha_1", "alpha_2", "beta") %in% names(sub))) { - stop(sprintf(err_msg, sub$model), call. = FALSE) + stop(sprintf(err_msg_, sub$model), call. = FALSE) + } + if (any(sub$pi_tcag < 0) || sum(sub$pi_tcag <= 0)) { + stop("\nNo values of `pi_tcag` should be < 0, and the vector needs at ", + "least one value > 0.", call. = FALSE) } Q <- TN93_rate_matrix(sub$pi_tcag, sub$alpha_1, sub$alpha_2, sub$beta) pi_tcag <- sub$pi_tcag } else if (sub$model == "JC69") { if (any(! c("lambda") %in% names(sub))) { - stop(sprintf(err_msg, sub$model), call. = FALSE) + stop(sprintf(err_msg_, sub$model), call. = FALSE) } Q <- JC69_rate_matrix(sub$lambda) pi_tcag <- rep(0.25, 4) } else if (sub$model == "K80") { if (any(! c("alpha", "beta") %in% names(sub))) { - stop(sprintf(err_msg, sub$model), call. = FALSE) + stop(sprintf(err_msg_, sub$model), call. = FALSE) } Q <- K80_rate_matrix(sub$alpha, sub$beta); pi_tcag <- rep(0.25, 4) } else if (sub$model == "F81") { if (any(! c("pi_tcag") %in% names(sub))) { - stop(sprintf(err_msg, sub$model), call. = FALSE) + stop(sprintf(err_msg_, sub$model), call. = FALSE) } Q <- F81_rate_matrix(sub$pi_tcag); pi_tcag <- sub$pi_tcag } else if (sub$model == "HKY85") { if (any(! c("pi_tcag", "alpha", "beta") %in% names(sub))) { - stop(sprintf(err_msg, sub$model), call. = FALSE) + stop(sprintf(err_msg_, sub$model), call. = FALSE) } Q <- HKY85_rate_matrix(sub$pi_tcag, sub$alpha, sub$beta) pi_tcag <- sub$pi_tcag } else if (sub$model == "F84") { if (any(! c("pi_tcag", "beta", "kappa") %in% names(sub))) { - stop(sprintf(err_msg, sub$model), call. = FALSE) + stop(sprintf(err_msg_, sub$model), call. = FALSE) } Q <- F84_rate_matrix(sub$pi_tcag, sub$beta, sub$kappa) pi_tcag <- sub$pi_tcag } else if (sub$model == "GTR") { if (any(! c("pi_tcag", "abcdef") %in% names(sub))) { - stop(sprintf(err_msg, sub$model), call. = FALSE) + stop(sprintf(err_msg_, sub$model), call. = FALSE) } Q <- GTR_rate_matrix(sub$pi_tcag, sub$abcdef) pi_tcag <- sub$pi_tcag } else if (sub$model == "UNREST") { if (any(! c("Q") %in% names(sub))) { - stop(sprintf(err_msg, sub$model), call. = FALSE) + stop(sprintf(err_msg_, sub$model), call. = FALSE) } q_pi_list <- UNREST_rate_matrix(sub$Q) Q <- q_pi_list$Q @@ -85,36 +133,47 @@ indels <- function(indel) { if (is.null(indel)) { rates <- numeric(0) } else { - err_msg <- paste("\nWhen specifying", which_type, "in `make_mevo`, ") - if (is.null(indel$rate)) { - stop(err_msg, "you must always provide a rate.", - call. = TRUE) - } else if (!single_number(indel$rate, 0)) { - stop(err_msg, "the rate must be a single number >= 0.", - call. = TRUE) + err_msg_ <- paste0("\nWhen specifying `", which_type, "` in `make_mevo`, ") + if (!single_number(indel$rate, 0)) { + err_msg("make_mevo", which_type, "a list with a \"rate\" field,", + "and that field must be a single number >= 0") } if (indel$rate == 0) return(numeric(0)) names_ <- sort(names(indel)) if (identical(names_, sort(c("rate", "max_length")))) { if (!single_integer(indel$max_length, 1)) { - stop(err_msg, "the max length must be a single whole number >= 1.", - call. = TRUE) + stop(err_msg_, "the `max_length` field within, if provided, must ", + "be a single whole number >= 1.", call. = FALSE) } rel_rates <- exp(-1 * 1:(indel$max_length)) } else if (identical(names_, sort(c("rate", "max_length", "a")))) { if (!single_integer(indel$max_length, 1)) { - stop(err_msg, "the max length must be a single whole number >= 1.", - call. = TRUE) + stop(err_msg_, "the `max_length` field within, if provided, must ", + "be a single whole number >= 1.", call. = FALSE) + } + if (!single_number(indel$a, 0)) { + stop(err_msg_, "the `a` field within, if provided, must ", + "be a single number >= 0.", call. = FALSE) } L <- 1:(indel$max_length) rel_rates <- {(L * indel$max_length) / (indel$max_length - L + 1)}^(-indel$a) } else if (identical(names_, sort(c("rate", "rel_rates")))) { + if (!is_type(indel$rel_rates, "numeric") || + any(indel$rel_rates < 0) || + all(indel$rel_rates == 0)) { + stop(sprintf(err_msg_, "pi_tcag", + paste("numeric vector of length 4, where", + "no values should be < 0 and", + "at least one value should be > 0")), call. = FALSE) + stop(err_msg_, "the `rel_rates` field within, if provided, must ", + "be a numeric vector, where no values should be < 0 and ", + "at least one value should be > 0.", call. = FALSE) + } rel_rates <- indel$rel_rates } else { - stop(err_msg, "it must contain names that coincide with one of the methods ", - "in the \"Indels\" section within `?make_mevo`.", - "Note that extra names return an error.", - call. = FALSE) + err_msg("make_mevo", which_type, "a list with names that coincide", + "with one of the methods in the \"Indels\" section within", + "`?make_mevo`. Note that extra names return an error.") } # So relative rates sum to 1: @@ -158,14 +217,14 @@ site_variability <- function(site_var, seq_sizes) { shape = site_var$shape) } else if ("mats" %in% names(site_var)) { - err_msg <- paste("\nThe `mats` field inside the `site_var`", + err_msg_ <- paste("\nThe `mats` field inside the `site_var`", "argument to the `make_mevo` function needs to", "be a list of matrices.") if (!inherits(site_var$mats, "list")) { - stop(err_msg, call. = FALSE) + stop(err_msg_, call. = FALSE) } else if (!all(sapply(site_var$mats, inherits, what = "matrix"))) { - stop(err_msg, call. = FALSE) + stop(err_msg_, call. = FALSE) } # Check matrices for proper end points and # columns: @@ -331,6 +390,7 @@ make_mevo <- function(reference, } sub$model <- match.arg(sub$model, c("JC69", "K80", "F81", "HKY85", "TN93", "F84", "GTR", "UNREST")) + validate_sub_pars(sub) sub_info <- substitutions(sub)