CMEnt Meta-Analysis on GSE244696

Building DMRs from the published DMP list and comparing with the supplied Bumphunter regions

CMEnt Package

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. The GEO deposition does not include the complete raw IDAT set used by the original study, so a head-to-head DMR caller benchmark on the locally available beta matrix would be unfair. This notebook therefore 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 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)
    }
}

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

try(devtools::load_all(), silent = TRUE)
suppressPackageStartupMessages({
    try(library(CMEnt), silent = TRUE)
    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)
}

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 Local 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$end,
        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(
        "Local samples",
        "Local Tumor samples",
        "Local Stroma samples",
        "Published DMPs",
        "Published DMPs present in local 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
Local samples	26
Local Tumor samples	13
Local Stroma samples	13
Published DMPs	25410
Published DMPs present in local 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 in the locally available samples")

Table 2: Direction of published DMPs in the locally available samples

Direction	Count
Higher in Tumor	17566
Higher in Stroma	5101

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") +
    theme_meta() +
    labs(
        title = "Local Effect Size for Published DMPs",
        subtitle = "Tumor minus Stroma beta difference, computed only from locally available samples",
        x = "Delta beta in local samples", y = "Published 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 535 seconds.

if (is.null(dmrs_cment)) {
    dmrs_cment <- GRanges()
}
methods_list <- list(CMEnt = dmrs_cment, `Published Bumphunter` = published_bumphunter)

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
Published Bumphunter	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
score	0.361	0.478	0.479	0.627
cv_accuracy	0.115	0.500	0.509	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)")

0.5 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
Published Bumphunter	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 (%)")

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.6 CMEnt-specific plots

dmr_rank_table <- as.data.frame(dmrs_cment)
top_dmr_idx <- if ("qval" %in% colnames(dmr_rank_table) && any(is.finite(dmr_rank_table$qval))) {
    which.min(dmr_rank_table$qval)
} else {
    delta_col <- intersect(c("median_delta_beta", "mean_delta_beta", "delta_beta"), colnames(dmr_rank_table))[1]
    if (is.na(delta_col)) 1L else order(abs(dmr_rank_table[[delta_col]]), decreasing = TRUE)[1]
}
cat("Top CMEnt DMR prioritized by DMR q-value when available; SVM score is complementary:\n")

Top CMEnt DMR prioritized by DMR q-value when available; SVM score is complementary:

CMEnt::plotDMRs(
    dmrs_cment,
    dmr_indices = top_dmr_idx,
    beta = beta_handler,
    pheno = pheno,
    sample_group_col = "tissue.ch1",
    array = "EPIC",
    genome = "hg19"
)

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

CMEnt DMR Manhattan plot:

CMEnt::plotDMRsManhattan(
    dmrs_cment,
    genome = "hg19"
)

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

CMEnt DMR Circos interactions plot:

CMEnt::plotDMRsCircos(
    dmrs_cment,
    genome = "hg19",
    array = "EPIC"
)

## Biological Interpretation

annotatr is used for CpG and gene-region context, while missMethyl is used 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.

context_summary <- NULL
annotatr_gene_ids_by_method <- list()
annotatr_types <- c(
    paste0(genome, "_cpg_islands"),
    paste0(genome, "_cpg_shores"),
    paste0(genome, "_cpg_shelves"),
    paste0(genome, "_cpg_inter"),
    paste0(genome, "_genes_promoters"),
    paste0(genome, "_genes_exons"),
    paste0(genome, "_genes_introns"),
    paste0(genome, "_genes_3UTRs"),
    paste0(genome, "_genes_5UTRs"),
    paste0(genome, "_genes_intergenic")
)
annotatr_types <- annotatr_types[annotatr_types %in% annotatr::builtin_annotations()]
annotatr_ann <- annotatr::build_annotations(genome = genome, annotations = annotatr_types)

annotation_type_labels <- c(
    setNames(
        c("CpG island", "CpG shore", "CpG shelf", "Open sea"),
        paste0(genome, c("_cpg_islands", "_cpg_shores", "_cpg_shelves", "_cpg_inter"))
    ),
    setNames(
        c("Promoter", "Exon", "Intron", "3' UTR", "5' UTR", "Intergenic"),
        paste0(genome, c("_genes_promoters", "_genes_exons", "_genes_introns", "_genes_3UTRs", "_genes_5UTRs", "_genes_intergenic"))
    )
)

context_summary <- base::do.call(rbind, lapply(names(methods_list), function(method) {
    gr <- methods_list[[method]]
    if (length(gr) == 0L) {
        return(NULL)
    }
    gr <- granges(gr)
    mcols(gr)$dmr_id <- seq_along(gr)
    names(gr) <- NULL
    ann_gr <- annotatr::annotate_regions(regions = gr, annotations = annotatr_ann, quiet = TRUE)
    ann_mcols <- as.data.frame(S4Vectors::mcols(ann_gr), optional = TRUE)
    rownames(ann_mcols) <- NULL
    ann_df <- data.frame(
        seqnames = as.character(seqnames(ann_gr)),
        start = start(ann_gr),
        end = end(ann_gr),
        ann_mcols,
        check.names = FALSE,
        stringsAsFactors = FALSE
    )
    if ("annot.gene_id" %in% colnames(ann_df)) {
        annotatr_gene_ids_by_method[[method]] <<- unique(stats::na.omit(as.character(ann_df$annot.gene_id)))
    }
    ann_df <- ann_df[ann_df$annot.type %in% names(annotation_type_labels), , drop = FALSE]
    if (nrow(ann_df) == 0L) {
        return(NULL)
    }
    ann_df <- unique(ann_df[, c("dmr_id", "annot.type"), drop = FALSE])
    counts <- as.data.frame(table(ann_df$annot.type), stringsAsFactors = FALSE)
    colnames(counts) <- c("annot.type", "DMRs")
    counts$Method <- method
    counts$Annotation <- unname(annotation_type_labels[as.character(counts$annot.type)])
    counts$Class <- ifelse(grepl("_cpg_", counts$annot.type), "CpG context", "Gene context")
    counts$Percentage <- 100 * counts$DMRs / length(gr)
    counts
}))
if (!is.null(context_summary) && nrow(context_summary) > 0L) {
    context_summary$Annotation <- factor(
        context_summary$Annotation,
        levels = c("CpG island", "CpG shore", "CpG shelf", "Open sea", "Promoter", "Exon", "Intron", "5' UTR", "3' UTR", "Intergenic")
    )
    knitr::kable(
        context_summary |>
            dplyr::select(Method, Class, Annotation, DMRs, Percentage) |>
            arrange(Class, Annotation, Method),
        digits = 1,
        row.names = FALSE,
        caption = "Region context summary. Categories can overlap."
    )
}

Table 8: Region context summary
Categories can overlap.

Method	Class	Annotation	DMRs	Percentage
CMEnt	CpG context	CpG island	971	46.2
Published Bumphunter	CpG context	CpG island	18	62.1
CMEnt	CpG context	CpG shore	1045	49.7
Published Bumphunter	CpG context	CpG shore	13	44.8
CMEnt	CpG context	CpG shelf	194	9.2
Published Bumphunter	CpG context	CpG shelf	2	6.9
CMEnt	CpG context	Open sea	900	42.8
Published Bumphunter	CpG context	Open sea	6	20.7
CMEnt	Gene context	Promoter	1588	75.5
Published Bumphunter	Gene context	Promoter	26	89.7
CMEnt	Gene context	Exon	1558	74.0
Published Bumphunter	Gene context	Exon	26	89.7
CMEnt	Gene context	Intron	1807	85.9
Published Bumphunter	Gene context	Intron	27	93.1
CMEnt	Gene context	5’ UTR	919	43.7
Published Bumphunter	Gene context	5’ UTR	19	65.5
CMEnt	Gene context	3’ UTR	244	11.6
Published Bumphunter	Gene context	3’ UTR	1	3.4
CMEnt	Gene context	Intergenic	91	4.3

if (!is.null(context_summary) && nrow(context_summary) > 0L) {
    ggplot(context_summary, aes(x = Annotation, y = Percentage, fill = Method)) +
        geom_col(position = position_dodge(width = 0.8), width = 0.7) +
        facet_wrap(~Class, scales = "free_x") +
        scale_fill_manual(values = method_palette, drop = FALSE) +
        theme_meta() +
        labs(title = "DMR Annotation Context", x = "", y = "DMRs overlapping annotation (%)", fill = "")
}

pdac_term_pattern <- paste(
    c(
        "immune", "immun", "macrophage", "myeloid", "leukocyte", "lymphocyte",
        "cytokine", "chemokine", "inflamm", "antigen", "interferon",
        "estrogen", "oestrogen", "hormone", "stromal", "extracellular matrix",
        "collagen", "fibroblast", "pancrea"
    ),
    collapse = "|"
)

normalize_missmethyl_result <- function(x, method, collection) {
    if (is.null(x) || nrow(x) == 0L) {
        return(NULL)
    }
    x <- as.data.frame(x)
    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)
    }
    data.frame(
        Method = method,
        Collection = collection,
        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)
cached <- getBenchmarkFile("missmethyl_enrichment_results.rds")
if (file.exists(cached)) {
    missmethyl_long <- readRDS(cached)
} else {
    missmethyl_long <- base::do.call(rbind, Filter(Negate(is.null), unlist(lapply(names(methods_list), function(method) {
        lapply(c("GO", "KEGG"), function(collection) {
            ret <- tryCatch(
                missMethyl::goregion(
                    regions = methods_list[[method]],
                    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
                }
            )
            normalize_missmethyl_result(ret, method, collection)
        })
    }), recursive = FALSE)))
    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_long$PDAC_Context <- ifelse(grepl(pdac_term_pattern, missmethyl_long$Term, ignore.case = TRUE), "Study-context term", "Other top term")

    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 = " - ")
    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 = "___")))
    )

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

Table 9: Top missMethyl GO/KEGG enrichments

Method	Collection	Ontology	Term	N	DE	PValue	FDR
CMEnt	GO	GO	DNA-binding transcription factor activity, RNA polymerase II-specific	1293	212	0	0.000
CMEnt	GO	GO	DNA-binding transcription factor activity	1390	223	0	0.000
CMEnt	GO	GO	RNA polymerase II transcription regulatory region sequence-specific DNA binding	1313	207	0	0.000
CMEnt	GO	GO	anatomical structure development	5163	642	0	0.000
CMEnt	GO	GO	developmental process	5743	690	0	0.000
CMEnt	GO	GO	system development	3423	472	0	0.000
CMEnt	GO	GO	DNA-binding transcription activator activity	453	103	0	0.000
CMEnt	GO	GO	sequence-specific double-stranded DNA binding	1523	225	0	0.000
CMEnt	GO	GO	multicellular organism development	3956	520	0	0.000
CMEnt	GO	GO	RNA polymerase II cis-regulatory region sequence-specific DNA binding	1112	178	0	0.000
Published Bumphunter	GO	GO	homophilic cell-cell adhesion	176	17	0	0.000
Published Bumphunter	GO	GO	cell adhesion molecule binding	640	18	0	0.000
Published Bumphunter	GO	GO	calcium ion binding	695	18	0	0.000
Published Bumphunter	GO	GO	cell-cell adhesion	919	17	0	0.000
Published Bumphunter	GO	GO	cell adhesion	1442	19	0	0.000
Published Bumphunter	GO	GO	metal ion binding	4326	28	0	0.000
Published Bumphunter	GO	GO	cation binding	4425	28	0	0.000
Published Bumphunter	GO	GO	nervous system development	2190	19	0	0.001
Published Bumphunter	GO	GO	multicellular organism development	3956	24	0	0.003
Published Bumphunter	GO	GO	system development	3423	22	0	0.004

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) sub("___.*$", "", x)) +
        theme_meta() +
        theme(axis.text.y = element_text(size = 7)) +
        labs(title = "missMethyl Enrichment", x = expression(-log[10]("FDR")), y = "", color = "", size = "DM genes")
} else {
    cat("missMethyl returned no enrichment results.\n")
}

focused_terms <- missmethyl_long |>
    filter(grepl(pdac_term_pattern, Term, ignore.case = TRUE)) |>
    group_by(Method, Collection) |>
    arrange(FDR, .by_group = TRUE) |>
    slice_head(n = 8) |>
    ungroup()

if (nrow(focused_terms) > 0L) {
    focused_terms$Term_Label <- short_label(focused_terms$Term, 62)
    focused_terms$Term_Method <- factor(
        paste(focused_terms$Term_Label, focused_terms$Method, focused_terms$Collection, sep = "___"),
        levels = rev(unique(paste(focused_terms$Term_Label, focused_terms$Method, focused_terms$Collection, sep = "___")))
    )
    knitr::kable(
        focused_terms |>
            dplyr::select(Method, Collection, Ontology, Term, N, DE, PValue, FDR) |>
            arrange(Collection, Method, FDR),
        digits = 3,
        row.names = FALSE,
        caption = "Top immune, macrophage, stromal, pancreatic, and hormone/estrogen-related terms"
    )
}

Table 10: Top immune, macrophage, stromal, pancreatic, and hormone/estrogen-related terms

Method	Collection	Ontology	Term	N	DE	PValue	FDR
CMEnt	GO	GO	extracellular matrix	507	69	0.002	0.299
CMEnt	GO	GO	regulation of thyroid hormone generation	8	4	0.004	0.425
CMEnt	GO	GO	thyroid hormone metabolic process	25	7	0.005	0.478
CMEnt	GO	GO	regulation of hormone levels	481	59	0.006	0.517
CMEnt	GO	GO	myeloid leukocyte differentiation	198	30	0.006	0.523
CMEnt	GO	GO	response to steroid hormone	315	43	0.009	0.674
CMEnt	GO	GO	regulation of myeloid leukocyte differentiation	107	17	0.015	0.872
CMEnt	GO	GO	negative regulation of myeloid leukocyte differentiation	46	9	0.016	0.895
Published Bumphunter	GO	GO	mitotic cytokinesis	103	0	1.000	1.000
Published Bumphunter	GO	GO	cytokinesis	190	0	1.000	1.000
Published Bumphunter	GO	GO	assembly of actomyosin apparatus involved in cytokinesis	9	0	1.000	1.000
Published Bumphunter	GO	GO	pancreatic polypeptide receptor activity	4	0	1.000	1.000
Published Bumphunter	GO	GO	growth hormone secretagogue receptor activity	1	0	1.000	1.000
Published Bumphunter	GO	GO	establishment of lymphocyte polarity	13	0	1.000	1.000
Published Bumphunter	GO	GO	immunological synapse formation	10	0	1.000	1.000
Published Bumphunter	GO	GO	immunological synapse	46	0	1.000	1.000

if (nrow(focused_terms) > 0L) {
    ggplot(focused_terms, aes(x = -log10(pmax(FDR, .Machine$double.xmin)), y = Term_Method, fill = Method)) +
        geom_col(width = 0.7) +
        facet_wrap(~Collection, scales = "free_y") +
        scale_y_discrete(labels = function(x) sub("___.*$", "", x)) +
        scale_fill_manual(values = method_palette, drop = FALSE) +
        theme_meta() +
        theme(axis.text.y = element_text(size = 8)) +
        labs(title = "Study-Relevant Enrichment Terms", x = expression(-log[10]("FDR")), y = "", fill = "")
} else {
    cat("No missMethyl terms matched the study-focused keyword set.\n")
}

0.7 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.8 Session Info

sessionInfo()

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