---
title: "CMEnt End-to-End Example"
author: "CMEnt Package"
date: "`r Sys.Date()`"
output:
  BiocStyle::html_document:
    code_folding: hide
vignette: >
  %\VignetteIndexEntry{CMEnt End-to-End Example}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r 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
)
```

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

## Introduction

CMEnt identifies differentially methylated regions (DMRs) from genomic seeds,
typically differentially methylated positions (DMPs), by expanding seed regions
with correlation-based local connectivity. The package is intended to sit after
established methylation preprocessing and DMP-calling workflows.

CMEnt belongs in Bioconductor because it uses common Bioconductor data
structures such as `GRanges` and `BSseq`, integrates with Bioconductor genome
annotation packages, and provides reproducible examples through packaged data.
It complements packages such as minfi, DSS, DMRcate, and bumphunter: those tools
can preprocess data or identify DMPs/DMRs, while CMEnt focuses on seed-driven
region expansion, optional discriminative DMR scoring, motif extraction,
motif-based DMR interaction summaries, and visualization.

This vignette showcases the core CMEnt workflow, including DMR scoring, motif
extraction, motif-based DMR interactions, and visualizations.

## Installation

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

## Setup

```{r setup}
suppressPackageStartupMessages({
    library(CMEnt)
    library(GenomicRanges)
    library(ggplot2)
})

genome <- "hg19"
```

## Load Example Inputs

```{r 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)

```

## DMR Assembly from Seeds
The `buildDMRs()` function builds DMRs from the provided seeds, which in this
case are DMPs. It performs correlation-based expansion around the seeds while
tracking statistical support and methylation effect size.

To keep the vignette fast and reproducible during package checks, we load a
bundled example result generated with this workflow. The resulting DMRs are
stored in a `GRanges` object with metadata columns for the number of seeds and
sites in each DMR, delta beta between sample groups, and DMR-level `pval`/`qval`
when seed p-values were supplied.
```{r 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
DMRs are primarily summarized by genomic support, DMR-level `pval`/`qval` when
available, and methylation effect size. The optional SVM-derived `score` and
`cv_accuracy` columns are complementary measures of how well the DMR methylation
pattern discriminates sample groups.

Score trajectories are also smoothed per chromosome and segmented, enabling
assignment of localized block IDs that capture local peaks in Manhattan plots.
```{r 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
Motifs are extracted from DMR seed neighborhoods with `extractDMRMotifs()`,
which is called by default from `buildDMRs()`. The function retrieves sequences
around the seeds, applies a flank size of 5 bp, and identifies consensus
sequences.
```{r 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-Based DMR Interactions
Once motifs are extracted, `computeDMRsInteraction()` can identify motif-similar
DMR pairs.

The rationale is that DMRs with similar motifs may be co-regulated or
functionally related, even if they are not in close genomic proximity.
The function computes pairwise similarities between the PWMs of the extracted
motifs,
excluding the center site positions, corrects for microarray design biases,
applies a minimum similarity threshold of 0.8, and identifies connected
components of interacting DMRs.
The resulting components are annotated with JASPAR motifs if available.
```{r 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"
        )
    }
}
```

## Visualizations

### DMR Score Manhattan Plot

A manhattan plot of the DMRs is drawn, showing complementary SVM scores,
their annotation as promoter or gene body DMRs, and translucent block squares
for localized score peaks identified by the scoring block extension. The DMRs
that belong to both promoter and gene body are omitted.
```{r 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
Each DMR can be plotted with `plotDMR()`, which shows DMR site beta values,
seeds, motif-based interactions with other DMRs, and available JASPAR
annotations.

Here we select a DMR that has more than one merged DMR and extension beyond both
ends of the seed region.

```{r 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
)

```

### Summary Circos Plot
A circos plot consisting of extracted CMEnt information is drawn along the
genome or a specified region.
The information includes:

- The DMRs and their complementary SVM scores (at most 40 DMRs per chromosome
  are shown)
- Selected beta values of the DMRs (at most 5 sites per DMR, scoring by
  absolute delta-beta, are shown)
- The motif-based interactions between DMRs, with links direction going from
  higher scoring to lower scoring (at most 120 interactions groups are shown)
- The JASPAR annotation of the motif-based groups, if available, on the legend.

It is also possible to plot circos output for one or more regions to assess
connectivity. The reported score in the legend is computed from the visible
DMRs of each component. The selected intervals lie on chromosomes 5 and 11,
which are covered by the example dataset.
```{r 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
)
```

### CMEnt interactive Shiny App
The CMEnt interactive Shiny app can be launched with the `launchCMEntViewer`
function,
which allows for interactive exploration of the DMRs and the aforementioned
plots. In order to use it, the `output_prefix` argument should be provided in
the `buildDMRs` function, as the app looks for the output files with this
prefix.
```{r launch_shiny_app, eval=FALSE}
launchCMEntViewer("tmp/CMEnt_example")
```



## Session Info

```{r session_info}
sessionInfo()
```