CMEnt End-to-End Example

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

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

Setup

suppressPackageStartupMessages({
    library(CMEnt)
    library(GenomicRanges)
    library(ggplot2)
})

genome <- "hg19"

Load Example 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)))

## Beta matrix dimensions: 43542 sites x 20 samples

cat(sprintf("Phenotype data dimensions: %d samples x %d columns\n", nrow(pheno), ncol(pheno)))

## Phenotype data dimensions: 20 samples x 3 columns

cat(sprintf("Array type: %s\n", array_type))

## Array type: 450K

cat(sprintf("Number of DMPs (seeds): %d\n", nrow(dmps)))

## Number of DMPs (seeds): 23538

cat(sprintf("Head of phenotype data:\n"))

## Head of phenotype data:

print(head(pheno))

##    Sample_Group Age Gender
## T1       cancer  65      m
## T2       cancer  48      m
## T3       cancer  75      m
## T4       cancer  70      f
## T5       cancer  55      m
## T6       cancer  73      m

cat(sprintf("Head of DMPs:\n"))

## Head of DMPs:

print(head(dmps))

##            intercept         pval          f         qval     pval_adj
## cg07461370  2.000713 1.792909e-03 -1.8944789 1.248781e-03 4.992027e-03
## cg04815577  1.866972 3.266108e-03 -2.4763589 2.074106e-03 8.291275e-03
## cg05648010 -1.263073 6.254793e-06  3.7266733 1.242563e-05 4.967167e-05
## cg22807585 -1.029384 3.612295e-04 -0.3214461 3.277151e-04 1.310047e-03
## cg26274580 -2.186916 1.993233e-02 -4.1934979 9.630916e-03 3.849976e-02
## cg25392700 -0.965860 1.313401e-02 -3.8039704 6.750557e-03 2.698547e-02

# 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.

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)))

## DMRs found: 2130

cat(sprintf("DMRs metadata columns: %s\n", paste(colnames(mcols(dmrs)), collapse = ", ")))

## DMRs metadata columns: id, start_seed, end_seed, start_seed_pos, end_seed_pos, seeds_num, stop_connection_reason, connection_corr_pval, seeds, start_site, end_site, upstream_expansion_stop_reason, upstream_sites, upstream_expansion_length, downstream_expansion_stop_reason, downstream_sites, downstream_expansion_length, merged_dmrs_num, sites, supporting_sites_num, sites_num, sites_num_bg, seeds_num_adj, cases_beta, controls_beta, delta_beta, cases_beta_sd, controls_beta_sd, cases_beta_min, cases_beta_max, controls_beta_min, controls_beta_max, sites_cases_beta, sites_controls_beta, sites_delta_beta, sites_cases_beta_sd, sites_controls_beta_sd, sites_cases_beta_min, sites_cases_beta_max, sites_controls_beta_min, sites_controls_beta_max, in_promoter_of, in_gene_body_of, score, cv_accuracy, score_smoothed, segment_id, segment_slope, block_id, pwm, consensus_seq

cat(sprintf(
    "DMRs genomic range example: %s\n",
    paste0(seqnames(dmrs)[1], ":", start(dmrs)[1], "-", end(dmrs)[1])
))

## DMRs genomic range example: chr11:192141-193016

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.

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.

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)))

## Motifs extracted for 2130/2130 DMRs.

cat(sprintf("Unique motif consensus sequences: %d\n", n_unique_motif_consensus))

## Unique motif consensus sequences: 2104

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.

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)))

## Interactions found: 11

cat(sprintf("Components found: %d\n", nrow(component_tbl)))

## Components found: 11

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)
        )
    )
}

## Largest component size: 2 / 2130 DMRs (0.09%)

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)))

## Components with JASPAR annotation: 0 / 11

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

plotDMRsManhattan(
    dmrs = dmrs,
    promoter_col = "in_promoter_of",
    gene_body_col = "in_gene_body_of"
)

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])
    )
)

## Selected DMR index 14 (score=0.657102738173804, seeds_num=9, sites_num=16)

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_dmrs <- dmrs
circos_summary_plotted <- FALSE

tryCatch(
    {
        plotDMRsCircos(
            dmrs = circos_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", # nolint
            degenerate_resolution = 1e6,
            verbose = 1
        )
        circos_summary_plotted <- TRUE
    },
    error = function(e) {
        message("Circos plot generation failed: ", conditionMessage(e))
        NULL
    }
)

if (!circos_summary_plotted) {
    plot.new()
    text(0.5, 0.5, "Circos summary plot unavailable for this dataset.", cex = 1.1)
}

circos_region_plotted <- FALSE

tryCatch(
    {
        plotDMRsCircos(
            dmrs = circos_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 = "chr17:2Mb-5Mb, chr17:72Mb-83Mb, chr19:1Mb-12Mb",
            degenerate_resolution = 10e2,
            verbose = 1
        )
        circos_region_plotted <- TRUE
    },
    error = function(e) {
        message("Circos plot generation failed: ", conditionMessage(e))
        NULL
    }
)

## NULL

if (!circos_region_plotted) {
    plot.new()
    text(0.5, 0.5, "Circos region plot unavailable for selected regions.", cex = 1.1)
}

launchCMEntViewer("tmp/CMEnt_example")

Session Info

sessionInfo()

## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 26.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] DMRsegaldata_1.1.0   ExperimentHub_3.3.2  AnnotationHub_4.3.2 
##  [4] BiocFileCache_3.3.0  dbplyr_2.6.0         ggplot2_4.0.3       
##  [7] GenomicRanges_1.65.4 Seqinfo_1.3.2        IRanges_2.47.5      
## [10] S4Vectors_0.51.10    BiocGenerics_0.59.12 generics_0.1.4      
## [13] CMEnt_0.99.12        BiocStyle_2.41.0    
## 
## loaded via a namespace (and not attached):
##   [1] BiocIO_1.23.3                           
##   [2] bitops_1.1-0                            
##   [3] filelock_1.0.3                          
##   [4] tibble_3.3.1                            
##   [5] R.oo_1.27.1                             
##   [6] XML_3.99-0.25                           
##   [7] DirichletMultinomial_1.55.1             
##   [8] lifecycle_1.0.5                         
##   [9] httr2_1.3.0                             
##  [10] pwalign_1.9.1                           
##  [11] doParallel_1.0.17                       
##  [12] lattice_0.23-1                          
##  [13] backports_1.5.1                         
##  [14] magrittr_2.0.5                          
##  [15] limma_3.99.1                            
##  [16] sass_0.4.10                             
##  [17] rmarkdown_2.32                          
##  [18] jquerylib_0.1.4                         
##  [19] yaml_2.3.12                             
##  [20] otel_0.2.0                              
##  [21] DBI_1.3.0                               
##  [22] buildtools_1.0.0                        
##  [23] RColorBrewer_1.1-3                      
##  [24] abind_1.4-8                             
##  [25] purrr_1.2.2                             
##  [26] R.utils_2.13.0                          
##  [27] RCurl_1.98-1.20                         
##  [28] rappdirs_0.3.4                          
##  [29] circlize_0.4.18                         
##  [30] maketools_1.3.2                         
##  [31] seqLogo_1.79.0                          
##  [32] testthat_3.3.2                          
##  [33] permute_0.9-10                          
##  [34] DelayedMatrixStats_1.35.0               
##  [35] codetools_0.2-20                        
##  [36] DelayedArray_0.39.8                     
##  [37] DT_0.34.0                               
##  [38] tidyselect_1.2.1                        
##  [39] shape_1.4.6.1                           
##  [40] futile.logger_1.4.9                     
##  [41] ggseqlogo_0.2.2                         
##  [42] UCSC.utils_1.9.0                        
##  [43] farver_2.1.2                            
##  [44] matrixStats_1.5.0                       
##  [45] showtext_0.9-8                          
##  [46] GenomicAlignments_1.49.2                
##  [47] jsonlite_2.0.0                          
##  [48] GetoptLong_1.1.1                        
##  [49] iterators_1.0.14                        
##  [50] foreach_1.5.2                           
##  [51] tools_4.6.1                             
##  [52] TxDb.Hsapiens.UCSC.hg19.knownGene_3.22.1
##  [53] TFMPvalue_1.0.0                         
##  [54] Rcpp_1.1.2                              
##  [55] glue_1.8.1                              
##  [56] gridExtra_2.3.1                         
##  [57] SparseArray_1.13.4                      
##  [58] BiocBaseUtils_1.15.1                    
##  [59] xfun_0.61                               
##  [60] MatrixGenerics_1.25.0                   
##  [61] GenomeInfoDb_1.49.1                     
##  [62] dplyr_1.2.1                             
##  [63] HDF5Array_1.41.3                        
##  [64] withr_3.0.3                             
##  [65] formatR_1.14                            
##  [66] BiocManager_1.30.27                     
##  [67] fastmap_1.2.0                           
##  [68] bedr_1.1.5                              
##  [69] rhdf5filters_1.25.4                     
##  [70] caTools_1.18.4                          
##  [71] digest_0.6.39                           
##  [72] R6_2.6.1                                
##  [73] colorspace_2.1-4                        
##  [74] gtools_3.9.5                            
##  [75] dichromat_2.0-1                         
##  [76] RSQLite_3.53.3                          
##  [77] cigarillo_1.3.1                         
##  [78] R.methodsS3_1.8.2                       
##  [79] h5mread_1.5.3                           
##  [80] data.table_1.18.6.1                     
##  [81] rtracklayer_1.73.0                      
##  [82] FNN_1.1.4.1                             
##  [83] httr_1.4.9                              
##  [84] htmlwidgets_1.6.4                       
##  [85] S4Arrays_1.13.2                         
##  [86] TFBSTools_1.51.0                        
##  [87] pkgconfig_2.0.3                         
##  [88] gtable_0.3.6                            
##  [89] blob_1.3.0                              
##  [90] ComplexHeatmap_2.29.0                   
##  [91] S7_0.2.2                                
##  [92] XVector_0.53.0                          
##  [93] sys_3.4.3                               
##  [94] brio_1.1.5                              
##  [95] htmltools_0.5.9                         
##  [96] sysfonts_0.8.9                          
##  [97] strex_2.1.0                             
##  [98] clue_0.3-68                             
##  [99] scales_1.4.0                            
## [100] Biobase_2.73.2                          
## [101] png_0.1-9                               
## [102] knitr_1.52                              
## [103] lambda.r_1.2.4                          
## [104] reshape2_1.4.5                          
## [105] rjson_0.2.23                            
## [106] checkmate_2.3.4                         
## [107] curl_8.0.0                              
## [108] org.Hs.eg.db_3.23.1                     
## [109] showtextdb_3.0                          
## [110] cachem_1.1.0                            
## [111] rhdf5_2.57.20                           
## [112] GlobalOptions_0.1.4                     
## [113] stringr_1.6.0                           
## [114] BiocVersion_3.24.0                      
## [115] parallel_4.6.1                          
## [116] AnnotationDbi_1.75.2                    
## [117] restfulr_0.0.17                         
## [118] pillar_1.11.1                           
## [119] grid_4.6.1                              
## [120] vctrs_0.7.3                             
## [121] beachmat_2.29.3                         
## [122] cluster_2.1.8.3                         
## [123] JASPAR2024_0.99.7                       
## [124] evaluate_1.0.5                          
## [125] bsseq_1.49.2                            
## [126] VennDiagram_1.8.2                       
## [127] GenomicFeatures_1.65.0                  
## [128] cli_3.6.6                               
## [129] locfit_1.5-9.12                         
## [130] compiler_4.6.1                          
## [131] futile.options_1.0.1                    
## [132] Rsamtools_2.29.0                        
## [133] rlang_1.3.0                             
## [134] crayon_1.5.3                            
## [135] labeling_0.4.3                          
## [136] plyr_1.8.9                              
## [137] stringi_1.8.9                           
## [138] gridBase_0.4-7                          
## [139] BiocParallel_1.47.0                     
## [140] Biostrings_2.81.9                       
## [141] Matrix_1.7-6                            
## [142] BSgenome_1.81.1                         
## [143] sparseMatrixStats_1.25.0                
## [144] bit64_4.8.6                             
## [145] Rhdf5lib_2.1.0                          
## [146] KEGGREST_1.53.6                         
## [147] statmod_1.5.2                           
## [148] SummarizedExperiment_1.43.0             
## [149] igraph_2.3.4                            
## [150] memoise_2.0.1                           
## [151] bslib_0.12.0                            
## [152] bit_4.6.0