Skip to content

Commit

Permalink
Adde more checks to mevo, made mevo_obj_to_ptr sep. fxn [skip ci]
Browse files Browse the repository at this point in the history
  • Loading branch information
@lucasnell
lucasnell committed Mar 20, 2019
1 parent ba9d1dd commit 3bb1ef9
Show file tree
Hide file tree
Showing 5 changed files with 135 additions and 77 deletions.
2 changes: 2 additions & 0 deletions .Rbuildignore
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
^Meta$
^doc$
^appveyor\.yml$
^.*\.Rproj$
^\.Rproj\.user$
Expand Down
2 changes: 2 additions & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
Meta
doc
.Rproj.user
.Rhistory
.RData
Expand Down
26 changes: 0 additions & 26 deletions R/aaa-classes.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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}.
Expand Down Expand Up @@ -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)
}

),
Expand Down
68 changes: 44 additions & 24 deletions R/create_variants.R
View file
Original file line number Diff line number Diff line change
@@ -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.
#'
Expand Down Expand Up @@ -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
Expand All @@ -93,7 +124,7 @@
create_variants <- function(reference,
method,
method_info,
mevo_obj,
mevo_obj = NULL,
n_cores = 1,
show_progress = FALSE) {

Expand All @@ -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")
}


Expand All @@ -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")

}

# -------+
Expand All @@ -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):
Expand Down
114 changes: 87 additions & 27 deletions R/mevo.R
View file
Original file line number Diff line number Diff line change
@@ -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).
#'
Expand All @@ -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
Expand All @@ -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:
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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)


Expand Down

0 comments on commit 3bb1ef9

Please sign in to comment.