## ----knitr-opts, include = FALSE----------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----modelstable, echo=FALSE--------------------------------------------------
library(knitr)

models <- data.frame(
  Model = c("DNAme ~ G + covars", "DNAme ~ E + covars", "DNAme ~ G + E + covars", "DNAme ~ G + E + G*E + covars"),
  Name = c("Genetics", "Environmental exposure", "Additive", "Interaction"),
  Abbreviation = c("G", "E", "G+E", "GxE")
)

kable(models, caption = "Fitted models")

## ----pipeline-fig, echo= FALSE, fig.cap="RAMEN pipeline"----------------------
include_graphics("RAMEN_pipeline.png")

## ----install, eval=FALSE------------------------------------------------------
# ## Install RAMEN
# if (!require("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# 
# BiocManager::install("RAMEN")

## ----setup, warning=FALSE, message=FALSE--------------------------------------
# Load the packages used throughout the vignette
library(RAMEN)
library(dplyr)
library(ggplot2)
library(tidyr)
library(doParallel)

# Set the parallel backend to use 2 workers
registerDoParallel(2)

## ----var-mad-demo-------------------------------------------------------------
set.seed(1)
sample <- c(
  rep(0.2, 110),
  sample(x = 0:10, size = 90, replace = TRUE) / 10
)
var(sample)
mad(sample)

## ----mad-var-fig, echo= FALSE, fig.cap="MAD vs var relation"------------------
include_graphics("mad_var.png")

## ----hvps-mad-var-fig, echo= FALSE, fig.cap="HVPs with mad vs var"------------
include_graphics("hvps_mad_var.png")

## ----methylation-input--------------------------------------------------------
# See the structure of the input DNA methylation data set
test_methylation_data[1:5, 1:5]

VML <- findVML(
  methylation_data = test_methylation_data,
  array_manifest = "IlluminaHumanMethylationEPICv1",
  cor_threshold = 0,
  var_method = "variance",
  var_distribution = "ultrastable",
  var_threshold_percentile = 0.99,
  max_distance = 1000
)

# Take a look at the resulting object
# check the specific threshold that was used to label HVPs
glimpse(VML$var_score_threshold)
# check the HVPs identified and their variability score
head(VML$highly_variable_probes)
# Take a look at the identified VML GRanges object
head(VML$VML)

## ----dopar-warning-note-------------------------------------------------------
#> Warning: executing %dopar% sequentially: no parallel backend registered

## ----explore-vml--------------------------------------------------------------
VML_gr <- VML$VML

# Example of an epxloration plot
as(VML_gr, "data.frame") |>
  filter(width > 1) |> # Only plot VMRs, since sVMPs all have a lenght of 1
  ggplot(aes(x = width)) +
  geom_histogram(binwidth = 50, fill = "#BAB4D8") +
  theme_classic() +
  ggtitle("VMRs width (bp)")

## ----run-summarizevml---------------------------------------------------------
summarized_methyl_VML <- summarizeVML(
  VML = VML_gr,
  methylation_data = test_methylation_data
)

# Look at the resulting object
summarized_methyl_VML[1:5, 1:5]

## ----genotype-input-----------------------------------------------------------
# Take a look at the structure of the genotype information data set, which
# contains the genomic location of the SNPs in our genotype data set
head(test_genotype_information)

VML_cis_snps <- findCisSNPs(
  VML = VML_gr,
  genotype_information = test_genotype_information,
  distance = 1e+06
)

# Take a look at the result
head(VML_cis_snps)

## ----cissnps, fig.cap="Disribution of SNPs in cis of VML."--------------------
VML_cis_snps |>
  as("data.frame") |>
  mutate(surrounding_SNPs = case_when(
    surrounding_SNPs > 3000 ~ 3000,
    TRUE ~ surrounding_SNPs
  )) |>
  ggplot(aes(x = surrounding_SNPs)) +
  geom_density() +
  facet_grid("type") +
  xlab("Number of cis SNPs") +
  theme_classic()

# Check the average number of cis snps in out VML data set
mean(VML_cis_snps$surrounding_SNPs)

## ----selectvars-input---------------------------------------------------------
# Take a look at the structure of the genotype, environment and covariate data
# sets, which will be used as input for the variable selection stage
test_genotype_matrix[1:5, 1:5]
test_environmental_matrix[1:5, 1:5]
head(test_covariates)

# Now let's run the function. Setting a seed beforehand makes the variable
# selection reproducible
set.seed(1)
selected_variables <- selectVariables(
  VML_wSNPs = VML_cis_snps,
  genotype_matrix = test_genotype_matrix,
  environmental_matrix = test_environmental_matrix,
  covariates = test_covariates,
  summarized_methyl_VML = summarized_methyl_VML
)

## ----explore-selectedvars-----------------------------------------------------
glimpse(selected_variables)

## ----selectedvars, fig.cap="Number of G and E selected variables."------------
selected_variables |>
  left_join(
    VML_cis_snps |>
      as("data.frame") |>
      select(c(VML_index, type)),
    by = "VML_index"
  ) |>
  transmute(
    VML_index = VML_index,
    type = type,
    Genome = lengths(selected_genot),
    Exposome = lengths(selected_env)
  ) |>
  pivot_longer(-c(VML_index, type)) |>
  dplyr::rename(
    group = name,
    variables = value
  ) |>
  ggplot(aes(x = type, y = variables)) +
  geom_violin() +
  geom_boxplot(width = 0.1, outlier.shape = NA) +
  facet_wrap(~group) +
  ggtitle("Selected variables") +
  theme_classic()

## ----run-lmge-----------------------------------------------------------------
lmge_res <- lmGE(
  selected_variables = selected_variables,
  summarized_methyl_VML = summarized_methyl_VML,
  genotype_matrix = test_genotype_matrix,
  environmental_matrix = test_environmental_matrix,
  covariates = test_covariates,
  model_selection = "AIC"
)

# Check the output
glimpse(lmge_res)

## ----vml-size-----------------------------------------------------------------
(vml_size <- length(VML_cis_snps)) # Number of VML in your data set
desired_obs <- 400000 # We want ~300k informative observations, and we are adding
# 100k extra to account for the VML that will be labelled as Basal (B) during
# selectVariables()
(permutations <- ceiling(desired_obs / vml_size))

## ----run-nulldistge-----------------------------------------------------------
# Compute the null distribution. Setting a seed beforehand makes the
# permutations reproducible
set.seed(1)
null_dist <- nullDistGE(
  VML_wSNPs = VML_cis_snps,
  genotype_matrix = test_genotype_matrix,
  environmental_matrix = test_environmental_matrix,
  summarized_methyl_VML = summarized_methyl_VML,
  permutations = 2,
  covariates = test_covariates,
  model_selection = "AIC"
)

# Take a look at the object
head(null_dist)

## ----nulldist-fig, fig.cap = "R2 difference (winner - basal) in a suffled data set."----
# See the distribution of R2_difference across different winning models
null_dist |>
  drop_na() |> # Remove Basal models from the results, where there is no
  # difference between chosen model and basal model
  ggplot(aes(x = R2_difference)) +
  geom_histogram() +
  facet_grid("model_group") +
  xlab("R2 difference") +
  theme_classic()

## ----nulldist-cutoff----------------------------------------------------------
# Get a cutoff of the 95th percentile of the null distribution for single and
# joint models
cutoff_single <- quantile(
  null_dist |>
    filter(model_group %in% c("G", "E")) |>
    pull(R2_difference),
  0.95
)
cutoff_joint <- quantile(
  null_dist |>
    filter(model_group %in% c("G+E", "GxE")) |>
    pull(R2_difference),
  0.95
)

# Get a data frame with the final results results
final_res <- lmge_res |>
  mutate(
    r2_difference_basal = tot_r_squared - basal_rsquared,
    # Label if the best explanatory model passes its corresponding threshold
    pass_cutoff_threshold = case_when(
      model_group %in% c("G", "E") ~ r2_difference_basal > cutoff_single,
      model_group %in% c("G+E", "GxE") ~ r2_difference_basal > cutoff_joint
    ),
    # Label the final model group, replacing bad performing winning models with
    # "B" (basal)
    model_group = case_when(
      pass_cutoff_threshold ~ model_group,
      TRUE ~ "B"
    )
  ) |>
  select(-pass_cutoff_threshold) # Drop temporary column

# Keep only VML that have informative models with out data
filtered_res <- final_res |>
  filter(!model_group == "B") # Filter based on the cutoff threshold

# Check the VML with informative models
glimpse(filtered_res)

## ----finalresults, fig.cap="Variable Methylated Loci best explanatory models"----
# Plot final results
final_res |>
  group_by(model_group) |>
  summarise(count = n()) |>
  ggplot(aes(x = model_group, y = count)) +
  geom_col() +
  xlab("Best explanatory model") +
  ylab("VML") +
  theme_classic()

table(final_res$model_group)

## ----faq-save-csv, eval=FALSE-------------------------------------------------
# # Example for saving and reading selected_variables object
# data.table::fwrite(selected_variables, file = "path/selected_variables.csv")
# 
# # Read the csv file and make lists the elements in the required columns
# selected_variables <- fread("path/selected_variables.csv", data.table = FALSE) |>
#   mutate(
#     selected_genot = str_split(selected_genot, pattern = "\\|"),
#     # fwrite saves lists as strings separated by |, so we need to split them
#     selected_env = str_split(selected_env, pattern = "\\|"),
#     VMR_index = as.character(VMR_index)
#   )

## ----faq-save-rds, eval=FALSE-------------------------------------------------
# # Example for saving the selected_variables object
# saveRDS(selected_variables, file = "path/selected_variables.Rds")
# 
# # Load the object
# readRDS(file = "path/selected_variables.Rds")

## ----faq-dopar-warning--------------------------------------------------------
#> Warning: executing %dopar% sequentially: no parallel backend registered

## ----faq-backend-seq, eval=FALSE----------------------------------------------
# foreach::registerDoSEQ()
# 
# # Run RAMEN functions here

## ----faq-backend-parallel, 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()

## ----faq-backend-future, 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

## ----sessionInfo--------------------------------------------------------------
sessionInfo()