## ----vignette_setup, include=FALSE--------------------------------------------
knitr::opts_chunk$set(
    echo = TRUE,
    warning = FALSE,
    message = FALSE,
    fig.width = 12,
    fig.height = 8,
    dpi = 200,
    fig.retina = 1
)

## ----html_styles, echo=FALSE, results='asis'----------------------------------
cat("
<style>
.main-container {
    max-width: 98vw !important;
    width: 98vw !important;
}
.container-fluid {
    max-width: 98vw !important;
}
</style>
")

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

## ----setup--------------------------------------------------------------------
suppressPackageStartupMessages({
    library(CMEnt)
    library(GenomicRanges)
    library(ggplot2)
})

genome <- "hg19"

## ----load_inputs--------------------------------------------------------------
# loadExampleInputData fetches example resources from DMRsegaldata via
# ExperimentHub and assigns them into the current environment.
loadExampleInputDataChr5And11("beta")
loadExampleInputDataChr5And11("pheno")
loadExampleInputDataChr5And11("array_type")
# DMPs are produced from minfi, and have the form of a dataframe
loadExampleInputDataChr5And11("dmps")
cat(sprintf("Beta matrix dimensions: %d sites x %d samples\n", nrow(beta), ncol(beta)))
cat(sprintf("Phenotype data dimensions: %d samples x %d columns\n", nrow(pheno), ncol(pheno)))
cat(sprintf("Array type: %s\n", array_type))
cat(sprintf("Number of DMPs (seeds): %d\n", nrow(dmps)))
cat(sprintf("Head of phenotype data:\n"))
print(head(pheno))
cat(sprintf("Head of DMPs:\n"))
print(head(dmps))

# The beta handler retrieves beta values for DMRs and their seeds.
# Reusing it avoids redundant preprocessing across DMR calling and scoring.
beta_handler <- getBetaHandler(beta, array = array_type, genome = genome)


## ----build_dmrs---------------------------------------------------------------
pheno[, "case_control"] <- pheno$Sample_Group == "cancer"
dmrs <- readRDS(system.file("extdata", "example_outputChr5And11.rds", package = "CMEnt"))
cat(sprintf("DMRs found: %d\n", length(dmrs)))
cat(sprintf("DMRs metadata columns: %s\n", paste(colnames(mcols(dmrs)), collapse = ", ")))
cat(sprintf(
    "DMRs genomic range example: %s\n",
    paste0(seqnames(dmrs)[1], ":", start(dmrs)[1], "-", end(dmrs)[1])
))

## ----dmr_summary--------------------------------------------------------------
dmr_meta <- as.data.frame(mcols(dmrs))
summary_tbl <- data.frame(
    chr = as.character(seqnames(dmrs)),
    start = start(dmrs),
    end = end(dmrs),
    width = width(dmrs),
    seeds_num = mcols(dmrs)$seeds_num,
    sites_num = mcols(dmrs)$sites_num,
    delta_beta = round(mcols(dmrs)$delta_beta, 4),
    pval = if ("pval" %in% colnames(dmr_meta)) signif(dmr_meta$pval, 3) else NA_real_,
    qval = if ("qval" %in% colnames(dmr_meta)) signif(dmr_meta$qval, 3) else NA_real_,
    score = if ("score" %in% colnames(dmr_meta)) round(dmr_meta$score, 4) else NA_real_,
    cv_accuracy = if ("cv_accuracy" %in% colnames(dmr_meta)) round(dmr_meta$cv_accuracy, 4) else NA_real_,
    block_id = if ("block_id" %in% colnames(dmr_meta)) dmr_meta$block_id else NA_character_
)
if (any(is.finite(summary_tbl$qval))) {
    summary_tbl <- summary_tbl[order(summary_tbl$qval, -abs(summary_tbl$delta_beta), -summary_tbl$score, na.last = TRUE), ]
} else {
    summary_tbl <- summary_tbl[order(-abs(summary_tbl$delta_beta), -summary_tbl$score, na.last = TRUE), ]
}
knitr::kable(
    head(summary_tbl, 5),
    caption = paste(
        "Top 5 DMRs by statistical evidence and effect size,",
        "with SVM score as a complementary measure"
    )
)

## ----motif_discovery----------------------------------------------------------
n_motif_consensus <- sum(!is.na(mcols(dmrs)$consensus_seq))
n_unique_motif_consensus <- length(unique(stats::na.omit(mcols(dmrs)$consensus_seq)))
cat(sprintf("Motifs extracted for %d/%d DMRs.\n", n_motif_consensus, length(dmrs)))
cat(sprintf("Unique motif consensus sequences: %d\n", n_unique_motif_consensus))

## ----motif_interactions-------------------------------------------------------
query_with_jaspar <- requireNamespace("JASPAR2024", quietly = TRUE)
interaction_results <- computeDMRsInteraction(
    dmrs = dmrs,
    genome = genome,
    array = array_type,
    min_similarity = 0.8,
    motif_site_flank_size = 5,
    min_component_size = 2,
    query_components_with_jaspar = query_with_jaspar
)

interaction_tbl <- interaction_results$interactions
component_tbl <- interaction_results$components
dmrs <- interaction_results$dmrs

cat(sprintf("Interactions found: %d\n", nrow(interaction_tbl)))
cat(sprintf("Components found: %d\n", nrow(component_tbl)))
if (nrow(component_tbl) > 0) {
    largest_component <- max(component_tbl$size, na.rm = TRUE)
    cat(
        sprintf(
            "Largest component size: %d / %d DMRs (%.2f%%)\n",
            largest_component, length(dmrs), 100 * largest_component / length(dmrs)
        )
    )
}
annotated_components_m <- !is.na(component_tbl$jaspar_names) & nzchar(component_tbl$jaspar_names)
annotated_components <- sum(annotated_components_m)
cat(sprintf("Components with JASPAR annotation: %d / %d\n", annotated_components, nrow(component_tbl)))

if (nrow(interaction_tbl) > 5) {
    knitr::kable(
        head(interaction_tbl[order(-interaction_tbl$sim), ], 5),
        caption = "Top 5 motif-similarity DMR interactions"
    )
}

if (annotated_components > 0) {
    components_tbl <- component_tbl[annotated_components_m, c("component_id", "size", "consensus_seq", "jaspar_names")]
    components_tbl <- components_tbl[order(components_tbl$size, decreasing = TRUE), , drop = FALSE]
    if (nrow(components_tbl) > 0) {
        knitr::kable(
            components_tbl,
            caption = "Motif interaction components with JASPAR annotation"
        )
    }
}

## ----manhattan_plot, fig.height=6, fig.width=14, dpi=150----------------------
plotDMRsManhattan(
    dmrs = dmrs,
    promoter_col = "in_promoter_of",
    gene_body_col = "in_gene_body_of"
)

## ----single_dmr_plot, fig.height=10-------------------------------------------
merged_counts <- suppressWarnings(as.numeric(as.character(mcols(dmrs)$merged_dmrs_num)))
has_extension <- (mcols(dmrs)$start_seed != mcols(dmrs)$start_site) &
    (mcols(dmrs)$end_seed != mcols(dmrs)$end_site)

complex_candidates <- which(is.finite(merged_counts) & merged_counts > 1 & has_extension)
if (length(complex_candidates) == 0) {
    complex_candidates <- which(has_extension)
}
if (length(complex_candidates) == 0) {
    complex_candidates <- 1L
}

selected_dmr_idx <- complex_candidates[1]
cat(
    sprintf(
        "Selected DMR index %d (score=%s, seeds_num=%s, sites_num=%s)\n",
        selected_dmr_idx,
        as.character(mcols(dmrs)$score[selected_dmr_idx]),
        as.character(mcols(dmrs)$seeds_num[selected_dmr_idx]),
        as.character(mcols(dmrs)$sites_num[selected_dmr_idx])
    )
)

plotDMR(
    dmrs = dmrs,
    dmr_index = selected_dmr_idx,
    beta = beta_handler,
    pheno = pheno,
    sample_group_col = "Sample_Group",
    genome = genome,
    array = array_type,
    motif_site_flank_size = 5
)


## ----circos_plot_region, eval=TRUE, fig.height=10, fig.width=14---------------
plotDMRsCircos(
    dmrs = dmrs,
    beta = beta_handler,
    pheno = pheno,
    genome = genome,
    array = array_type,
    sample_group_col = "Sample_Group",
    min_similarity = 0.8,
    motif_site_flank_size = 5,
    max_num_samples_per_group = 5,
    max_dmrs_per_chr = 40,
    max_sites_per_dmr = 5,
    max_components = 120,
    region = "chr11:2Mb-5Mb, chr11:72Mb-83Mb, chr5:1Mb-12Mb",
    degenerate_resolution = 10e2,
    verbose = 1
)

## ----launch_shiny_app, eval=FALSE---------------------------------------------
# launchCMEntViewer("tmp/CMEnt_example")

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