Contents

0.1 Introduction

GSE244696 (doi: https://doi.org/10.1016/j.trsl.2024.05.005) is a PDAC MethylationEPIC study comparing tumor-cell enriched (“Tumor”) and stroma-cell enriched (“Stroma”) compartments. This notebook uses CMEnt in a meta-analysis mode: the published DMP list (Supplementary table 1) is treated as the external evidence, and CMEnt uses the available beta matrix from GEO to assemble coherent DMRs around those known signals. The only external DMR comparison retained here is the Bumphunter DMR table (Supplementary table 2) provided by the original researchers.

0.2 Setup

ensure_cran_packages <- function(pkgs) {
    repos <- getOption("repos")
    if (is.null(repos) || is.na(repos["CRAN"]) || repos["CRAN"] %in% c("@CRAN@", "")) {
        options(repos = c(CRAN = "https://cloud.r-project.org"))
    }
    missing <- pkgs[!vapply(pkgs, requireNamespace, logical(1), quietly = TRUE)]
    if (length(missing) > 0L) {
        install.packages(missing)
    }
}

ensure_bioc_packages <- function(pkgs) {
    if (!requireNamespace("BiocManager", quietly = TRUE)) {
        ensure_cran_packages("BiocManager")
    }
    missing <- pkgs[!vapply(pkgs, requireNamespace, logical(1), quietly = TRUE)]
    if (length(missing) > 0L) {
        BiocManager::install(missing, update = FALSE, ask = FALSE)
    }
}
suppressWarnings(
    {
        ensure_cran_packages(c("readxl", "data.table", "devtools", "dplyr", "ggplot2", "tidyr", "scales", "svglite"))
        ensure_bioc_packages(c(
            "BiocStyle", "minfi", "GenomicRanges", "annotatr", "missMethyl",
            "org.Hs.eg.db", "TxDb.Hsapiens.UCSC.hg19.knownGene",
            "IlluminaHumanMethylationEPICanno.ilm10b4.hg19", "clusterProfiler",
            "ReactomePA", "pathview", "KEGGREST"
        ))
    }
)
.loadCMEnt <- function() {
    root <- if (file.exists("NAMESPACE")) "." else if (file.exists(file.path("..", "NAMESPACE"))) ".." else NULL
    if (!is.null(root) && requireNamespace("pkgload", quietly = TRUE)) {
        ok <- tryCatch({
            pkgload::load_all(root, quiet = TRUE)
            TRUE
        }, error = function(e) {
            message("Could not load local CMEnt source: ", conditionMessage(e))
            FALSE
        })
        if (ok) return(invisible(TRUE))
    }
    suppressPackageStartupMessages(library(CMEnt))
}
.loadCMEnt()
suppressPackageStartupMessages({
    library(data.table)
    library(dplyr)
    library(ggplot2)
    library(GenomicRanges)
    library(minfi)
})

getBenchmarkFile <- function(filename) {
    folder <- "metastudyEpic"
    if (getwd() == "notebooks") {
        benchmark_file <- file.path("cached_results", folder, filename)
    } else if ("NAMESPACE" %in% dir()) {
        benchmark_file <- file.path("notebooks", "cached_results", folder, filename)
    } else {
        benchmark_file <- file.path("cached_results", folder, filename)
    }
    dir.create(dirname(benchmark_file), showWarnings = FALSE, recursive = TRUE)
    benchmark_file
}

short_label <- function(x, width = 65L) {
    x <- as.character(x)
    ifelse(nchar(x) > width, paste0(substr(x, 1L, width - 3L), "..."), x)
}

has_text <- function(x) {
    !is.na(x) & nzchar(as.character(x))
}

theme_meta <- function() {
    theme_minimal(base_size = 11) +
        theme(
            panel.grid.minor = element_blank(),
            axis.text.x = element_text(angle = 35, hjust = 1),
            legend.position = "bottom"
        )
}

method_palette <- c(CMEnt = "#0072B2", `Published Bumphunter` = "#D55E00")
array_type <- "EPIC"
genome <- "hg19"

0.3 Load Data and Published Evidence

beta_file <- getBenchmarkFile("GSE244696_beta_values.tsv.gz")
samplesheet <- getBenchmarkFile("GSE244696_samplesheet.tsv")
dmp_file <- getBenchmarkFile("GSE244696_DMPs.xlsx")
published_dmr_file <- getBenchmarkFile("GSE244696_DMRs.xlsx")


zenodo_base_url <- "https://zenodo.org/records/20592944/files/"
if (!file.exists(beta_file)) {
    curl::curl_download(
        url = paste0(zenodo_base_url, basename(beta_file), "?download=1"),
        destfile = beta_file
    )
}
if (!file.exists(samplesheet)) {
    curl::curl_download(
        url = paste0(zenodo_base_url, basename(samplesheet), "?download=1"),
        destfile = samplesheet
    )
}
if (!file.exists(dmp_file)) {
    curl::curl_download(
        url = paste0(zenodo_base_url, basename(dmp_file), "?download=1"),
        destfile = dmp_file
    )
}
if (!file.exists(published_dmr_file)) {
    curl::curl_download(
        url = paste0(zenodo_base_url, basename(published_dmr_file), "?download=1"),
        destfile = published_dmr_file
    )
}


pheno <- utils::read.table(samplesheet, header = TRUE, stringsAsFactors = FALSE, sep = "\t", check.names = FALSE)
if (all(rownames(pheno) == as.character(seq_len(nrow(pheno)))) && nzchar(names(pheno)[1])) {
    rownames(pheno) <- pheno[[1]]
}
if (all(rownames(pheno) == as.character(seq_len(nrow(pheno)))) && !nzchar(names(pheno)[1])) {
    rownames(pheno) <- pheno[[1]]
}
pheno$casecontrol <- pheno$tissue.ch1 == "Tumor"

beta_handler <- CMEnt::getBetaHandler(beta_file, array = array_type, genome = genome)
beta_mat <- as.matrix(beta_handler$getBeta())
locs <- beta_handler$getBetaLocs()
beta_row_names <- beta_handler$getBetaRowNames()

common_samples <- intersect(colnames(beta_mat), rownames(pheno))
beta_mat <- beta_mat[, common_samples, drop = FALSE]
pheno <- pheno[common_samples, , drop = FALSE]
case_samples <- rownames(pheno)[pheno$tissue.ch1 == "Tumor"]
control_samples <- rownames(pheno)[pheno$tissue.ch1 == "Stroma"]

reported_dmps_raw <- suppressWarnings(suppressMessages(readxl::read_xlsx(
    dmp_file,
    sheet = "Suppl_Table 1",
    range = readxl::cell_cols(1:21),
    .name_repair = "unique"
)))
reported_dmps <- data.frame(
    site_id = as.character(reported_dmps_raw[[1]]),
    pval = reported_dmps_raw[["Stroma_to_Tumor.P.Value"]],
    qval = reported_dmps_raw[["Stroma_to_Tumor.adj.P.Val"]],
    logFC_reported = reported_dmps_raw[["Stroma_to_Tumor.logFC"]],
    delta_beta_reported = reported_dmps_raw[["Stroma_to_Tumor.deltaBeta"]],
    stroma_avg_reported = reported_dmps_raw[["Stroma_to_Tumor.Stroma_AVG"]],
    tumor_avg_reported = reported_dmps_raw[["Stroma_to_Tumor.Tumor_AVG"]],
    chr_reported = paste0("chr", reported_dmps_raw[["Stroma_to_Tumor.CHR"]]),
    mapinfo_reported = reported_dmps_raw[["Stroma_to_Tumor.MAPINFO"]],
    gene_reported = reported_dmps_raw[["Stroma_to_Tumor.gene"]],
    feature_reported = reported_dmps_raw[["Stroma_to_Tumor.feature"]],
    cgi_reported = reported_dmps_raw[["Stroma_to_Tumor.cgi"]],
    stringsAsFactors = FALSE
)
reported_dmps <- reported_dmps[!is.na(reported_dmps$site_id) & nzchar(reported_dmps$site_id), ]
reported_dmps <- reported_dmps[!duplicated(reported_dmps$site_id), ]

sig_dmps <- reported_dmps[reported_dmps$site_id %in% beta_row_names, , drop = FALSE]
sig_dmps$delta_beta_available <- rowMeans(beta_mat[sig_dmps$site_id, case_samples, drop = FALSE], na.rm = TRUE) -
    rowMeans(beta_mat[sig_dmps$site_id, control_samples, drop = FALSE], na.rm = TRUE)

dmp_locs <- locs[sig_dmps$site_id, , drop = FALSE]
dmp_gr <- makeGRangesFromDataFrame(
    data.frame(
        chr = dmp_locs$chr,
        start = dmp_locs$start,
        end = dmp_locs$start + 1L,
        site_id = sig_dmps$site_id,
        pval = sig_dmps$pval,
        qval = sig_dmps$qval
    ),
    keep.extra.columns = TRUE
)
names(dmp_gr) <- sig_dmps$site_id

published_dmr_raw <- suppressWarnings(suppressMessages(readxl::read_xlsx(
    published_dmr_file,
    sheet = "Suppl_DMR",
    range = readxl::cell_cols(1:32),
    .name_repair = "unique"
)))
published_bumphunter <- makeGRangesFromDataFrame(
    as.data.frame(published_dmr_raw, check.names = FALSE),
    seqnames.field = "seqnames",
    start.field = "start",
    end.field = "end",
    keep.extra.columns = TRUE
)
names(published_bumphunter) <- make.unique(as.character(published_dmr_raw[[1]]))
genome(published_bumphunter) <- genome
mcols(published_bumphunter)$dmr_id <- names(published_bumphunter)

input_summary <- data.frame(
    Metric = c(
        "Samples",
        "Tumor samples",
        "Stroma samples",
        "Published DMPs",
        "Published DMPs present in beta matrix",
        "Published Bumphunter DMRs"
    ),
    Value = c(
        nrow(pheno),
        length(case_samples),
        length(control_samples),
        nrow(reported_dmps),
        nrow(sig_dmps),
        length(published_bumphunter)
    )
)
knitr::kable(input_summary, row.names = FALSE, caption = "Meta-study input summary")

Table 1: Meta-study input summary
Metric Value
Samples 26
Tumor samples 13
Stroma samples 13
Published DMPs 25410
Published DMPs present in beta matrix 22667
Published Bumphunter DMRs 29 
dmp_direction <- data.frame(
    Direction = c("Higher in Tumor", "Higher in Stroma"),
    Count = c(sum(sig_dmps$delta_beta_available > 0, na.rm = TRUE), sum(sig_dmps$delta_beta_available < 0, na.rm = TRUE))
)
knitr::kable(dmp_direction, row.names = FALSE, caption = "Direction of published DMPs")

Table 2: Direction of published DMPs
Direction Count
Higher in Tumor 17566
Higher in Stroma 5101 
sig_dmps_delta_beta_medians <- sig_dmps |>
    group_by(Direction = ifelse(delta_beta_available > 0, "Higher in Tumor", "Higher in Stroma")) |>
    summarise(
        Median = median(delta_beta_available, na.rm = TRUE),
        y = 1.1 * max(hist(delta_beta_available, breaks = 60, plot = FALSE)$counts, na.rm = TRUE)
    )
ggplot(sig_dmps, aes(x = delta_beta_available)) +
    geom_histogram(bins = 60, fill = "#0072B2", color = "white", linewidth = 0.15) +
    geom_vline(xintercept = 0, linetype = "dashed", color = "grey35") +
    geom_text(
        data = sig_dmps_delta_beta_medians,
        aes(x = Median, y = y, label = paste0("Median: ", round(Median, 3))),
        inherit.aes = FALSE,
        color = "black",
        fontface = "bold",
        size = 3.5
    ) +
    theme_meta() +
    labs(
        title = "Effect Size for Published DMPs",
        x = "Delta beta", y = "Published DMPs"
    )
Histogram shows published DMP counts across delta beta values. A dashed vertical line marks zero; text labels give medians for DMPs with higher methylation in Tumor or Stroma.

Figure 1: Distribution of available Tumor-minus-Stroma delta beta values for the published GSE244696 DMPs

0.4 CMEnt Meta-Analysis DMRs

cment_file <- getBenchmarkFile("dmrs.CMEnt.meta_reported_dmps.rds")
if (!file.exists(cment_file)) {
    cment_njobs <- max(1L, min(8L, parallel::detectCores(logical = TRUE) - 1L))
    options("CMEnt.verbose" = 1)
    options("CMEnt.njobs" = cment_njobs)
    cat(sprintf("CMEnt parallel jobs: %d\n", cment_njobs))
    start_time <- Sys.time()
    dmrs_cment <- CMEnt::buildDMRs(
        beta = beta_handler,
        seeds = sig_dmps,
        pheno = pheno,
        seeds_id_col = "site_id",
        sample_group_col = "tissue.ch1",
        covariates = "predictedSex",
        max_bridge_seeds_gaps = 1,
        max_bridge_extension_gaps = 3,
        min_seeds = 2,
        min_sites = 3,
        ext_site_delta_beta = 0.2,
        max_lookup_dist = 1000,
        max_pval = 0.05,
        entanglement = "weak",
        score_dmrs = TRUE,
        annotate_with_genes = TRUE,
        extract_motifs = TRUE,
        njobs = cment_njobs,
        verbose = 1
    )
    time_cment <- Sys.time() - start_time
    saveRDS(list(dmrs = dmrs_cment, time = time_cment), cment_file)
} else {
    ret <- readRDS(cment_file)
    dmrs_cment <- ret$dmrs
    time_cment <- ret$time
}

cat(sprintf("CMEnt found %d DMRs in %d seconds.\n", length(dmrs_cment), round(as.numeric(time_cment, units = "secs"))))

CMEnt found 2104 DMRs in 177 seconds.

if (is.null(dmrs_cment)) {
    dmrs_cment <- GRanges()
}
if (length(dmrs_cment) > 0L && !"in_promoter_of" %in% colnames(mcols(dmrs_cment))) {
    dmrs_cment <- CMEnt::annotateDMRsWithGenes(dmrs_cment, genome = genome)
}
published_bumphunter <- CMEnt::annotateDMRsWithGenes(published_bumphunter, genome = genome)
methods_list <- list(CMEnt = dmrs_cment, `BH` = published_bumphunter)
dmr_effect_df <- bind_rows(
    data.frame(
        Region = "Promoter-overlapping DMRs",
        Width = width(dmrs_cment),
        Delta_beta = as.numeric(mcols(dmrs_cment)$delta_beta_promoter)
    ),
    data.frame(
        Region = "Gene-body-overlapping DMRs",
        Width = width(dmrs_cment),
        Delta_beta = as.numeric(mcols(dmrs_cment)$delta_beta_gene_body)
    )
) |>
    filter(is.finite(Delta_beta)) |>
    mutate(
        Region = factor(Region, levels = unique(Region)),
        Direction = factor(
            ifelse(Delta_beta > 0, "Hyper", "Hypo"),
            levels = c("Hyper", "Hypo")
        )
    )

library(dplyr)
library(ggplot2)
library(scales)

ignore_abs_delta <- !is.finite(dmr_effect_df$Delta_beta) | abs(dmr_effect_df$Delta_beta) < 0.05
dmr_effect_df <- dmr_effect_df[!ignore_abs_delta, ]
# Summary per Region and Direction
dmr_summary <- dmr_effect_df |>
    group_by(Region, Direction) |>
    summarise(
        median_delta_beta = median(Delta_beta, na.rm = TRUE),
        median_width_bp = median(Width, na.rm = TRUE),
        n_dmrs = n(),
        .groups = "drop"
    )

mood_median_p_value <- function(width, direction) {
    widths <- split(width, direction)
    if (!all(c("Hyper", "Hypo") %in% names(widths)) || any(lengths(widths[c("Hyper", "Hypo")]) < 2L)) {
        return(NA_real_)
    }
    median_width <- median(width, na.rm = TRUE)
    keep <- is.finite(width) & width != median_width
    tab <- table(
        factor(direction[keep], levels = c("Hypo", "Hyper")),
        factor(width[keep] > median_width, levels = c(FALSE, TRUE))
    )
    if (any(rowSums(tab) == 0L) || any(colSums(tab) == 0L)) {
        return(NA_real_)
    }
    suppressWarnings(stats::chisq.test(tab, correct = FALSE)$p.value)
}

p_to_stars <- function(p) {
    case_when(
        is.na(p) ~ "",
        p < 0.001 ~ "***",
        p < 0.01 ~ "**",
        p < 0.05 ~ "*",
        TRUE ~ ""
    )
}

x_range <- range(dmr_effect_df$Delta_beta, na.rm = TRUE)
x_min <- x_range[[1]]

region_levels <- unique(dmr_effect_df$Region)
dmr_mood <- dmr_summary |>
    group_by(Region) |>
    summarise(
        star_y = exp(mean(log(median_width_bp))), .groups = "drop",
        segy_min = min(median_width_bp), segy_max = max(median_width_bp)
    ) |>
    left_join(
        dmr_effect_df |>
            group_by(Region) |>
            summarise(mood_p = mood_median_p_value(Width, Direction), .groups = "drop"),
        by = "Region"
    ) |>
    mutate(
        stars = p_to_stars(mood_p),
        star_x = x_min - 0.05 * abs(x_min),
        hjust = 0.5
    )

# Define two y-levels for the two summary bars
y_max <- max(dmr_effect_df$Width, na.rm = TRUE)

dmr_summary <- dmr_summary |>
    mutate(
        bar_y = case_when(
            Direction == "Hyper" ~ y_max * 2.0,
            Direction == "Hypo" ~ y_max * 3.2
        ),
        label_y = bar_y * 1.12
    )

# Expand y-range so the bars and labels fit
y_upper <- max(dmr_summary$label_y) * 1.15

ggplot(dmr_effect_df, aes(x = Delta_beta, y = Width, color = Direction)) +
    geom_hline(
        data = dmr_summary,
        aes(yintercept = median_width_bp, color = Direction),
        inherit.aes = FALSE,
        linewidth = 0.2
    ) +
    geom_vline(xintercept = 0, linewidth = 0.2, color = "grey70") +
    geom_point(alpha = 0.55, size = 1.2) +
    geom_segment(
        data = dmr_mood,
        aes(x = star_x + 0.05 * abs(x_min), xend = star_x + 0.05 * abs(x_min), y = segy_min, yend = segy_max),
        inherit.aes = FALSE,
        linewidth = 0.4,
        arrow = arrow(length = unit(0.15, "cm"), type = "closed", ends = "both"),
        color = "black",
        linetype = "dashed"
    ) +
    geom_text(
        data = dmr_mood,
        aes(x = star_x, y = star_y, label = stars, hjust = hjust),
        inherit.aes = FALSE,
        angle = 90,
        color = "black",
        fontface = "bold",
        size = 5
    ) +

    # Separate horizontal bars for hypo and hyper
    geom_segment(
        data = dmr_summary,
        aes(
            x = 0,
            xend = median_delta_beta,
            y = bar_y,
            yend = bar_y,
            color = Direction
        ),
        inherit.aes = FALSE,
        linewidth = 5,
        lineend = "butt"
    ) +

    # Labels for each bar
    geom_text(
        data = dmr_summary,
        aes(
            x = median_delta_beta,
            y = label_y,
            label = paste0(
                "median Δβ = ", round(median_delta_beta, 3),
                "\nmedian width = ", round(median_width_bp, 1), " bp",
                "\nDMRs = ", n_dmrs
            )
        ),
        inherit.aes = FALSE,
        color = "black",
        size = 4,
        hjust = ifelse(dmr_summary$median_delta_beta >= 0, -0.05, 1.05)
    ) +
    facet_wrap(~Region, nrow = 1) +
    scale_y_log10(
        limits = c(min(dmr_effect_df$Width[dmr_effect_df$Width > 0], na.rm = TRUE), y_upper),
        labels = comma
    ) +
    scale_color_manual(
        values = c(
            "Hyper" = "#B8405E",
            "Hypo" = "#247BA0"
        )
    ) +
    coord_cartesian(clip = "off") +
    theme_minimal(base_size = 12) +
    theme(legend.position = "bottom", strip.text.x = element_text(size = 12)) +
    labs(
        title = "DMR Effect Size and Width by Gene Context",
        x = "Context-specific delta beta (Tumor - Stroma)",
        y = "DMR width (bp, log scale)",
        color = ""
    )
Two scatter panels plot context-specific delta beta against DMR width on a logarithmic axis. Red points indicate hypermethylation and blue points hypomethylation; bars and labels summarize medians and counts, and stars mark Mood median-test significance.

(#fig:cment_dmr_scatter_plot)DMR width versus Tumor-minus-Stroma delta beta in promoter and gene-body contexts, restricted to absolute delta beta of at least 0.05. 

metadata_numeric_summary <- function(gr, fields) {
    fields <- intersect(fields, colnames(mcols(gr)))
    if (length(fields) == 0L || length(gr) == 0L) {
        return(data.frame())
    }
    base::do.call(rbind, lapply(fields, function(field) {
        x <- suppressWarnings(as.numeric(mcols(gr)[[field]]))
        data.frame(
            Field = field,
            Min = min(x, na.rm = TRUE),
            Median = stats::median(x, na.rm = TRUE),
            Mean = mean(x, na.rm = TRUE),
            Max = max(x, na.rm = TRUE)
        )
    }))
}

dmr_summary <- data.frame(
    Set = names(methods_list),
    DMRs = vapply(methods_list, length, integer(1)),
    Median_width_bp = vapply(methods_list, function(gr) stats::median(width(gr)), numeric(1)),
    Total_covered_bp = vapply(methods_list, function(gr) sum(width(reduce(gr))), numeric(1)),
    stringsAsFactors = FALSE
)
knitr::kable(dmr_summary, row.names = FALSE, caption = "DMR set summary")

Table 3: DMR set summary
Set DMRs Median_width_bp Total_covered_bp
CMEnt 2104 1095.5 3097452
BH 29 654.0 19371 
cment_fields <- c("n_sites", "n_seeds", "pval", "qval", "score", "cv_accuracy", "mean_delta_beta", "median_delta_beta")
cment_meta_summary <- metadata_numeric_summary(dmrs_cment, cment_fields)
if (nrow(cment_meta_summary) > 0L) {
    knitr::kable(cment_meta_summary, digits = 3, row.names = FALSE, caption = "CMEnt numeric metadata summary")
}

Table 4: CMEnt numeric metadata summary
Field Min Median Mean Max
pval 0.000 0.000 0.000 0.000
qval 0.000 0.000 0.000 0.000
score 0.357 0.475 0.475 0.612
cv_accuracy 0.077 0.500 0.494 0.846 
width_df <- base::do.call(rbind, lapply(names(methods_list), function(method) {
    data.frame(Method = method, Width = width(methods_list[[method]]))
}))
ggplot(width_df, aes(x = Method, y = Width, fill = Method)) +
    geom_boxplot(outlier.alpha = 0.35, width = 0.55) +
    scale_y_log10(labels = scales::comma) +
    scale_fill_manual(values = method_palette, drop = FALSE) +
    theme_meta() +
    labs(title = "DMR Width Distribution", x = "", y = "Width (bp, log10 scale)")
Boxplots compare DMR widths in base pairs on a logarithmic axis for CMEnt and the supplied Bumphunter regions.

Figure 2: Width distributions of CMEnt and supplied Bumphunter DMRs in GSE244696

0.5 CpG and Gene Context

gene_context_df <- data.frame(
    Context = c("Promoter", "Gene body", "Promoter and gene body", "No annotated gene overlap"),
    DMRs = c(
        sum(has_text(mcols(dmrs_cment)$in_promoter_of)),
        sum(has_text(mcols(dmrs_cment)$in_gene_body_of)),
        sum(has_text(mcols(dmrs_cment)$in_promoter_of) & has_text(mcols(dmrs_cment)$in_gene_body_of)),
        sum(!has_text(mcols(dmrs_cment)$in_promoter_of) & !has_text(mcols(dmrs_cment)$in_gene_body_of))
    )
) |>
    mutate(Percentage = round(100 * DMRs / length(dmrs_cment), 1))

knitr::kable(gene_context_df, caption = "Gene context of CMEnt DMRs")
(#tab:gene_context)Gene context of CMEnt DMRs
Context DMRs Percentage
Promoter 1455 69.2
Gene body 1044 49.6
Promoter and gene body 620 29.5
No annotated gene overlap 225 10.7 
ggplot(gene_context_df, aes(x = reorder(Context, DMRs), y = DMRs, fill = Context)) +
    geom_col(width = 0.7) +
    coord_flip() +
    scale_y_continuous(labels = scales::comma) +
    theme_minimal() +
    theme(legend.position = "none") +
    labs(title = "DMR Gene Context", x = "", y = "DMRs")
Horizontal bars show DMR counts for promoter overlap, gene-body overlap, both overlaps, and no annotated gene overlap.

(#fig:gene_context_plot)Gene context of CMEnt DMRs. Promoter and gene-body overlap categories can include the same DMR. 

split_ids <- function(x) {
    ids <- unlist(strsplit(stats::na.omit(as.character(x)), ","), use.names = FALSE)
    unique(ids[nzchar(ids)])
}

if (requireNamespace("IlluminaHumanMethylationEPICanno.ilm10b4.hg19", quietly = TRUE)) {
    library(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
    anno <- as.data.frame(minfi::getAnnotation(IlluminaHumanMethylationEPICanno.ilm10b4.hg19))
    anno$probe_id <- rownames(anno)

    dmr_site_ids <- split_ids(mcols(dmrs_cment)$sites)
    dmr_site_context <- anno[intersect(dmr_site_ids, rownames(anno)), c("probe_id", "Relation_to_Island")]
    dmp_site_context <- anno[intersect(sig_dmps$site_id, rownames(anno)), c("probe_id", "Relation_to_Island")]

    site_context_df <- bind_rows(
        dmp_site_context |> mutate(Source = "Significant DMPs"),
        dmr_site_context |> mutate(Source = "Sites inside DMRs")
    ) |>
        dplyr::count(Source, Relation_to_Island, name = "Sites") |>
        group_by(Source) |>
        mutate(Percentage = round(100 * Sites / sum(Sites), 1)) |>
        ungroup()

    knitr::kable(site_context_df, caption = "CpG island context of significant DMPs and sites inside DMRs")
}
(#tab:cpg_island_context)CpG island context of significant DMPs and sites inside DMRs
Source Relation_to_Island Sites Percentage
Significant DMPs Island 3101 13.7
Significant DMPs N_Shelf 684 3.0
Significant DMPs N_Shore 1934 8.5
Significant DMPs OpenSea 14907 65.8
Significant DMPs S_Shelf 558 2.5
Significant DMPs S_Shore 1483 6.5
Sites inside DMRs Island 5917 37.4
Sites inside DMRs N_Shelf 358 2.3
Sites inside DMRs N_Shore 2884 18.2
Sites inside DMRs OpenSea 4355 27.5
Sites inside DMRs S_Shelf 298 1.9
Sites inside DMRs S_Shore 2014 12.7 
if (exists("site_context_df")) {
    context_levels <- site_context_df |>
        group_by(Relation_to_Island) |>
        summarise(Max_sites = max(Sites), .groups = "drop") |>
        arrange(Max_sites) |>
        pull(Relation_to_Island)
    plot_context_df <- site_context_df |>
        mutate(
            Relation_to_Island = factor(Relation_to_Island, levels = context_levels),
            Source = factor(Source, levels = c("Significant DMPs", "Sites inside DMRs"))
        )
    diff_df <- plot_context_df |>
        dplyr::select(Source, Relation_to_Island, Sites) |>
        tidyr::pivot_wider(names_from = Source, values_from = Sites, values_fill = 0) |>
        dplyr::mutate(
            Difference = `Sites inside DMRs` - `Significant DMPs`,
            Label = paste0("Delta=", scales::comma(Difference)),
            Label_x = pmax(`Significant DMPs`, `Sites inside DMRs`)
        )

    ggplot(plot_context_df, aes(x = Sites, y = Relation_to_Island, fill = Source)) +
        geom_col(data = dplyr::filter(plot_context_df, Source == "Significant DMPs"), width = 0.78, alpha = 0.45) +
        geom_col(data = dplyr::filter(plot_context_df, Source == "Sites inside DMRs"), width = 0.48, alpha = 0.9) +
        geom_text(data = diff_df, aes(x = Label_x, y = Relation_to_Island, label = Label), inherit.aes = FALSE, hjust = -0.05, size = 3) +
        scale_x_continuous(labels = scales::comma, expand = expansion(mult = c(0, 0.2))) +
        scale_fill_manual(values = c("Significant DMPs" = "grey65", "Sites inside DMRs" = "#0072B2")) +
        theme_minimal() +
        theme(legend.position = "bottom") +
        labs(title = "CpG Island Context", x = "Sites", y = "", fill = "")
}
Overlaid horizontal bars compare site counts by CpG island context: grey for significant DMPs and blue for sites inside DMRs. Labels show the count differences.

(#fig:cpg_island_context_plot)CpG island context of significant DMPs and sites inside CMEnt DMRs.

0.6 Comparison With the Supplied Bumphunter DMRs

overlap_hits <- findOverlaps(dmrs_cment, published_bumphunter, ignore.strand = TRUE)
cment_overlapping <- unique(queryHits(overlap_hits))
published_overlapping <- unique(subjectHits(overlap_hits))

intersection_bp <- sum(width(intersect(reduce(dmrs_cment), reduce(published_bumphunter), ignore.strand = TRUE)))
union_bp <- sum(width(union(reduce(dmrs_cment), reduce(published_bumphunter), ignore.strand = TRUE)))
overlap_summary <- data.frame(
    Metric = c(
        "CMEnt DMRs overlapping supplied Bumphunter",
        "Supplied Bumphunter DMRs overlapping CMEnt",
        "Fraction of CMEnt DMRs overlapping supplied Bumphunter",
        "Fraction of supplied Bumphunter DMRs overlapping CMEnt",
        "Base-pair Jaccard index"
    ),
    Value = c(
        length(cment_overlapping),
        length(published_overlapping),
        length(cment_overlapping) / max(1, length(dmrs_cment)),
        length(published_overlapping) / max(1, length(published_bumphunter)),
        intersection_bp / max(1, union_bp)
    )
)
knitr::kable(overlap_summary, digits = 3, row.names = FALSE, caption = "CMEnt vs supplied Bumphunter overlap")

Table 5: CMEnt vs supplied Bumphunter overlap
Metric Value
CMEnt DMRs overlapping supplied Bumphunter 29.000
Supplied Bumphunter DMRs overlapping CMEnt 29.000
Fraction of CMEnt DMRs overlapping supplied Bumphunter 0.014
Fraction of supplied Bumphunter DMRs overlapping CMEnt 1.000
Base-pair Jaccard index 0.006 
capture_by_set <- data.frame(
    Set = names(methods_list),
    Captured_DMPs = vapply(methods_list, function(gr) sum(overlapsAny(dmp_gr, gr, ignore.strand = TRUE)), integer(1)),
    Total_DMPs = length(dmp_gr),
    stringsAsFactors = FALSE
)
capture_by_set$Captured_Pct <- 100 * capture_by_set$Captured_DMPs / capture_by_set$Total_DMPs
knitr::kable(capture_by_set, digits = 1, row.names = FALSE, caption = "Published DMPs captured by each DMR set")

Table 6: Published DMPs captured by each DMR set
Set Captured_DMPs Total_DMPs Captured_Pct
CMEnt 6739 22667 29.7
BH 238 22667 1.0 
ggplot(capture_by_set, aes(x = Set, y = Captured_Pct, fill = Set)) +
    geom_col(width = 0.6) +
    geom_text(aes(label = paste0(scales::comma(Captured_DMPs), " probes")), vjust = -0.4, size = 3.3) +
    scale_fill_manual(values = method_palette, drop = FALSE) +
    scale_y_continuous(labels = scales::percent_format(scale = 1), limits = c(0, max(capture_by_set$Captured_Pct) * 1.18)) +
    theme_meta() +
    labs(title = "Capture of Published DMPs", x = "", y = "Published DMPs captured (%)")
Bars compare the percentage of published DMPs overlapping each DMR set. Labels above the bars give the captured probe counts.

Figure 3: Percentage of published GSE244696 DMPs captured by each DMR set

best_overlap_table <- data.frame()
if (length(overlap_hits) > 0L) {
    qi <- queryHits(overlap_hits)
    si <- subjectHits(overlap_hits)
    inter_width <- pmax(
        0,
        pmin(end(dmrs_cment)[qi], end(published_bumphunter)[si]) -
            pmax(start(dmrs_cment)[qi], start(published_bumphunter)[si]) + 1
    )
    pair_jaccard <- inter_width / (width(dmrs_cment)[qi] + width(published_bumphunter)[si] - inter_width)
    best_overlap_table <- data.frame(
        CMEnt_region = paste0(as.character(seqnames(dmrs_cment))[qi], ":", start(dmrs_cment)[qi], "-", end(dmrs_cment)[qi]),
        Bumphunter_region = paste0(as.character(seqnames(published_bumphunter))[si], ":", start(published_bumphunter)[si], "-", end(published_bumphunter)[si]),
        Bumphunter_id = names(published_bumphunter)[si],
        Overlap_bp = inter_width,
        Pair_Jaccard = pair_jaccard,
        Bumphunter_fwer = suppressWarnings(as.numeric(mcols(published_bumphunter)$fwer[si])),
        Bumphunter_gene = as.character(mcols(published_bumphunter)$genename[si]),
        stringsAsFactors = FALSE
    )
    best_overlap_table <- best_overlap_table |>
        arrange(CMEnt_region, desc(Pair_Jaccard)) |>
        group_by(CMEnt_region) |>
        slice_head(n = 1) |>
        ungroup() |>
        arrange(desc(Pair_Jaccard))
    knitr::kable(
        head(best_overlap_table, 15),
        digits = 3,
        row.names = FALSE,
        caption = "Best supplied Bumphunter match for each overlapping CMEnt DMR"
    )
} else {
    cat("No direct genomic overlap was found between CMEnt and the supplied Bumphunter DMRs.\n")
}

Table 7: Best supplied Bumphunter match for each overlapping CMEnt DMR
CMEnt_region Bumphunter_region Bumphunter_id Overlap_bp Pair_Jaccard Bumphunter_fwer Bumphunter_gene
chr6:30139478-30140325 chr6:30139478-30140231 DMR_19 754 0.889 0 TRIM15
chr8:17270347-17271215 chr8:17270604-17271215 DMR_22 612 0.704 0 MTMR7
chr4:4859935-4860655 chr4:4859772-4860547 DMR_6 613 0.693 0 MSX1
chr19:57183016-57183342 chr19:57182816-57183342 DMR_14 327 0.620 0 ZNF835
chr19:58545001-58546307 chr19:58545122-58545837 DMR_13 716 0.548 0 ZSCAN1
chr16:809476-811347 chr16:810365-811347 DMR_23 983 0.525 0 MSLN
chr5:127873228-127874587 chr5:127873283-127873980 DMR_3 698 0.513 0 FBN2
chr5:140810051-140811642 chr5:140810106-140810920 DMR_5 815 0.512 0 PCDHGA8
chr7:27280585-27282112 chr7:27280914-27281687 DMR_20 774 0.507 0 EVX1
chr7:70596308-70598282 chr7:70597058-70597921 DMR_12 864 0.437 0 GALNT17
chr16:51182904-51186266 chr16:51184355-51185772 DMR_1 1418 0.422 0 SALL1
chr19:35605533-35607221 chr19:35606534-35607209 DMR_25 676 0.400 0 FXYD3
chr22:46480891-46482023 chr22:46481603-46482023 DMR_21 421 0.372 0 MIRLET7BHG
chr6:28602513-28603779 chr6:28602705-28603173 DMR_17 469 0.370 0 SCAND3
chr6:30130109-30132715 chr6:30131283-30132133 DMR_4 851 0.326 0 TRIM15

0.7 CMEnt-specific plots

dmr_rank_table <- as.data.frame(dmrs_cment)
CMEnt::plotDMRs(
    dmrs_cment,
    top_n = 1,
    score_by = "score",
    plot_motif = FALSE,
    beta = beta_handler,
    pheno = pheno,
    sample_group_col = "tissue.ch1",
    array = "EPIC",
    genome = "hg19"
)
Aligned panels show the top-scoring DMR, seed and surrounding-site positions, annotated gene context, and a heatmap of sample beta values along the locus.

Figure 4: Structure, gene context, and sample methylation of the top CMEnt DMR ranked by score in GSE244696

cat("CMEnt DMR Manhattan plot:\n")

CMEnt DMR Manhattan plot:

CMEnt::plotDMRsManhattan(
    dmrs_cment,
    genome = "hg19"
)
Manhattan plot places DMR scores along genomic position across chromosomes, with point colors indicating primary gene annotation class.

Figure 5: Genome-wide CMEnt DMR scores in GSE244696, colored by primary gene annotation class

cat("CMEnt DMR Circos interactions plot:\n")

CMEnt DMR Circos interactions plot:

CMEnt::plotDMRsCircos(
    dmrs_cment,
    genome = "hg19",
    array = "EPIC"
)
Circular chromosome tracks show DMR positions and methylation direction. Links connect interacting DMRs; the adjacent legends explain track colors and motif-supported interaction groups.

Figure 6: Circos summary of CMEnt DMR locations, methylation direction, and detected interactions in GSE244696

0.8 Biological Interpretation

annotatr is used for CpG and gene-region context, while missMethyl is applied to promoter-overlapping DMRs for array-aware GO/KEGG enrichment because it accounts for probe-per-gene bias in methylation arrays, and it was also included in the original publication.

top_10_cment_dmrs <- as.data.frame(dmrs_cment) |>
    arrange(by = -score) |>
    head(10) |>
    dplyr::select(id, score, cv_accuracy, in_promoter_of, in_gene_body_of)
knitr::kable(
    top_10_cment_dmrs,
    digits = 3,
    row.names = FALSE,
    caption = "Top 10 CMEnt DMRs with gene annotations"
)
(#tab:top_scoring_dmrs_genes)Top 10 CMEnt DMRs with gene annotations
id score cv_accuracy in_promoter_of in_gene_body_of
chr5:cg06715097-cg27197785 0.612 0.731 CAMK2A NA
chr15:cg11681321-cg21398489 0.597 0.846 ANXA2 ANXA2
chr22:cg04329382-cg07546360 0.590 0.808 MIRLET7BHG NA
chr13:cg11970797-cg04595004 0.584 0.692 NA CRYL1
chr8:cg16163338-cg09920804 0.579 0.769 PLAT NA
chr12:cg04796818-cg06992831 0.578 0.769 NA LOC283387
chr8:cg04311582-cg07188077 0.576 0.769 NA PEBP4
chr3:cg11268077-cg07552538 0.575 0.731 NA CMSS1,LOC105374010,FILIP1L
chr22:cg22887880-cg23349242 0.573 0.769 SUSD2 NA
chr1:cg17079374-cg20800117 0.571 0.808 NA CAMTA1 
normalize_missmethyl_result <- function(x, method, collection) {
    if (is.null(x) || nrow(x) == 0L) {
        return(NULL)
    }
    x <- as.data.frame(x)
    id_col <- intersect(c("ID", "PathwayID", "Pathway.ID", "KEGGID"), colnames(x))[1]
    term_col <- intersect(c("Term", "TERM", "Description", "Pathway"), colnames(x))[1]
    fdr_col <- intersect(c("FDR", "adj.P.Val", "P.DE", "P.Value"), colnames(x))[1]
    p_col <- intersect(c("P.DE", "P.Value", "P.DE.weighted"), colnames(x))[1]
    de_col <- intersect(c("DE", "N_DE", "Count"), colnames(x))[1]
    n_col <- intersect(c("N", "Total", "Size"), colnames(x))[1]
    ont_col <- intersect(c("Ont", "Ontology", "Category"), colnames(x))[1]
    if (is.na(term_col) || is.na(fdr_col)) {
        return(NULL)
    }
    pathway_id <- if (!is.na(id_col)) as.character(x[[id_col]]) else rownames(x)
    pathway_id <- ifelse(collection == "KEGG" & grepl("^(path:)?hsa[0-9]{5}$|^[0-9]{5}$", pathway_id), pathway_id, NA_character_)
    data.frame(
        Method = method,
        Collection = collection,
        PathwayID = pathway_id,
        Term = as.character(x[[term_col]]),
        Ontology = if (!is.na(ont_col)) as.character(x[[ont_col]]) else collection,
        N = if (!is.na(n_col)) suppressWarnings(as.numeric(x[[n_col]])) else NA_real_,
        DE = if (!is.na(de_col)) suppressWarnings(as.numeric(x[[de_col]])) else NA_real_,
        PValue = if (!is.na(p_col)) suppressWarnings(as.numeric(x[[p_col]])) else NA_real_,
        FDR = suppressWarnings(as.numeric(x[[fdr_col]])),
        stringsAsFactors = FALSE
    )
}
library(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
missmethyl_long <- NULL
missmethyl_array <- if (toupper(array_type) == "EPICV2") "EPIC_V2" else toupper(array_type)
missmethyl_regions <- lapply(names(methods_list), function(method) {
    gr <- methods_list[[method]]
    if (!"in_promoter_of" %in% colnames(mcols(gr))) {
        return(gr[integer()])
    }
    gr[has_text(mcols(gr)$in_promoter_of)]
})
names(missmethyl_regions) <- names(methods_list)
missmethyl_region_summary <- data.frame(
    Method = names(methods_list),
    Promoter_overlapping_DMRs = vapply(missmethyl_regions, length, integer(1)),
    Total_DMRs = vapply(methods_list, length, integer(1)),
    stringsAsFactors = FALSE
)
missmethyl_region_summary$Percentage <- 100 * missmethyl_region_summary$Promoter_overlapping_DMRs / pmax(1L, missmethyl_region_summary$Total_DMRs)
knitr::kable(
    missmethyl_region_summary,
    digits = 1,
    row.names = FALSE,
    caption = "CMEnt-annotated promoter-overlapping DMRs supplied to missMethyl"
)

Table 8: CMEnt-annotated promoter-overlapping DMRs supplied to missMethyl
Method Promoter_overlapping_DMRs Total_DMRs Percentage
CMEnt 1455 2104 69.2
BH 25 29 86.2 
cached <- getBenchmarkFile("missmethyl_cment_gene_annotation_promoter_enrichment_results.rds")
if (file.exists(cached)) {
    missmethyl_long <- readRDS(cached)
    if (!all(c("PathwayID", "RegionScope") %in% colnames(missmethyl_long))) {
        unlink(cached)
        missmethyl_long <- NULL
    }
}
if (is.null(missmethyl_long)) {
    missmethyl_results <- Filter(Negate(is.null), unlist(lapply(names(missmethyl_regions), function(method) {
        regions <- missmethyl_regions[[method]]
        if (length(regions) == 0L) {
            cat(sprintf("Skipping missMethyl enrichment for %s: no promoter-overlapping DMRs.\n", method))
            return(NULL)
        }
        lapply(c("GO", "KEGG"), function(collection) {
            cat(sprintf("Running missMethyl %s enrichment for %s CMEnt-annotated promoter-overlapping DMRs (n = %d)...\n", collection, method, length(regions)))
            ret <- tryCatch(
                missMethyl::goregion(
                    regions = regions,
                    all.cpg = beta_row_names,
                    collection = collection,
                    array.type = missmethyl_array,
                    plot.bias = FALSE,
                    prior.prob = TRUE
                ),
                error = function(e) {
                    warning(sprintf("missMethyl %s enrichment failed for %s: %s", collection, method, e$message))
                    NULL
                }
            )
            ret <- normalize_missmethyl_result(ret, method, collection)
            if (!is.null(ret)) {
                ret$RegionScope <- "CMEnt-annotated promoter-overlapping DMRs"
            }
            ret
        })
    }), recursive = FALSE))
    missmethyl_long <- if (length(missmethyl_results) > 0L) {
        base::do.call(rbind, missmethyl_results)
    } else {
        NULL
    }
    saveRDS(missmethyl_long, cached)
}
if (!is.null(missmethyl_long) && nrow(missmethyl_long) > 0L) {
    missmethyl_long <- missmethyl_long[is.finite(missmethyl_long$FDR), , drop = FALSE]
    missmethyl_long$NegLog10FDR <- -log10(pmax(missmethyl_long$FDR, .Machine$double.xmin))

    missmethyl_top <- missmethyl_long |>
        group_by(Method, Collection) |>
        arrange(FDR, .by_group = TRUE) |>
        slice_head(n = 10) |>
        ungroup()
    missmethyl_top$Term_Label <- short_label(missmethyl_top$Term, 70)
    missmethyl_top$Facet <- paste(missmethyl_top$Method, missmethyl_top$Collection, sep = " - ")
    shared_missmethyl_terms <- missmethyl_top |>
        distinct(Collection, Term, Method) |>
        dplyr::count(Collection, Term, name = "Methods") |>
        filter(Methods > 1L) |>
        transmute(Term_Key = paste(Collection, Term, sep = "___")) |>
        pull(Term_Key)
    missmethyl_top$Shared_Term <- paste(missmethyl_top$Collection, missmethyl_top$Term, sep = "___") %in% shared_missmethyl_terms
    missmethyl_top$Term_Facet <- factor(
        paste(missmethyl_top$Term_Label, missmethyl_top$Facet, sep = "___"),
        levels = rev(unique(paste(missmethyl_top$Term_Label, missmethyl_top$Facet, sep = "___")))
    )
    missmethyl_axis_term_labels <- gsub("\\\\", "\\\\\\\\", missmethyl_top$Term_Label)
    missmethyl_axis_term_labels <- gsub('"', '\\"', missmethyl_axis_term_labels, fixed = TRUE)
    missmethyl_axis_term_labels <- paste0('"', missmethyl_axis_term_labels, '"')
    missmethyl_axis_labels <- setNames(
        ifelse(missmethyl_top$Shared_Term, paste0("bold(", missmethyl_axis_term_labels, ")"), missmethyl_axis_term_labels),
        as.character(missmethyl_top$Term_Facet)
    )

    knitr::kable(
        missmethyl_top |>
            dplyr::select(Method, Collection, PathwayID, Ontology, Term, N, DE, PValue, FDR) |>
            arrange(Collection, Method, FDR),
        digits = 3,
        row.names = FALSE,
        caption = "Top missMethyl GO/KEGG enrichments from promoter-overlapping DMRs"
    )
}

Table 9: Top missMethyl GO/KEGG enrichments from promoter-overlapping DMRs
Method Collection PathwayID Ontology Term N DE PValue FDR
BH GO NA GO homophilic cell-cell adhesion 170 17 0.000 0.000
BH GO NA GO calcium ion binding 677 18 0.000 0.000
BH GO NA GO cell adhesion molecule binding 624 17 0.000 0.000
BH GO NA GO cell-cell adhesion 889 17 0.000 0.000
BH GO NA GO cell adhesion 1393 19 0.000 0.000
BH GO NA GO metal ion binding 4170 27 0.000 0.000
BH GO NA GO cation binding 4261 27 0.000 0.000
BH GO NA GO nervous system development 2087 18 0.000 0.001
BH GO NA GO multicellular organism development 3798 23 0.000 0.005
BH GO NA GO system development 3278 21 0.000 0.007
CMEnt GO NA GO DNA-binding transcription factor activity, RNA polymerase II-specific 1254 182 0.000 0.000
CMEnt GO NA GO DNA-binding transcription factor activity 1351 191 0.000 0.000
CMEnt GO NA GO developmental process 5516 596 0.000 0.000
CMEnt GO NA GO anatomical structure development 4968 551 0.000 0.000
CMEnt GO NA GO RNA polymerase II transcription regulatory region sequence-specific DNA binding 1265 178 0.000 0.000
CMEnt GO NA GO homophilic cell-cell adhesion 170 49 0.000 0.000
CMEnt GO NA GO system development 3278 395 0.000 0.000
CMEnt GO NA GO DNA-binding transcription activator activity 437 87 0.000 0.000
CMEnt GO NA GO RNA polymerase II cis-regulatory region sequence-specific DNA binding 1072 153 0.000 0.000
CMEnt GO NA GO cis-regulatory region sequence-specific DNA binding 1095 156 0.000 0.000
BH KEGG hsa04519 KEGG Cadherin signaling 394 17 0.000 0.000
BH KEGG hsa00010 KEGG Glycolysis / Gluconeogenesis 64 0 1.000 1.000
BH KEGG hsa00020 KEGG Citrate cycle (TCA cycle) 28 0 1.000 1.000
BH KEGG hsa00030 KEGG Pentose phosphate pathway 26 0 1.000 1.000
BH KEGG hsa00040 KEGG Pentose and glucuronate interconversions 36 0 1.000 1.000
BH KEGG hsa00051 KEGG Fructose and mannose metabolism 35 0 1.000 1.000
BH KEGG hsa00052 KEGG Galactose metabolism 30 0 1.000 1.000
BH KEGG hsa00053 KEGG Ascorbate and aldarate metabolism 29 0 1.000 1.000
BH KEGG hsa00061 KEGG Fatty acid biosynthesis 16 0 1.000 1.000
BH KEGG hsa00062 KEGG Fatty acid elongation 28 0 1.000 1.000
CMEnt KEGG hsa04519 KEGG Cadherin signaling 394 67 0.000 0.000
CMEnt KEGG hsa00053 KEGG Ascorbate and aldarate metabolism 29 6 0.002 0.294
CMEnt KEGG hsa04382 KEGG Cornified envelope formation 206 25 0.002 0.294
CMEnt KEGG hsa00920 KEGG Sulfur metabolism 10 4 0.003 0.303
CMEnt KEGG hsa00140 KEGG Steroid hormone biosynthesis 63 9 0.008 0.501
CMEnt KEGG hsa00980 KEGG Metabolism of xenobiotics by cytochrome P450 71 9 0.008 0.501
CMEnt KEGG hsa04950 KEGG Maturity onset diabetes of the young 26 6 0.011 0.591
CMEnt KEGG hsa00565 KEGG Ether lipid metabolism 51 8 0.028 0.758
CMEnt KEGG hsa00591 KEGG Linoleic acid metabolism 30 5 0.033 0.758
CMEnt KEGG hsa00592 KEGG alpha-Linolenic acid metabolism 26 5 0.031 0.758 
if (!is.null(missmethyl_long) && nrow(missmethyl_long) > 0L) {
    ggplot(missmethyl_top, aes(x = NegLog10FDR, y = Term_Facet, size = DE)) +
        geom_point(alpha = 0.9) +
        facet_wrap(~Facet, scales = "free_y", ncol = 2) +
        scale_y_discrete(labels = function(x) parse(text = unname(missmethyl_axis_labels[x]))) +
        theme_minimal(base_size = 17) +
        labs(
            x = expression(-log[10]("FDR")), y = "", color = "", size = "DM genes"
        )
} else {
    cat("missMethyl returned no enrichment results.\n")
}
Faceted dot plots show enriched terms against negative log10 FDR by method and collection. Dot size gives the number of differentially methylated genes, and bold labels identify terms shared across methods.

Figure 7: Top missMethyl GO and KEGG enrichments from promoter-overlapping DMRs by method
Bold term labels indicate terms shared across methods. 

if (!is.null(missmethyl_long) && nrow(missmethyl_long) > 0L &&
    requireNamespace("pathview", quietly = TRUE) &&
    requireNamespace("KEGGREST", quietly = TRUE)) {
    library(pathview)
    kegg_candidates <- missmethyl_long |>
        filter(Method == "CMEnt", Collection == "KEGG", !is.na(PathwayID), is.finite(PValue)) |>
        arrange(FDR, PValue) |>
        mutate(
            PathwayID = sub("^path:", "", PathwayID),
            PathwayID = ifelse(grepl("^hsa", PathwayID), PathwayID, paste0("hsa", PathwayID))
        )
    promoter_gene_delta <- as.data.frame(dmrs_cment) |>
        arrange(desc(score)) |>
        filter(is.finite(as.numeric(delta_beta_promoter)), !is.na(in_promoter_of), nzchar(in_promoter_of)) |>
        mutate(delta_beta_promoter = as.numeric(delta_beta_promoter)) |>
        dplyr::select(id, score, delta_beta_promoter, in_promoter_of) |>
        tidyr::separate_rows(in_promoter_of, sep = "\\s*,\\s*") |>
        mutate(Gene = trimws(in_promoter_of)) |>
        filter(nzchar(Gene)) |>
        arrange(desc(abs(delta_beta_promoter))) |>
        distinct(Gene, .keep_all = TRUE)

    if (nrow(kegg_candidates) > 0L && nrow(promoter_gene_delta) > 0L) {
        gene_map <- suppressMessages(AnnotationDbi::select(
            org.Hs.eg.db::org.Hs.eg.db,
            keys = promoter_gene_delta$Gene,
            keytype = "SYMBOL",
            columns = c("SYMBOL", "ENTREZID")
        )) |>
            filter(!is.na(ENTREZID)) |>
            distinct(SYMBOL, ENTREZID)

        promoter_gene_delta <- promoter_gene_delta |>
            inner_join(gene_map, by = c("Gene" = "SYMBOL")) |>
            arrange(desc(abs(delta_beta_promoter))) |>
            distinct(ENTREZID, .keep_all = TRUE)

        if (nrow(promoter_gene_delta) == 0L) {
            cat("No promoter genes could be mapped to Entrez IDs for pathview.\n")
        } else {
            kegg_gene_map_file <- getBenchmarkFile("kegg_hsa_pathway_gene_map.rds")
            kegg_gene_map <- if (file.exists(kegg_gene_map_file)) {
                readRDS(kegg_gene_map_file)
            } else {
                links <- tryCatch(
                    KEGGREST::keggLink("hsa", "pathway"),
                    error = function(e) {
                        warning(sprintf("Could not retrieve KEGG pathway-gene links: %s", e$message))
                        NULL
                    }
                )
                if (is.null(links)) {
                    data.frame(PathwayID = character(), ENTREZID = character())
                } else {
                    data.frame(
                        PathwayID = sub("^path:", "", names(links)),
                        ENTREZID = sub("^hsa:", "", unname(links)),
                        stringsAsFactors = FALSE
                    ) |>
                        distinct()
                }
            }
            if (!file.exists(kegg_gene_map_file)) {
                saveRDS(kegg_gene_map, kegg_gene_map_file)
            }

            candidate_overlaps <- kegg_candidates |>
                left_join(kegg_gene_map, by = "PathwayID") |>
                filter(ENTREZID %in% promoter_gene_delta$ENTREZID) |>
                group_by(PathwayID) |>
                summarise(OverlapEntrez = list(unique(ENTREZID)), .groups = "drop")

            matched_pathway <- match(kegg_candidates$PathwayID, candidate_overlaps$PathwayID)
            matched_pathway <- which(!is.na(matched_pathway))[1]

            if (is.na(matched_pathway)) {
                cat("No missMethyl KEGG pathway contained the top promoter-overlapping DMR genes for pathview coloring.\n")
            } else {
                top_kegg <- kegg_candidates |> dplyr::slice(matched_pathway)
                matched_genes <- candidate_overlaps$OverlapEntrez[[match(top_kegg$PathwayID[1], candidate_overlaps$PathwayID)]]
                promoter_pathway_delta <- promoter_gene_delta |>
                    filter(ENTREZID %in% matched_genes)
                gene_delta <- promoter_pathway_delta$delta_beta_promoter
                names(gene_delta) <- promoter_pathway_delta$ENTREZID
                kegg_id <- sub("^hsa", "", top_kegg$PathwayID[1])
                pathview_dir <- getBenchmarkFile("pathview")
                dir.create(pathview_dir, showWarnings = FALSE, recursive = TRUE)
                pathview_suffix <- "promoter_delta_beta"
                pathview_png <- file.path(pathview_dir, paste0("hsa", kegg_id, ".", pathview_suffix, ".png"))
                gene_limit <- max(0.001, max(abs(gene_delta), na.rm = TRUE))
                unlink(c(pathview_png, basename(pathview_png)))

                pathview_result <- tryCatch(
                    pathview::pathview(
                        gene.data = gene_delta,
                        gene.idtype = "entrez",
                        pathway.id = kegg_id,
                        species = "hsa",
                        out.suffix = pathview_suffix,
                        kegg.dir = pathview_dir,
                        limit = list(gene = gene_limit, cpd = 1),
                        low = list(gene = "#2F80ED", cpd = "#2F80ED"),
                        mid = list(gene = "#FFFFFF", cpd = "#FFFFFF"),
                        high = list(gene = "#FF3B30", cpd = "#FF3B30")
                    ),
                    error = function(e) {
                        warning(sprintf("pathview failed for KEGG pathway %s: %s", kegg_id, e$message))
                        NULL
                    }
                )
                if (file.exists(basename(pathview_png))) {
                    file.rename(basename(pathview_png), pathview_png)
                }
                pathview_mapped <- if (!is.null(pathview_result$plot.data.gene)) {
                    sum(!is.na(pathview_result$plot.data.gene[["mol.data"]]))
                } else {
                    NA_integer_
                }

                knitr::kable(
                    promoter_pathway_delta |>
                        dplyr::select(Gene, ENTREZID, delta_beta_promoter, score, id) |>
                        arrange(desc(abs(delta_beta_promoter))) |>
                        head(20),
                    digits = 3,
                    row.names = FALSE,
                    caption = paste0(
                        "Promoter delta beta values supplied to pathview for the top missMethyl KEGG pathway with promoter-gene overlap: ",
                        top_kegg$Term[1],
                        ". Genes overlapping this KEGG pathway: ", nrow(promoter_pathway_delta),
                        ". Genes mapped by pathview: ", pathview_mapped,
                        ". Red indicates higher methylation in Tumor; blue indicates lower methylation in Tumor."
                    )
                )

                if (file.exists(pathview_png)) {
                    knitr::include_graphics(pathview_png)
                } else {
                    cat("pathview did not produce a KEGG pathway image.\n")
                }
            }
        }
    } else {
        cat("No CMEnt KEGG pathway with ID or promoter gene delta-beta values were available for pathview.\n")
    }
} else {
    cat("pathview and KEGGREST are required for the KEGG pathway diagram.\n")
}
KEGG pathway diagram colors mapped genes by promoter delta beta: red for higher methylation in Tumor, blue for lower methylation in Tumor, and white near zero.

Figure 8: Top missMethyl KEGG pathway with promoter-gene overlap, colored by CMEnt promoter delta beta in GSE244696

0.9 Conclusions

In the GSE244696 PDAC tumor-stroma comparison, CMEnt was applied in meta-analysis mode by using the published DMPs as external evidence and the available beta matrix to assemble spatially coherent DMRs around those signals. The resulting DMRs recovered biological themes that are consistent with the original study and with PDAC tumor microenvironment biology. Study-context terms highlighted extracellular matrix, myeloid leukocyte differentiation, steroid-hormone response and hormone-regulatory processes. These signals are plausible because PDAC is characterized by a dense desmoplastic extracellular matrix and because the original study reported compartment-specific immune and estrogen receptor-related epigenetic dysregulation, with tumor-associated macrophages enriched in tumor-cell regions and associated with patient outcome. The annotation profile further supports a regulatory interpretation: CMEnt DMRs overlap promoters, exons, introns and 5′ UTRs, indicating that the method organizes CpG-level evidence into regions connected to gene-regulatory architecture rather than producing isolated CpG hits. The genome-wide circos view shows both hypermethylated and hypomethylated DMRs, consistent with a compartment-specific remodeling process. A motif-supported DMR interaction involving AHR::ARNT provides an additional hypothesis-generating link to xenobiotic and cytochrome P450-related metabolism, although this motif result should not be overinterpreted given the small number of motif-supported interactions. Overall, these results suggest that CMEnt preserves the biological context of the published DMP evidence while emphasizing coherent regulatory, extracellular-matrix, myeloid and hormone-associated tumor-stroma programs.

0.10 Session Info

sessionInfo()

R version 4.6.1 (2026-06-24) Platform: x86_64-pc-linux-gnu Running under: Ubuntu 24.04.4 LTS

Matrix products: default BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0 LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 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: Europe/Brussels tzcode source: system (glibc)

attached base packages: [1] parallel stats4 stats graphics grDevices datasets utils
[8] methods base

other attached packages: [1] pathview_1.52.0
[2] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0 [3] scales_1.4.0
[4] minfi_1.58.0
[5] bumphunter_1.54.0
[6] locfit_1.5-9.12
[7] iterators_1.0.14
[8] foreach_1.5.2
[9] Biostrings_2.80.2
[10] XVector_0.52.0
[11] CMEnt_0.99.14
[12] DSS_2.60.0
[13] BiocParallel_1.46.0
[14] bsseq_1.48.0
[15] SummarizedExperiment_1.42.0
[16] Biobase_2.72.0
[17] MatrixGenerics_1.24.0
[18] matrixStats_1.5.0
[19] ggplot2_4.0.3
[20] GenomicRanges_1.64.0
[21] Seqinfo_1.2.0
[22] IRanges_2.46.0
[23] S4Vectors_0.50.3
[24] BiocGenerics_0.58.1
[25] generics_0.1.4
[26] dplyr_1.2.1
[27] data.table_1.18.6.1
[28] testthat_3.3.2
[29] BiocStyle_2.40.0

loaded via a namespace (and not attached): [1] graph_1.90.0
[2] igraph_2.3.4
[3] plotly_4.12.1
[4] Formula_1.2-6
[5] devtools_2.5.2
[6] tidyselect_1.2.1
[7] bit_4.6.0
[8] doParallel_1.0.17
[9] clue_0.3-68
[10] lattice_0.23-1
[11] rjson_0.2.23
[12] nor1mix_1.3-3
[13] blob_1.3.0
[14] stringr_1.6.0
[15] rngtools_1.5.2
[16] S4Arrays_1.12.1
[17] base64_2.0.2
[18] dichromat_2.0-1
[19] scrime_1.3.7
[20] seqLogo_1.78.0
[21] png_0.1-9
[22] tinytex_0.61
[23] cli_3.6.6
[24] ggplotify_0.1.4
[25] bedr_1.1.5
[26] ProtGenerics_1.44.0
[27] askpass_1.2.1
[28] outliers_0.15
[29] multtest_2.68.0
[30] openssl_2.4.2
[31] pkgdown_2.2.1
[32] textshaping_1.0.5
[33] BiocIO_1.22.0
[34] purrr_1.2.2
[35] curl_8.1.0
[36] mime_0.13
[37] evaluate_1.0.5
[38] tidytree_0.4.8
[39] ComplexHeatmap_2.28.0
[40] stringi_1.8.9
[41] backports_1.5.1
[42] desc_1.4.3
[43] XML_3.99-0.25
[44] httpuv_1.6.17
[45] AnnotationDbi_1.74.0
[46] magrittr_2.0.5
[47] clusterProfiler_4.20.2
[48] rappdirs_0.3.4
[49] splines_4.6.1
[50] mclust_6.1.3
[51] DelayedDataFrame_1.28.0
[52] jpeg_0.1-11
[53] doRNG_1.8.6.3
[54] ggraph_2.2.2
[55] rentrez_1.2.4
[56] org.Hs.eg.db_3.23.1
[57] DT_0.34.0
[58] sessioninfo_1.2.4
[59] DBI_1.3.0
[60] HDF5Array_1.40.0
[61] reactome.db_1.96.0
[62] jquerylib_0.1.4
[63] genefilter_1.94.0
[64] withr_3.0.3
[65] class_7.3-24
[66] systemfonts_1.3.2
[67] rprojroot_2.1.1
[68] enrichplot_1.32.1
[69] ggnewscale_0.5.2
[70] benchmarkme_1.0.8
[71] brio_1.1.5
[72] tidygraph_1.3.1
[73] formatR_1.14
[74] rtracklayer_1.72.0
[75] BiocManager_1.30.27
[76] htmlwidgets_1.6.4
[77] fs_2.1.0
[78] biomaRt_2.68.0
[79] missMethyl_1.46.0
[80] ggrepel_0.9.8
[81] labeling_0.4.3
[82] SparseArray_1.12.3
[83] cellranger_1.1.0
[84] shinycssloaders_1.1.0
[85] h5mread_1.4.1
[86] annotate_1.90.0
[87] VariantAnnotation_1.58.0
[88] knitr_1.52
[89] TFBSTools_1.50.0
[90] UCSC.utils_1.8.0
[91] beanplot_1.3.1
[92] TFMPvalue_1.0.0
[93] annotatr_1.38.3
[94] patchwork_1.3.2
[95] caTools_1.18.4
[96] grid_4.6.1
[97] ggtree_4.2.0
[98] rhdf5_2.56.1
[99] pwalign_1.8.0
[100] R.oo_1.27.1
[101] ggiraph_0.9.6
[102] regioneR_1.44.0
[103] gridGraphics_0.5-1
[104] ellipsis_0.3.3
[105] lazyeval_0.2.3
[106] yaml_2.3.12
[107] survival_3.8-12
[108] BiocVersion_3.23.1
[109] crayon_1.5.3
[110] Rgraphviz_2.56.0
[111] tweenr_2.0.3
[112] RColorBrewer_1.1-3
[113] tidyr_1.3.2
[114] later_1.4.8
[115] codetools_0.2-20
[116] base64enc_0.1-6
[117] GlobalOptions_0.1.4
[118] KEGGREST_1.52.2
[119] shape_1.4.6.1
[120] ReactomePA_1.56.0
[121] limma_3.68.5
[122] gdtools_0.5.2
[123] Rsamtools_2.28.0
[124] filelock_1.0.3
[125] showtext_0.9-8
[126] foreign_0.8-91
[127] pkgconfig_2.0.3
[128] KEGGgraph_1.72.0
[129] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.1 [130] xml2_1.6.0
[131] GenomicAlignments_1.48.0
[132] aplot_0.3.2
[133] hummingbird_1.22.0
[134] BSgenome_1.80.0
[135] ape_5.8-1
[136] viridisLite_0.4.3
[137] biovizBase_1.60.0
[138] TxDb.Hsapiens.UCSC.hg19.knownGene_3.22.1
[139] gridBase_0.4-7
[140] xtable_1.8-8
[141] interp_1.1-6
[142] plyr_1.8.9
[143] httr_1.4.9
[144] tools_4.6.1
[145] pkgbuild_1.4.8
[146] htmlTable_2.5.0
[147] checkmate_2.3.4
[148] nlme_3.1-170
[149] futile.logger_1.4.9
[150] lambda.r_1.2.6
[151] dbplyr_2.6.0
[152] ExperimentHub_3.2.2
[153] digest_0.6.39
[154] permute_0.9-10
[155] bookdown_0.48
[156] Matrix_1.7-6
[157] farver_2.1.2
[158] tzdb_0.5.0
[159] AnnotationFilter_1.36.0
[160] reshape2_1.4.5
[161] viridis_0.6.5
[162] yulab.utils_0.2.5
[163] DirichletMultinomial_1.54.1
[164] rpart_4.1.27
[165] glue_1.8.1
[166] cachem_1.1.0
[167] bspm_0.5.8
[168] VennDiagram_1.8.2
[169] BiocFileCache_3.2.0
[170] polyclip_1.10-7
[171] Hmisc_5.3-0
[172] sysfonts_0.8.9
[173] pkgload_1.5.3
[174] statmod_1.5.2
[175] fontBitstreamVera_0.1.1
[176] GEOquery_2.80.0
[177] httr2_1.3.0
[178] showtextdb_3.0
[179] gson_0.2.2
[180] dmrseq_1.32.0
[181] graphlayouts_1.2.5
[182] siggenes_1.86.0
[183] gtools_3.9.5
[184] readxl_1.6.0
[185] preprocessCore_1.74.0
[186] gridExtra_2.3.1
[187] shiny_1.14.0
[188] tidydr_0.0.7
[189] R.utils_2.13.0
[190] rhdf5filters_1.24.1
[191] RCurl_1.98-1.20
[192] memoise_2.0.1
[193] rmarkdown_2.32
[194] R.methodsS3_1.8.2
[195] svglite_2.2.2
[196] reshape_0.8.10
[197] fontLiberation_0.1.0
[198] illuminaio_0.54.0
[199] rstudioapi_0.19.0
[200] cluster_2.1.8.3
[201] hms_1.1.4
[202] DMRcate_3.8.0
[203] colorspace_2.1-4
[204] FNN_1.1.4.1
[205] rlang_1.3.0
[206] quadprog_1.5-8
[207] GenomeInfoDb_1.48.0
[208] DelayedMatrixStats_1.34.0
[209] sparseMatrixStats_1.24.0
[210] ggforce_0.5.0
[211] circlize_0.4.18
[212] ggtangle_0.1.4
[213] xfun_0.62
[214] ggseqlogo_0.2.2
[215] e1071_1.7-17
[216] aisdk_1.4.12
[217] abind_1.4-8
[218] GOSemSim_2.38.3
[219] tibble_3.3.1
[220] treeio_1.36.1
[221] Rhdf5lib_2.0.0
[222] readr_2.2.0
[223] futile.options_1.0.1
[224] bitops_1.1-0
[225] ps_1.9.3
[226] promises_1.5.0
[227] qvalue_2.44.0
[228] scatterpie_0.2.6
[229] inline_0.3.21
[230] RSQLite_3.53.3
[231] DelayedArray_0.38.2
[232] proxy_0.4-29
[233] GO.db_3.23.1
[234] compiler_4.6.1
[235] prettyunits_1.2.0
[236] beachmat_2.28.0
[237] graphite_1.58.0
[238] Rcpp_1.1.2
[239] benchmarkmeData_2.0.0
[240] fontquiver_0.2.1
[241] edgeR_4.10.6
[242] AnnotationHub_4.2.2
[243] usethis_3.2.2
[244] MASS_7.3-66
[245] progress_1.2.3
[246] R6_2.6.1
[247] fastmap_1.2.0
[248] ensembldb_2.36.1
[249] enrichit_0.2.6
[250] nnet_7.3-21
[251] gtable_0.3.6
[252] strex_2.1.0
[253] latticeExtra_0.6-31
[254] deldir_2.0-4
[255] htmltools_0.5.9
[256] bit64_4.8.6
[257] lifecycle_1.0.5
[258] S7_0.2.2
[259] processx_3.9.0
[260] callr_3.8.0
[261] Gviz_1.56.0
[262] restfulr_0.0.17
[263] sass_0.4.10
[264] vctrs_0.7.3
[265] DOSE_4.6.0
[266] ggfun_0.2.2
[267] bslib_0.12.0
[268] pillar_1.11.1
[269] GenomicFeatures_1.64.0
[270] magick_2.9.1
[271] otel_0.2.0
[272] jsonlite_2.0.0
[273] cigarillo_1.2.1
[274] GetoptLong_1.1.1