Ropensci sofware peer review second round of revisions

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: RAMEN
 Title: RAMEN: Regional Association of Methylome variability with the Exposome and geNome
-Version: 2.1.0
+Version: 2.1.1
 Authors@R: c(
     person("Erick I.", "Navarro-Delgado", 
             email = "erick.delgado@ubc.ca", 
@@ -33,9 +33,6 @@ Suggests:
     knitr,
     magick,
     rmarkdown,
-    IlluminaHumanMethylation450kanno.ilmn12.hg19,
-    IlluminaHumanMethylationEPICanno.ilm10b4.hg19,
-    IlluminaHumanMethylationEPICv2anno.20a1.hg38,
     ggplot2,
     tidyr,
     doParallel,
@@ -51,6 +48,9 @@ Imports:
     methods,
     relaimpo,
     S4Vectors,
+    IlluminaHumanMethylation450kanno.ilmn12.hg19,
+    IlluminaHumanMethylationEPICanno.ilm10b4.hg19,
+    IlluminaHumanMethylationEPICv2anno.20a1.hg38,
     stats
 VignetteBuilder: knitr
 Depends: 

NEWS.md

@@ -1,3 +1,26 @@
+# RAMEN 2.1.1
+
+This patch fixes some minor bugs in RAMEN including: 
+
+  - Improved argument check - empty objects now throw errors, and all functions
+  check for matched IDs across all objects. 
+  - Fixed bug in `nullDistGE()` where the seed was fixed to 1 during the 
+  `selectVariables()` step instead of using the seed argument provided. 
+  - Fixed a bug in `selectVariables()` where the function would throw an error 
+  when, in the absence of covariates, a VML had only 1 SNP. In that case, the 
+  variable selection for G is now not conducted and the same SNP is returned.
+  - Fixed a bug in `lmGE()` where the function would throw an error when no 
+  covariates were fed into the function. 
+
+In addition to that, we added a section on the vignette explaining how to set
+up a parallel backend, and expanding on the data assumptions/expectations
+from RAMEN.
+
+Finally, the Illumina annotation packages are now a dependency, rather than a 
+suggestion, as they are needed during the package tests. 
+
+This patch presents no changes on the user's end. 
+
 # RAMEN 2.1.0 - July 14, 2026
 
 This minor version presents a couple of changes that improve the integration of

R/findCisSNPs.R

@@ -78,7 +78,7 @@ findCisSNPs <- function(VML, genotype_information, distance = 1e6) {
 
   # Add a VML index to each region if not already existing
   if (!"VML_index" %in% colnames(S4Vectors::mcols(VML))) {
-    S4Vectors::mcols(VML)$VML_index <- paste0("VML", seq_len(length(VML)))
+    S4Vectors::mcols(VML)$VML_index <- paste0("VML", seq_along(VML))
   }
 
   #### Identify the SNPs that are present in each VML ####

R/findVML.R

@@ -11,7 +11,6 @@
 # @return a vector with the names of the probes that conform one reduced region
 
 # @examples
-# \dontrun{
 # target <- data.frame(row.names = c("a", "b", "c", "d"),
 #                      values = c(1, 1, 1, 1))
 # query <- c(2, 1)
@@ -173,9 +172,8 @@ findVML <- function(methylation_data,
                     var_threshold_percentile = 0.99,
                     max_distance = 1000) {
   #### Binding of variables used within the tidyverse framework ####
-  Methyl450_Loci <- epicv2_probes <- TargetID <- chr <- pos <- strand <- NULL
-  var_score <- probes <- . <- median_correlation <- n_VMPs <- VML_index <- NULL
-  type <- seqnames <- start <- end <- width <- NULL
+  Methyl450_Loci <- TargetID <- chr <- pos <- strand <- NULL
+  var_score <- NULL
   #### Argument checks ####
   # methylation_data
   argument_check(methylation_data, "data.frame")
@@ -469,7 +467,7 @@ findVML <- function(methylation_data,
   VML <- c(VMRs, sVMPs)
   S4Vectors::mcols(VML)$type <- c(rep("VMR", length(VMRs)),
                        rep("sVMP", length(sVMPs)))
-  S4Vectors::mcols(VML)$VML_index <- paste0("VML", seq_len(length(VML)))
+  S4Vectors::mcols(VML)$VML_index <- paste0("VML", seq_along(VML))
 
   return(list(
     var_score_threshold = var_threshold,

R/lmGE.R

@@ -194,10 +194,16 @@ lmGE <- function(selected_variables,
   if (!is.null(covariates)) argument_check(covariates, "matrix")
   argument_char_options(object = model_selection, options = c("AIC", "BIC"))
   # All objects have matching IDs
-  vectors_match(rownames(summarized_methyl_VML), colnames(genotype_matrix))
-  vectors_match(rownames(summarized_methyl_VML), rownames(environmental_matrix))
+  vectors_match(rownames(summarized_methyl_VML), colnames(genotype_matrix),
+               object_1_name = "rownames(summarized_methyl_VML)",
+               object_2_name = "colnames(genotype_matrix)")
+  vectors_match(rownames(summarized_methyl_VML), rownames(environmental_matrix),
+               object_1_name = "rownames(summarized_methyl_VML)",
+               object_2_name = "rownames(environmental_matrix)")
   if (!is.null(covariates)) {
-    vectors_match(rownames(summarized_methyl_VML), rownames(covariates))
+    vectors_match(rownames(summarized_methyl_VML), rownames(covariates),
+                 object_1_name = "rownames(summarized_methyl_VML)",
+                 object_2_name = "rownames(covariates)")
   }
   # Matrices have only finite numeric values
   finite_numeric_check(genotype_matrix)
@@ -237,15 +243,22 @@ lmGE <- function(selected_variables,
       } else {
         genot_i <- NULL
       }
-    if (!is.null(covariates)) {
+    if (is.null(covariates)) {
+      covariates_i <- NULL
+    } else {
       covariates_i <- covariates[rownames(summarized_methyl_VML), , drop = FALSE]
     }
     full_data_vml_i <- cbind(summ_vml_i, env_i, genot_i, covariates_i)
     colnames(full_data_vml_i) <- make.names(colnames(full_data_vml_i))
-    # Set the basal model (only covariates)
-    basal_model_formula <- colnames(covariates) |>
-      make.names() |>
-      paste(collapse = " + ")
+    # Set the basal model (only covariates, or the intercept-only model if
+    # there are no covariates)
+    if (!is.null(covariates)) {
+      basal_model_formula <- colnames(covariates) |>
+        make.names() |>
+        paste(collapse = " + ")
+    } else {
+      basal_model_formula <- "1"
+    }
     #### Fit G Models ####
     ## Fit models involving G if G has selected variables
     if (!VML_i$selected_genot %in% empty_lists) {
@@ -433,12 +446,20 @@ lmGE <- function(selected_variables,
                                               "+",
                                               basal_model_formula)
                               )
-      r_decomp <- relaimpo::calc.relimp.lm(
-        object = winning_lm,
-        rela = FALSE,
-        type = "last"
-      ) # This would be the equivalent to using lmg and setting always = covariates.
-      winning_model_VML_i$g_r_squared <- r_decomp$last[make.names(unlist(winning_model_VML_i$variables))[1]]
+      if (is.null(covariates)) {
+        # With no covariates, this model has a single predictor (the SNP), so
+        # its unique contribution to R^2 is simply the model's own R^2.
+        # relaimpo::calc.relimp.lm() requires 2+ regressors and errors out
+        # otherwise.
+        winning_model_VML_i$g_r_squared <- summary(winning_lm)$r.squared
+      } else {
+        r_decomp <- relaimpo::calc.relimp.lm(
+          object = winning_lm,
+          rela = FALSE,
+          type = "last"
+        ) # This would be the equivalent to using lmg and setting always = covariates.
+        winning_model_VML_i$g_r_squared <- r_decomp$last[make.names(unlist(winning_model_VML_i$variables))[1]]
+      }
       winning_model_VML_i$e_r_squared <- NA_real_
       winning_model_VML_i$gxe_r_squared <- NA_real_
     } else if (winning_model_VML_i$model_group == "E") {
@@ -448,13 +469,21 @@ lmGE <- function(selected_variables,
                                               "+",
                                               basal_model_formula)
                               )
-      r_decomp <- relaimpo::calc.relimp.lm(
-        object = winning_lm,
-        rela = FALSE,
-        type = "last"
-      ) # This would be the equivalent to using lmg and setting always = covariates.
+      if (is.null(covariates)) {
+        # With no covariates, this model has a single predictor (the env.
+        # variable), so its unique contribution to R^2 is simply the model's
+        # own R^2. relaimpo::calc.relimp.lm() requires 2+ regressors and
+        # errors out otherwise.
+        winning_model_VML_i$e_r_squared <- summary(winning_lm)$r.squared
+      } else {
+        r_decomp <- relaimpo::calc.relimp.lm(
+          object = winning_lm,
+          rela = FALSE,
+          type = "last"
+        ) # This would be the equivalent to using lmg and setting always = covariates.
+        winning_model_VML_i$e_r_squared <- r_decomp$last[make.names(unlist(winning_model_VML_i$variables))[1]]
+      }
       winning_model_VML_i$g_r_squared <- NA_real_
-      winning_model_VML_i$e_r_squared <- r_decomp$last[make.names(unlist(winning_model_VML_i$variables))[1]]
       winning_model_VML_i$gxe_r_squared <- NA_real_
     } else if (winning_model_VML_i$model_group == "G+E") {
       winning_lm <- stats::lm(data = as.data.frame(full_data_vml_i),

R/selectVariables.R

@@ -160,10 +160,16 @@ selectVariables <- function(VML_wSNPs,
   argument_check(VML_wSNPs$SNP, "list")
   argument_check(VML_wSNPs$VML_index, "character")
   # Check that IDs match across data sets
-  vectors_match(rownames(summarized_methyl_VML), colnames(genotype_matrix))
-  vectors_match(rownames(summarized_methyl_VML), rownames(environmental_matrix))
+  vectors_match(rownames(summarized_methyl_VML), colnames(genotype_matrix),
+               object_1_name = "rownames(summarized_methyl_VML)",
+               object_2_name = "colnames(genotype_matrix)")
+  vectors_match(rownames(summarized_methyl_VML), rownames(environmental_matrix),
+               object_1_name = "rownames(summarized_methyl_VML)",
+               object_2_name = "rownames(environmental_matrix)")
   if (!is.null(covariates)) {
-    vectors_match(rownames(summarized_methyl_VML), rownames(covariates))
+    vectors_match(rownames(summarized_methyl_VML), rownames(covariates),
+                 object_1_name = "rownames(summarized_methyl_VML)",
+                 object_2_name = "rownames(covariates)")
   }
   ## Check that genotype_matrix, environmental_matrix, and covariates (in case
   ## it is provided) have only numeric values and no NA, NaN, Inf
@@ -216,7 +222,7 @@ selectVariables <- function(VML_wSNPs,
     ### Run variable selection
     #### Genotype only model ####
     # Get coefficients with the optimal lambda found by k-fold cross-validation
-    if (any_snp) { # Catch cases when VML dont have surrounding genotyped SNPs
+    if (any_snp && ncol(genot_VMLi) >= 2) { # glmnet requires 2+ predictor columns
       coef_genot <- stats::coef(
         glmnet::cv.glmnet(
           x = genot_VMLi, # Variables
@@ -237,7 +243,11 @@ selectVariables <- function(VML_wSNPs,
       selected_vars_genot <- names(coef_genot)[-1]
       # Remove covariates from selected variables
       selected_vars_genot <- selected_vars_genot[!selected_vars_genot %in% colnames(covariates)]
-    } else {
+    } else if (any_snp) {
+      # Only one SNP and no covariates: glmnet needs 2+ columns to run, and
+      # with a single candidate there is nothing to select against, so keep it
+      selected_vars_genot <- colnames(genot_VMLi)
+    } else { # Catch cases when VML dont have surrounding genotyped SNPs
       selected_vars_genot <- character(0)
     }
 

R/summarizeVML.R

@@ -48,7 +48,7 @@ summarizeVML <- function(VML,
   argument_check(methylation_data, "data.frame")
   # Add a VML index to each region if not already existing
   if (!"VML_index" %in% colnames(S4Vectors::mcols(VML))) {
-    S4Vectors::mcols(VML)$VML_index <- paste0("VML", seq_len(length(VML)))
+    S4Vectors::mcols(VML)$VML_index <- paste0("VML", seq_along(VML))
   }
   if (!all(unique(unlist(VML$probes)) %in% rownames(methylation_data))) {
     warning(paste("Some probes listed in the VML data frame are not found in",
@@ -61,7 +61,7 @@ summarizeVML <- function(VML,
     stop("Please make sure the 'probes' column in the VML object is a list")
   }
   mcols_vml <- S4Vectors::mcols(VML)
-  VML_1 <- i <- NULL # To avoid R CMD check notes
+  i <- NULL # To avoid R CMD check notes
   #### Summarize VML ####
   summarized_VML <- foreach::foreach(i = mcols_vml$VML_index,
                                      .combine = "cbind") %dopar% {

R/test_covariates.R

@@ -9,4 +9,6 @@
 #'   \item{covar1}{Concomitant variable drawn from a normal distribution with
 #'   mean = 0 and sd = 1}
 #' }
+#' @examples
+#' summary(test_covariates)
 "test_covariates"

R/test_environmental_matrix.R

@@ -10,4 +10,7 @@
 #'   \item{E1:100}{100 example environmental exposures; column names correspond
 #'   to exposure IDs}
 #' }
+#' @examples
+#' summary(test_environmental_matrix)
+#'
 "test_environmental_matrix"

R/test_genotype_information.R

@@ -10,4 +10,7 @@
 #'   \item{POS}{Probe genomic position (h19)}
 #'   \item{ID}{SNP ID}
 #' }
+#' @examples
+#' summary(test_genotype_information)
+#'
 "test_genotype_information"

R/test_genotype_matrix.R

@@ -12,4 +12,7 @@
 #'   \item{ID1:30}{Individual's 1 to 30 genotypes; column names correspond to individual IDs}
 #'   ...
 #' }
+#' @examples
+#' summary(test_genotype_matrix)
+#'
 "test_genotype_matrix"

R/test_methylation_data.R

@@ -11,4 +11,7 @@
 #'   \item{ID1:30}{DNAme profile of individuals 1 to 30; column names correspond to individual IDs}
 #'   ...
 #' }
+#' @examples
+#' summary(test_methylation_data)
+#'
 "test_methylation_data"

R/ultrastable_cpgs.R

@@ -13,5 +13,8 @@
 #' A vector with the name of the 15,224 ultrastable probes identified by Edgar
 #' et al. (2014). The name of the probes are based on the Illumina 450k manifest.
 #'
+#' @examples
+#' head(ultrastable_cpgs)
+#'
 #' @source https://static-content.springer.com/esm/art%3A10.1186%2F1756-8935-7-28/MediaObjects/13072_2014_333_MOESM2_ESM.txt
 "ultrastable_cpgs"

R/utils.R

@@ -16,6 +16,21 @@ argument_check <- function(object, data_type, extra_msg = NULL,
                "class.",
                extra_msg))
   }
+
+  is_empty <- switch(data_type,
+                 list = length(object) == 0,
+                 data.frame = nrow(object) == 0 || ncol(object) == 0,
+                 character = length(object) == 0,
+                 matrix = nrow(object) == 0 || ncol(object) == 0,
+                 numeric = length(object) == 0,
+                 GRanges = length(object) == 0
+                 )
+  if (is_empty) {
+    stop(paste("Please make sure the input",
+               obj_name,
+               "is not empty.",
+               extra_msg))
+  }
 }
 
 argument_char_options <- function(object, options, extra_msg = NULL) {
@@ -47,12 +62,15 @@ columns_exist <- function(data.frame, columns, extra_msg = NULL) {
   }
 }
 
-vectors_match <- function(object_1, object_2, extra_msg = NULL) {
-  if (!all(object_1 %in% object_2)) {
+vectors_match <- function(object_1, object_2,
+                          object_1_name = "object_1",
+                          object_2_name = "object_2",
+                          extra_msg = NULL) {
+  if (!setequal(object_1, object_2)) {
     stop(paste("The objects",
-               deparse(substitute(object_1)),
+               object_1_name,
                "and",
-               deparse(substitute(object_2)),
+               object_2_name,
                "must match.",
                extra_msg))
   }