diff --git a/DESCRIPTION b/DESCRIPTION index d4d1507..784b3b0 100644 --- a/DESCRIPTION +++ b/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: diff --git a/NEWS.md b/NEWS.md index 31b1500..efa4fb8 100644 --- a/NEWS.md +++ b/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 diff --git a/R/findCisSNPs.R b/R/findCisSNPs.R index 36248b5..64fb88f 100644 --- a/R/findCisSNPs.R +++ b/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 #### diff --git a/R/findVML.R b/R/findVML.R index 06d1a32..86454fa 100644 --- a/R/findVML.R +++ b/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, diff --git a/R/lmGE.R b/R/lmGE.R index 33638ed..f823459 100644 --- a/R/lmGE.R +++ b/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), diff --git a/R/selectVariables.R b/R/selectVariables.R index 2bfcaf1..1180ebe 100644 --- a/R/selectVariables.R +++ b/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) } diff --git a/R/summarizeVML.R b/R/summarizeVML.R index 3196b2d..9a5135d 100644 --- a/R/summarizeVML.R +++ b/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% { diff --git a/R/test_covariates.R b/R/test_covariates.R index ea556f9..c964d6d 100644 --- a/R/test_covariates.R +++ b/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" diff --git a/R/test_environmental_matrix.R b/R/test_environmental_matrix.R index 72b43e3..d123abf 100644 --- a/R/test_environmental_matrix.R +++ b/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" diff --git a/R/test_genotype_information.R b/R/test_genotype_information.R index cf4094a..f34266a 100644 --- a/R/test_genotype_information.R +++ b/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" diff --git a/R/test_genotype_matrix.R b/R/test_genotype_matrix.R index 0214078..2732308 100644 --- a/R/test_genotype_matrix.R +++ b/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" diff --git a/R/test_methylation_data.R b/R/test_methylation_data.R index 8a74a6e..8fd365c 100644 --- a/R/test_methylation_data.R +++ b/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" diff --git a/R/ultrastable_cpgs.R b/R/ultrastable_cpgs.R index b33fea9..b9a1a6c 100644 --- a/R/ultrastable_cpgs.R +++ b/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" diff --git a/R/utils.R b/R/utils.R index 7a355ad..cd960f4 100644 --- a/R/utils.R +++ b/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)) } diff --git a/man/test_covariates.Rd b/man/test_covariates.Rd index 48df777..11feb20 100644 --- a/man/test_covariates.Rd +++ b/man/test_covariates.Rd @@ -21,4 +21,7 @@ test_covariates \description{ Concomitant variable drawn from a normal distribution } +\examples{ +summary(test_covariates) +} \keyword{datasets} diff --git a/man/test_environmental_matrix.Rd b/man/test_environmental_matrix.Rd index 8cc55c1..9ebd4d6 100644 --- a/man/test_environmental_matrix.Rd +++ b/man/test_environmental_matrix.Rd @@ -21,5 +21,9 @@ test_environmental_matrix \description{ One hundred environmental exposure variables drawn from a normal distribution with mean=0 and sd=1 for 30 individuals. +} +\examples{ +summary(test_environmental_matrix) + } \keyword{datasets} diff --git a/man/test_genotype_information.Rd b/man/test_genotype_information.Rd index 0f150f7..67f296f 100644 --- a/man/test_genotype_information.Rd +++ b/man/test_genotype_information.Rd @@ -21,5 +21,9 @@ test_genotype_information \description{ Genotype chromosome 21 position metadata obtained from a private genotyping data set. +} +\examples{ +summary(test_genotype_information) + } \keyword{datasets} diff --git a/man/test_genotype_matrix.Rd b/man/test_genotype_matrix.Rd index 4d4fa8c..65371c1 100644 --- a/man/test_genotype_matrix.Rd +++ b/man/test_genotype_matrix.Rd @@ -23,5 +23,9 @@ Genotype matrix example using a gene-dosage model, which encodes the SNPs ordinally depending on the genotype allele charge, such as 2 (AA), 1 (AB) and 0 (BB). Valus were drawn from a binomial distribution with size 2 and probability 0.5. +} +\examples{ +summary(test_genotype_matrix) + } \keyword{datasets} diff --git a/man/test_methylation_data.Rd b/man/test_methylation_data.Rd index f3b474a..97b729a 100644 --- a/man/test_methylation_data.Rd +++ b/man/test_methylation_data.Rd @@ -22,5 +22,9 @@ test_methylation_data Simulated M values of the 3000 probes selected in \emph{test_array_manifest} for 30 individuals. Values were converted from beta values, which were drawn from a bimodal Beta distribution. +} +\examples{ +summary(test_methylation_data) + } \keyword{datasets} diff --git a/man/ultrastable_cpgs.Rd b/man/ultrastable_cpgs.Rd index cbe53f7..885b563 100644 --- a/man/ultrastable_cpgs.Rd +++ b/man/ultrastable_cpgs.Rd @@ -26,5 +26,9 @@ CpGs consistently methylated or unmethylated in every sample (1,737 samples from 30 publically available studies). These CpGs are used to create a "null DNAme variance" distribution in the RAMEN package, from which a threshold is taken to identify Highly Variable Probes. +} +\examples{ +head(ultrastable_cpgs) + } \keyword{datasets} diff --git a/tests/testthat/test-lmGE.R b/tests/testthat/test-lmGE.R index 2964749..1438940 100644 --- a/tests/testthat/test-lmGE.R +++ b/tests/testthat/test-lmGE.R @@ -28,6 +28,28 @@ test_that("lmGE works with BIC", { expect_equal(nrow(lmge_res), 5) }) +test_that("lmGE works when there are no covariates", { + foreach::registerDoSEQ() + lmge_res <- RAMEN::lmGE( + selected_variables = selected_variables_test[1:5, ], + summarized_methyl_VML = summarized_methyl_VML_test, + genotype_matrix = RAMEN::test_genotype_matrix, + environmental_matrix = RAMEN::test_environmental_matrix, + covariates = NULL, + model_selection = "AIC" + ) + expect_true(is.data.frame(lmge_res)) + expect_equal(nrow(lmge_res), 5) + expect_false(any(is.na(lmge_res$tot_r_squared))) + + # VML4's winning model only involves one variable (E). With no covariates + # to pad out the model, this used to break relaimpo::calc.relimp.lm(), + # which requires a model with 2+ regressors + vml_e <- lmge_res[lmge_res$model_group == "E", ][1,] + expect_false(is.na(vml_e$e_r_squared)) + expect_true(is.na(vml_e$g_r_squared)) +}) + test_that("lmGE throws errors when expected", { expect_error( RAMEN::lmGE( diff --git a/tests/testthat/test-selectVariables.R b/tests/testthat/test-selectVariables.R index e820a8e..4e462ac 100644 --- a/tests/testthat/test-selectVariables.R +++ b/tests/testthat/test-selectVariables.R @@ -172,3 +172,30 @@ test_that("selectVariables throws errors when expected", { fixed = TRUE ) }) + +# Test that a VML with a single candidate SNP does not throw an error +test_that("selectVariables does not error when a VML has only one candidate SNP", { + # glmnet::cv.glmnet() requires a predictor matrix with 2+ columns. VML1 in + # the test fixture has exactly one candidate SNP, so running it without + # covariates (which would otherwise pad genot_VMLi to 2+ columns) used to + # throw "x should be a matrix with 2 or more columns". + VML_one_snp <- VML_cis_snps_test[VML_cis_snps_test$VML_index == "VML1"] + expect_equal(length(VML_one_snp$SNP[[1]]), 1) + + result <- expect_no_error( + suppressWarnings(RAMEN::selectVariables( + VML_wSNPs = VML_one_snp, + genotype_matrix = RAMEN::test_genotype_matrix, + environmental_matrix = RAMEN::test_environmental_matrix, + covariates = NULL, + summarized_methyl_VML = summarized_methyl_VML_test, + seed = 1 + )) + ) + + # The lone candidate SNP has nothing to be selected against, so it should + # be kept as-is + expect_true( + VML_one_snp$SNP[[1]] %in% result$selected_genot[[1]] + ) +}) diff --git a/vignettes/RAMEN.Rmd b/vignettes/RAMEN.Rmd index c1a53d9..107c017 100644 --- a/vignettes/RAMEN.Rmd +++ b/vignettes/RAMEN.Rmd @@ -122,11 +122,50 @@ cohort: ## Data expectations RAMEN expects all data sets (genome, exposome, methylome and covariates) to have -undergone quality control, pre-processing and normalization when required. The -choice of methods to do so is beyond the scope of this vignette, but we provide -some broad recommendations and resources here. - -### Genomic data +undergone quality control, pre-processing and normalization when required. The +choice of methods to do so is beyond the scope of this vignette, but we provide +some broad recommendations and resources here. + +### RAMEN's underlying assumptions + +Beyond general QC, there are assumptions specific to how RAMEN models the data +internally that pre-processing choices should be made compatible with: + + - Variance stability in DNAme data: `RAMEN::findVML()` identifies variably + methylated loci (VML) by comparing each probe's variance (or MAD) against a + percentile threshold derived from a reference distribution of probes expected + to show little biological variability (`RAMEN::ultrastable_cpgs`). This + comparison assumes that variance is on a comparable scale across probes and + samples, so pre-processing steps must not compress or inflate variance + differently for the ultrastable (i.e. noise variance) probes versus the rest + of the probes. Therefore, as noted below, we recommend users to: + - Use M-values over Beta-values, since the former are unbounded and roughly + homoscedastic, while Beta-values are bounded (0-1) and compressed near the + extremes. + - Avoid methods relying solely on quantile normalization, since the + rank-based mapping compresses variance more for probes with more extreme + values. + - Normality: `RAMEN::lmGE()` relies on `stats::lm()` to + compare candidate models via AIC/BIC; this comparison is most valid when + residuals are approximately normally distributed, so strongly skewed or + zero-inflated input variables (genomic, methylomic or exposomic) should be + transformed to reduce skew (e.g. log, Box-Cox or Yeo-Johnson transforms) + beforehand. + - No missing values: RAMEN performs no internal handling of missing + data. Every function checks that its numeric inputs (genotype, DNAme, + environmental and covariate matrices) are complete and finite, and will + throw an error if any `NA`, `NaN` or `Inf` value is found. If you use RAMEN + with data types or platforms that can produce missing values, you must resolve + them through imputation or complete-case filtering. Along those lines, please + choose your imputation strategy carefully, and keep in mind that zero has a + meaning across all data modalities, and should not be used as a placeholder + (e.g. a homozygous reference genotype dosage, a methylation value for a CpG, + or no exposure). + +With these assumptions in mind, the following sections provide general +guidance on pre-processing steps for each data type. + +### Genomic data The QC of genotyping data is a crucial step in any genetic study, as it ensures the accuracy and reliability of the results. We recommend the following steps @@ -228,9 +267,17 @@ data set. However, we recommend to take the following into account: - Missing values are very likely to be present. Consider [imputation methods](https://stefvanbuuren.name/fimd/) or complete case analyses. -RAMEN standardizes all exposome variables, so there is no need to do so before. -We also make use of LASSO during the variable selection stage, so the method can -handle correlation in the data set, which is expected. +There is no need to standardize exposome variables before providing them to +RAMEN, even though they typically come in very different scales. During +variable selection, `RAMEN::selectVariables()` uses `glmnet` for its LASSO-based +selection, which by default (`standardize = TRUE`) z-scores all predictors +before fitting and un-standardizes the resulting coefficients, so scale +differences across variables do not bias which ones get selected. Model +comparison in `RAMEN::lmGE()` does not require standardization either, since +AIC/BIC/R2 are invariant to a linear rescaling of a predictor. + +We also make use of LASSO during the variable selection stage, so the method can +handle correlation in the data set, which is expected. ### Covariates to include @@ -1041,7 +1088,7 @@ the package's function can help users in other DNAme analyses, such as: # Frequently Asked Questions -**How can I save the RAMEN data frames that have columns with lists as observations?** +## How can I save the RAMEN data frames that have columns with lists as observations? Saving the data frames produced by RAMEN might seem difficult because it has lists as observations in several columns, which is not supported by some @@ -1074,6 +1121,61 @@ saveRDS(selected_variables, file = "path/selected_variables.Rds") readRDS(file = "path/selected_variables.Rds") ``` +## How can I set up the parallel backend? + +Most functions in RAMEN are compatible with parallel computing (`findVML()`, +`medCorVMR()`, `summarizeVML()`, `selectVariables()`, `lmGE()`, `nullDistGE()`). +You can run these functions sequentially, which is the default in a fresh R +session. This will trigger the following warning, which is just a message: + +```{r} +#> Warning: executing %dopar% sequentially: no parallel backend registered +``` + +To silence it, you can explicitly set up a sequential evaluation before running +any of the above functions: + +```{r eval=FALSE} +foreach::registerDoSEQ() + +# Run RAMEN functions here +``` + +To run the function in parallel, you just need to specify the parallel backend +of your preference. Just to provide examples, you can use the `doParallel` +package: + +```{r eval=FALSE} +library(parallel) +library(doParallel) +# Find the number of CPU cores to use (leaving one core free is common practice) +num_cores <- detectCores() - 1 +#Create the cluster object +cl <- makePSOCKcluster(num_cores) +# Register the parallel backend +registerDoParallel(cl) + +# Run RAMEN functions here + +#Stop the cluster and return to sequential processing when finished: +stopCluster(cl) +registerDoSEQ() +``` + +You can also use the future framework: + +```{r, eval=FALSE} +library(future) + +# Set the backend to use multiple background R sessions +plan(multisession) + +# Explicitly limit the number of worker cores (optional) +plan(multisession, workers = 4) + +# Run RAMEN functions here +``` + # Session info ```{r sessionInfo}