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

## ----library------------------------------------------------------------------
library(seqRdna)

## ----installation, eval=FALSE-------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE)) {
#     install.packages("BiocManager")
# }
# 
# BiocManager::install("seqRdna")

## ----load-package-------------------------------------------------------------
library(seqRdna)

## ----load-data, message = FALSE, warning = FALSE------------------------------
library(seqRdna)
library(methylKit)

file_path <- system.file(
  "extdata",
  "filtered_meth_raw_chr1_chrRDNAm_small.rds",
  package = "seqRdna",
  mustWork = TRUE
)

loaded_data <- readRDS(file_path)

## ----define-treatment---------------------------------------------------------
treatment <- c(
  0, 0, 0, 0, 0, 0,
  1, 1, 1, 1
)

num_samples <- length(treatment)

## ----define-samples-----------------------------------------------------------
samples <- list(
  'OLD' = c(
    "SRR28431752", "SRR28431753", "SRR28431754", "SRR28431755", 
    "SRR28431756", "SRR28431757"
  ),
  'YOUNG' = c(
    "SRR28431748", "SRR28431749", "SRR28431750", "SRR28431751"
  )
)

## ----create-methylrawlist-----------------------------------------------------
stopifnot(length(loaded_data) == length(treatment))

loaded_data <- methylKit::methylRawList(
  loaded_data,
  treatment = treatment
)

## ----extract-rdna-------------------------------------------------------------
meth_rDNA <- get_rDNA_data(
  loaded_data,
  num_samples,
  rDNA_chr_name = "chrRDNAm"
)

## ----unique-chromosomes-------------------------------------------------------
unique_names <- get_unique_chr_names(
    loaded_data,
    num_samples
)

## ----aged-coverage------------------------------------------------------------
knitr::include_graphics(c(
  "plots/hist_coverage_chr_1.png",
  "plots/hist_coverage_chr_rDNA.png"
))

## ----define-cutoff------------------------------------------------------------
cutoff <- 35000

## ----shift-rdna---------------------------------------------------------------
meth_rDNA_shifted <- shifting(
  meth_rdna_raw = meth_rDNA,
  cutoff = cutoff,
  treatment = treatment
)

## ----absolute-copy-number-----------------------------------------------------
absolute_CN <- absolute_rDNA_CN(
  Meth_object = loaded_data,
  num_samples = num_samples,
  autosomes = TRUE,
  rDNA_chr = "chrRDNAm",
  diploid = TRUE
)

## ----inspect-absolute-copy-number---------------------------------------------
absolute_CN

## ----relative-copy-number-----------------------------------------------------
relative_CN <- relative_rDNA_CN(
  Meth_object = loaded_data,
  num_samples = num_samples,
  autosomes = TRUE,
  rDNA_chr = "chrRDNAm",
  diploid = TRUE
)

## ----inspect-relative-copy-number---------------------------------------------
relative_CN

## ----absolute-v-relative------------------------------------------------------
knitr::include_graphics("plots/relative_vs_absolute_rDNA_copy_number.svg")

## ----sort-rdna-regions--------------------------------------------------------
region_res <- sort_by_rDNA_region(
  meth_rDNA,
  num_samples,
  species = "mm39",
  min_cov = 10
)

## ----aggregate-methylation----------------------------------------------------
region <- "18S"

aggregate_18S <- calculate_aggregate_meth(
  region_res,
  region,
  conditions = samples,
  med = FALSE
)

## ----copy-number-methylation--------------------------------------------------
knitr::include_graphics(
  "plots/absolute_vs_aggregate_rDNA_methylation.svg"
)

## ----volcano-methylation------------------------------------------------------
knitr::include_graphics(
  "plots/volcano_plot_exercised_bin_regress.svg"
)

## ----differential-methylation-sites-------------------------------------------
site_results <- diff_meth(
  meth_data = meth_rDNA_shifted,
  samples = samples,
  species = "mm39",
  num_samples = num_samples,
  min_cov = 10,
  region = "18S",
  model_type = "binomial",
  p_adjust_method = "BH",
  site = TRUE
)

## ----differential-methylation-region------------------------------------------
region_results <- diff_meth(
  meth_data = meth_rDNA_shifted,
  samples = samples,
  species = "mm39",
  num_samples = num_samples,
  min_cov = 10,
  region = "18S",
  model_type = "binomial",
  p_adjust_method = "BH",
  site = FALSE
)

## ----beta-binomial-diagnostics, eval=FALSE------------------------------------
# # Inspect fitted model
# summary(fit)
# 
# # Check convergence code
# fit$fit$convergence
# 
# # Check whether the Hessian is positive definite
# fit$sdr$pdHess
# 
# # Run additional diagnostics
# glmmTMB::diagnose(fit)

## ----session-info-------------------------------------------------------------
sessionInfo()