Contents

0.1 Introduction

This vignette analyzes GSE63409, an Illumina 450K methylation dataset from acute myeloid leukemia (AML). The dataset contains 15 AML patients represented by leukemia stem and blast cell samples, and normal hematopoietic stem and progenitor cell controls from healthy bone marrow.

The original study, Hu et al. 2019 (doi:10.21037/atm.2019.11.122), used GSE63409 together with an expression dataset to identify AML-relevant methylation/expression signals, immune and differentiation pathways, and candidate prognostic genes.

0.2 Setup and Dependencies

try(devtools::load_all(), silent = TRUE)
suppressPackageStartupMessages({
    try(library(CMEnt), silent = TRUE)
    library(minfi)
    library(data.table)
    library(ggplot2)
    library(dplyr)
    library(GenomicRanges)
})

getBenchmarkFile <- function(filename) {
    benchmark_file <- ""
    folder <- "microarray450k"
    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
}

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

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

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

0.2.1 Install External Dependencies

r <- getOption("repos")
if (is.null(r) || is.na(r["CRAN"]) || r["CRAN"] %in% c("@CRAN@", "")) {
  options(repos = c(CRAN = "https://cloud.r-project.org"))
}
if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}

if (!requireNamespace("tidyr", quietly = TRUE)) {
    install.packages("tidyr")
}

if (!requireNamespace("scales", quietly = TRUE)) {
    install.packages("scales")
}

if (!requireNamespace("IlluminaHumanMethylation450kanno.ilmn12.hg19", quietly = TRUE)) {
    BiocManager::install("IlluminaHumanMethylation450kanno.ilmn12.hg19", update = FALSE, ask = FALSE)
}

if (!requireNamespace("TxDb.Hsapiens.UCSC.hg19.knownGene", quietly = TRUE)) {
    BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene", update = FALSE, ask = FALSE)
}

if (!requireNamespace("org.Hs.eg.db", quietly = TRUE)) {
    BiocManager::install("org.Hs.eg.db", update = FALSE, ask = FALSE)
}

if (!requireNamespace("missMethyl", quietly = TRUE)) {
    BiocManager::install("missMethyl", update = FALSE, ask = FALSE)
}

0.3 Load Data

Covariates are: - predictedSex: the sex predicted from methylation data, which can be a strong confounder in methylation analyses. - subject.id.ch1: the patient ID - source_name_ch1: the cell type

covariates <- c("predictedSex", "subject.id.ch1", "source_name_ch1")
sample_group_col <- "subject.status.ch1"

beta <- getBenchmarkFile("GSE63409_beta_values.tsv.gz")
samplesheet <- getBenchmarkFile("GSE63409_samplesheet.tsv")

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


pheno <- fread(samplesheet)
pheno[[sample_group_col]] <- factor(pheno[[sample_group_col]], levels = c("normal", "AML"))
pheno <- as.data.frame(pheno)
rownames(pheno) <- pheno$V1
array_type <- "450K"
genome <- "hg19"

beta_handler <- CMEnt::getBetaHandler(beta, array = array_type, genome = genome)
beta_mat <- as.matrix(beta_handler$getBeta())
locs <- beta_handler$getBetaLocs()
pheno$casecontrol <- pheno[[sample_group_col]] == "AML"

0.4 Find Differentially Methylated Positions

dmps_file <- getBenchmarkFile("dmps.rds")
if (!file.exists(dmps_file)) {
    dmps <- CMEnt::findDMPsArray(
        beta_handler,
        pheno,
        sample_group_col = sample_group_col,
        id_col = "V1",
        covariates = covariates,
        array = array_type,
        genome = genome,
        njobs = max(1L, min(8L, parallel::detectCores(logical = TRUE) - 1L))
    )
    # sort DMPs by location
    rownames(dmps) <- dmps$site_id
    dmps <- dmps[rownames(locs), , drop = FALSE]
    dmps <- dmps[!is.na(dmps$pval), ]
    dmps$qval <- p.adjust(dmps$pval, method = "BH")

    case_idx <- which(pheno$casecontrol)
    ctrl_idx <- which(!pheno$casecontrol)
    dmps$delta_beta <- rowMeans(beta_mat[rownames(dmps), case_idx, drop = FALSE]) -
        rowMeans(beta_mat[rownames(dmps), ctrl_idx, drop = FALSE])

    dmps_gr <- makeGRangesFromDataFrame(
        data.frame(
            chr = locs$chr,
            start = locs$start,
            end = locs$start + 1L,
            pval = dmps$pval,
            qval = dmps$qval,
            stat = dmps$qval,
            delta_beta = dmps$delta_beta
        ),
        keep.extra.columns = TRUE
    )
    names(dmps_gr) <- rownames(dmps)
    saveRDS(dmps_gr, dmps_file)
} else {
    dmps_gr <- readRDS(dmps_file)
}
dmps <- as.data.frame(dmps_gr)
sig_dmps <- dmps[dmps$qval < 0.05 & abs(dmps$delta_beta) >= 0.05, ]
cat(sprintf("Found %d significant DMPs (FDR < 0.05, |delta beta| >= 0.05)\n", nrow(sig_dmps)))

Found 72556 significant DMPs (FDR < 0.05, |delta beta| >= 0.05)

0.5 Method 1: CMEnt

cment_file <- getBenchmarkFile("dmrs.CMEnt.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))
    cat(sprintf("parallel::detectCores(): %d\n", parallel::detectCores(logical = TRUE)))
    cat("Starting CMEnt::buildDMRs()...\n")
    flush.console()
    start_time <- Sys.time()
    dmrs_cment <- CMEnt::buildDMRs(
        beta = beta_handler,
        seeds = sig_dmps,
        pheno = pheno,
        sample_group_col = sample_group_col,
        casecontrol_col = "casecontrol",
        testing_mode  = "parametric",
        covariates = covariates,
        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",
        njobs = cment_njobs,
        verbose = 3
    )
    end_time <- Sys.time()
    time_cment <- end_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 8654 DMRs in 794 seconds

0.6 CMEnt Region Summary

if (!all(c("in_promoter_of", "in_gene_body_of", "delta_beta_promoter", "delta_beta_gene_body") %in% colnames(mcols(dmrs_cment))) &&
    requireNamespace("TxDb.Hsapiens.UCSC.hg19.knownGene", quietly = TRUE) &&
    requireNamespace("org.Hs.eg.db", quietly = TRUE)) {
    dmrs_cment <- CMEnt::annotateDMRsWithGenes(dmrs_cment, genome = genome)
}
for (gene_col in c("in_promoter_of", "in_gene_body_of")) {
    if (!gene_col %in% colnames(mcols(dmrs_cment))) {
        mcols(dmrs_cment)[[gene_col]] <- NA_character_
    }
}
for (delta_col in c("delta_beta_promoter", "delta_beta_gene_body")) {
    if (!delta_col %in% colnames(mcols(dmrs_cment))) {
        mcols(dmrs_cment)[[delta_col]] <- NA_real_
    }
}

dmr_delta_col <- if ("sites_delta_beta" %in% colnames(mcols(dmrs_cment))) {
    "sites_delta_beta"
} else {
    "delta_beta"
}
dmr_delta <- as.numeric(mcols(dmrs_cment)[[dmr_delta_col]])
dmr_direction <- ifelse(dmr_delta > 0, "Hyper", "Hypo")
dmr_summary <- data.frame(
    DMRs = length(dmrs_cment),
    Median_width_bp = median(width(dmrs_cment)),
    Mean_width_bp = round(mean(width(dmrs_cment)), 1),
    Total_coverage_bp = sum(width(dmrs_cment)),
    Median_delta_beta = round(median(dmr_delta, na.rm = TRUE), 4),
    Mean_abs_delta_beta = round(mean(abs(dmr_delta), na.rm = TRUE), 4),
    Hypermethylated_DMRs = sum(dmr_delta > 0.05, na.rm = TRUE),
    Hypomethylated_DMRs = sum(dmr_delta < -0.05, na.rm = TRUE)
)

knitr::kable(dmr_summary, caption = "CMEnt DMR summary for AML versus normal samples")
(#tab:cment_summary)CMEnt DMR summary for AML versus normal samples
DMRs Median_width_bp Mean_width_bp Total_coverage_bp Median_delta_beta Mean_abs_delta_beta Hypermethylated_DMRs Hypomethylated_DMRs
8654 1315 1740.2 15059673 0.0082 0.1087 2832 1633 
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)
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") %>%
  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 = ifelse(Region == region_levels[1], x_range[1], x_range[2]),
    hjust = ifelse(Region == region_levels[1], 0, 1)
  )

# 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_text(
    data = dmr_mood,
    aes(x = star_x, y = star_y, label = stars, hjust = hjust),
    inherit.aes = FALSE,
    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 = 3,
    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() +
  theme(legend.position = "bottom") +
  labs(
    title = "DMR Effect Size and Width by Gene Context",
    x = "Context-specific delta beta (AML - normal)",
    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 AML-minus-normal delta beta in promoter and gene-body contexts, restricted to absolute delta beta of at least 0.05.

0.7 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 5701 65.9
Gene body 7379 85.3
Promoter and gene body 5560 64.2
No annotated gene overlap 1134 13.1 
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. 

if (requireNamespace("IlluminaHumanMethylation450kanno.ilmn12.hg19", quietly = TRUE)) {
    library(IlluminaHumanMethylation450kanno.ilmn12.hg19)
    anno <- as.data.frame(minfi::getAnnotation(IlluminaHumanMethylation450kanno.ilmn12.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(rownames(sig_dmps), 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 15527 21.4
Significant DMPs N_Shelf 3821 5.3
Significant DMPs N_Shore 12188 16.8
Significant DMPs OpenSea 28474 39.2
Significant DMPs S_Shelf 3289 4.5
Significant DMPs S_Shore 9257 12.8
Sites inside DMRs Island 34457 44.3
Sites inside DMRs N_Shelf 2049 2.6
Sites inside DMRs N_Shore 15149 19.5
Sites inside DMRs OpenSea 13572 17.5
Sites inside DMRs S_Shelf 1403 1.8
Sites inside DMRs S_Shore 11136 14.3 
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.8 GO and KEGG Enrichment

missMethyl performs methylation-array-aware enrichment for CMEnt DMRs while correcting for probe-number bias across genes.

aml_term_pattern <- paste(
    c(
        "acute myeloid", "leukemia", "leukaemia", "myeloid", "hematopoiet", "haematopoiet",
        "stem cell", "differentiation", "immune", "leukocyte", "lymphocyte", "granulocyte",
        "monocyte", "macrophage", "cytokine", "chemokine", "phagosome", "JAK", "STAT",
        "HOX", "FLT3", "NPM1", "DNA methylation", "apoptosis", "cell cycle",
        "transcriptional misregulation"
    ),
    collapse = "|"
)

normalize_missmethyl_result <- function(x, 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(
        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
    )
}

missmethyl_long <- NULL
if (requireNamespace("missMethyl", quietly = TRUE) &&
    requireNamespace("IlluminaHumanMethylation450kanno.ilmn12.hg19", quietly = TRUE)) {
    library(IlluminaHumanMethylation450kanno.ilmn12.hg19)
    cached <- getBenchmarkFile("missmethyl_cment_enrichment_results.rds")
    if (file.exists(cached)) {
        missmethyl_long <- readRDS(cached)
    } else {
        missmethyl_long <- base::do.call(rbind, Filter(Negate(is.null), lapply(c("GO", "KEGG"), function(collection) {
            ret <- tryCatch(
                missMethyl::goregion(
                    regions = dmrs_cment,
                    all.cpg = rownames(locs),
                    collection = collection,
                    array.type = "450K",
                    plot.bias = FALSE,
                    prior.prob = TRUE
                ),
                error = function(e) {
                    warning(sprintf("missMethyl %s enrichment failed for CMEnt DMRs: %s", collection, e$message))
                    NULL
                }
            )
            normalize_missmethyl_result(ret, collection)
        })))
        saveRDS(missmethyl_long, cached)
    }
} else {
    cat("missMethyl or the 450K annotation package is required for GO/KEGG enrichment.\n")
}
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$AML_Context <- ifelse(grepl(aml_term_pattern, missmethyl_long$Term, ignore.case = TRUE), "AML-context term", "Other top term")

    missmethyl_top <- missmethyl_long |>
        group_by(Collection) |>
        arrange(FDR, .by_group = TRUE) |>
        slice_head(n = 10) |>
        ungroup()
    missmethyl_top$Term_Label <- short_label(missmethyl_top$Term, 70)
    missmethyl_top$Term_Collection <- factor(
        paste(missmethyl_top$Term_Label, missmethyl_top$Collection, sep = "___"),
        levels = rev(unique(paste(missmethyl_top$Term_Label, missmethyl_top$Collection, sep = "___")))
    )

    knitr::kable(
        missmethyl_top |>
            dplyr::select(Collection, Ontology, Term, N, DE, PValue, FDR, AML_Context) |>
            arrange(Collection, FDR),
        digits = 3,
        row.names = FALSE,
        caption = "Top missMethyl GO/KEGG enrichments from CMEnt DMRs"
    )
}
(#tab:enrichment_table)Top missMethyl GO/KEGG enrichments from CMEnt DMRs
Collection Ontology Term N DE PValue FDR AML_Context
GO GO cell periphery 5528 2122 0 0.000 Other top term
GO GO system development 3248 1426 0 0.000 Other top term
GO GO multicellular organismal process 6251 2353 0 0.000 Other top term
GO GO anatomical structure development 4906 1970 0 0.000 Other top term
GO GO multicellular organism development 3760 1586 0 0.000 Other top term
GO GO developmental process 5445 2140 0 0.000 Other top term
GO GO nervous system development 2074 978 0 0.000 Other top term
GO GO plasma membrane 5086 1931 0 0.000 Other top term
GO GO cell-cell signaling 1146 546 0 0.000 Other top term
GO GO gated channel activity 298 183 0 0.000 Other top term
KEGG KEGG Neuroactive ligand-receptor interaction 341 180 0 0.000 Other top term
KEGG KEGG Neuroactive ligand signaling 192 120 0 0.000 Other top term
KEGG KEGG Cadherin signaling 392 185 0 0.000 Other top term
KEGG KEGG Nicotine addiction 36 30 0 0.000 Other top term
KEGG KEGG Calcium signaling pathway 238 129 0 0.000 Other top term
KEGG KEGG Cell adhesion molecule (CAM) interaction 144 77 0 0.000 Other top term
KEGG KEGG Type I diabetes mellitus 42 27 0 0.004 Other top term
KEGG KEGG Hormone signaling 209 93 0 0.007 Other top term
KEGG KEGG Signaling pathways regulating pluripotency of stem cells 138 74 0 0.007 AML-context term
KEGG KEGG Morphine addiction 88 52 0 0.007 Other top term 
if (!is.null(missmethyl_long) && nrow(missmethyl_long) > 0L) {
    ggplot(missmethyl_top, aes(x = NegLog10FDR, y = Term_Collection, size = DE, color = AML_Context)) +
        geom_point(alpha = 0.9) +
        facet_wrap(~Collection, scales = "free_y") +
        scale_y_discrete(labels = function(x) sub("___.*$", "", x)) +
        scale_color_manual(values = c("AML-context term" = "#B8405E", "Other top term" = "grey50")) +
        theme_minimal() +
        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")
}
Faceted dot plots show enriched terms against negative log10 FDR. Dot size gives the number of differentially methylated genes; red dots mark AML-context terms and grey dots mark other top terms.

Figure 1: Top missMethyl GO and KEGG enrichments from CMEnt DMRs in GSE63409

focused_terms <- if (!is.null(missmethyl_long) && nrow(missmethyl_long) > 0L) {
    missmethyl_long |>
        filter(grepl(aml_term_pattern, Term, ignore.case = TRUE)) |>
        group_by(Collection) |>
        arrange(FDR, .by_group = TRUE) |>
        slice_head(n = 8) |>
        ungroup()
} else {
    data.frame()
}

if (nrow(focused_terms) > 0L) {
    focused_terms$Term_Label <- short_label(focused_terms$Term, 62)
    focused_terms$Term_Collection <- factor(
        paste(focused_terms$Term_Label, focused_terms$Collection, sep = "___"),
        levels = rev(unique(paste(focused_terms$Term_Label, focused_terms$Collection, sep = "___")))
    )
    knitr::kable(
        focused_terms |>
            dplyr::select(Collection, Ontology, Term, N, DE, PValue, FDR) |>
            arrange(Collection, FDR),
        digits = 3,
        row.names = FALSE,
        caption = "Top AML-relevant GO/KEGG terms"
    )
}
(#tab:focused_terms)Top AML-relevant GO/KEGG terms
Collection Ontology Term N DE PValue FDR
GO GO cell differentiation 3702 1479 0.000 0.000
GO GO neuron differentiation 1112 534 0.000 0.000
GO GO positive regulation of cell differentiation 738 334 0.000 0.000
GO GO central nervous system neuron differentiation 115 76 0.000 0.000
GO GO regulation of cell differentiation 1331 545 0.000 0.000
GO GO muscle cell differentiation 350 166 0.000 0.000
GO GO striated muscle cell differentiation 263 123 0.000 0.002
GO GO regulation of neuron differentiation 114 65 0.000 0.004
KEGG KEGG Signaling pathways regulating pluripotency of stem cells 138 74 0.000 0.007
KEGG KEGG Autoimmune thyroid disease 39 21 0.003 0.063
KEGG KEGG Hematopoietic cell lineage 90 37 0.013 0.138
KEGG KEGG Cytokine-cytokine receptor interaction 258 84 0.020 0.186
KEGG KEGG Intestinal immune network for IgA production 45 20 0.032 0.282
KEGG KEGG Viral protein interaction with cytokine and cytokine receptor 89 26 0.161 0.716
KEGG KEGG Chemokine signaling pathway 181 54 0.834 1.000
KEGG KEGG Cell cycle 153 22 1.000 1.000 
if (nrow(focused_terms) > 0L) {
    ggplot(focused_terms, aes(x = -log10(pmax(FDR, .Machine$double.xmin)), y = Term_Collection, fill = Collection)) +
        geom_col(width = 0.7) +
        facet_wrap(~Collection, scales = "free_y") +
        scale_y_discrete(labels = function(x) sub("___.*$", "", x)) +
        theme_minimal() +
        theme(axis.text.y = element_text(size = 8), legend.position = "none") +
        labs(title = "AML-Relevant Enrichment Terms", x = expression(-log[10]("FDR")), y = "")
} else {
    cat("No missMethyl terms matched the AML-focused keyword set.\n")
}
Faceted horizontal bars show negative log10 FDR for AML-relevant enrichment terms, with separate panels for GO and KEGG collections.

(#fig:focused_terms_plot)Top AML-relevant GO and KEGG enrichment terms from CMEnt DMRs in GSE63409.

0.9 Candidate CMEnt Gene Overlapping Loci

dmr_table <- as.data.frame(dmrs_cment)
dmr_table$dmr_index <- seq_along(dmrs_cment)
dmr_table$DMR <- paste0(dmr_table$seqnames, ":", dmr_table$start, "-", dmr_table$end)
dmr_table$delta_beta_for_ranking <- dmr_delta
dmr_table$abs_delta_beta <- abs(dmr_delta)
dmr_table$svm_score_for_ranking <- if ("score" %in% colnames(dmr_table)) dmr_table$score else 0
if (!"qval" %in% colnames(dmr_table)) {
    dmr_table$qval <- NA_real_
}
if (!"pval" %in% colnames(dmr_table)) {
    dmr_table$pval <- NA_real_
}
if (!"score" %in% colnames(dmr_table)) {
    dmr_table$score <- NA_real_
}

top_dmr_index <- dmr_table |>
    arrange(desc(svm_score_for_ranking), qval, desc(abs_delta_beta)) |>
    slice_head(n = 1) |>
    pull(dmr_index)

top_dmrs <- dmr_table |>
    filter(has_text(in_promoter_of) | has_text(in_gene_body_of)) |>
    arrange(desc(svm_score_for_ranking), qval, desc(abs_delta_beta)) |>
    dplyr::select(
        DMR,
        width,
        score,
        qval,
        pval,
        delta_beta_for_ranking,
        any_of(c("delta_beta_promoter", "delta_beta_gene_body", "seeds_num", "sites_num", "in_promoter_of", "in_gene_body_of"))
    ) |>
    head(20)

knitr::kable(top_dmrs, digits = 4, caption = "Top 20 promoter or gene-body overlapping CMEnt DMRs ranked by score")
(#tab:candidate_dmrs)Top 20 promoter or gene-body overlapping CMEnt DMRs ranked by score
DMR width score qval pval delta_beta_for_ranking delta_beta_promoter delta_beta_gene_body seeds_num sites_num in_promoter_of in_gene_body_of
chr1:cg14278300-cg09554443 chr1:167486978-167487762 785 0.6472 0e+00 0e+00 -0.0627 -0.3938 -0.0627 2 5 CD247 CD247
chr21:cg05126444-cg16704958 chr21:38629859-38630728 870 0.6286 0e+00 0e+00 -0.1452 NA -0.1452 2 5 NA VPS26C
chr8:cg25247520-cg00780520 chr8:128808017-128808554 538 0.6271 0e+00 0e+00 -0.3510 -0.3603 -0.3395 3 3 PVT1 PVT1
chr16:cg05673293-cg07605834 chr16:71475279-71475926 648 0.6231 0e+00 0e+00 0.4046 0.0978 0.4250 4 5 TLE7 TLE7
chr17:cg09489844-cg10168457 chr17:79880647-79882876 2230 0.6211 0e+00 0e+00 -0.0552 -0.0020 -0.0020 3 10 MAFG-AS1,MAFG,SIRT7,MILIP MAFG-AS1,MAFG,MILIP
chr10:cg10770832-cg09622269 chr10:76803270-76803925 656 0.6196 0e+00 0e+00 0.2904 NA 0.2755 5 5 NA DUSP29
chr19:cg24291974-cg19417168 chr19:4531638-4531961 324 0.6144 0e+00 0e+00 -0.4472 NA -0.3740 4 4 NA PLIN5
chr11:cg18431910-cg13228701 chr11:78614239-78614522 284 0.6116 1e-04 1e-04 0.2831 NA 0.2831 2 3 NA TENM4
chr11:cg26512142-cg05898618 chr11:2554198-2555576 1379 0.6050 0e+00 0e+00 -0.0686 NA -0.0520 4 11 NA KCNQ1
chr7:cg01572810-cg16810272 chr7:150901654-150902685 1032 0.5996 0e+00 0e+00 0.3312 -0.0022 0.3526 2 3 DRC11L DRC11L
chr1:cg16527057-cg25678231 chr1:157069242-157069611 370 0.5983 0e+00 0e+00 -0.2523 -0.2506 -0.2489 3 4 ETV3L ETV3L
chr17:cg26118142-cg16810621 chr17:3289363-3290164 802 0.5966 0e+00 0e+00 0.1574 0.0763 0.0763 3 6 OR3A2,OR1R1 OR3A2,OR1R1
chr1:cg16104450-cg06882058 chr1:243645911-243646787 877 0.5955 0e+00 0e+00 -0.1044 NA -0.1044 2 7 NA SDCCAG8
chr17:cg22350578-cg18960803 chr17:33470664-33470733 70 0.5937 2e-04 2e-04 0.0438 0.0438 NA 2 3 NLE1 NA
chr1:cg16100063-cg07357157 chr1:10690011-10690588 578 0.5931 0e+00 0e+00 -0.7241 NA -0.6760 3 4 NA PEX14
chr3:cg21463380-cg21526749 chr3:187389127-187389513 387 0.5929 0e+00 0e+00 -0.0854 -0.0854 NA 2 3 SST NA
chr15:cg21819984-cg08428878 chr15:101084507-101085341 835 0.5925 0e+00 0e+00 0.2879 0.3293 0.3515 7 8 CERS3 CERS3
chr1:cg11003133-cg17515347 chr1:159046391-159047163 773 0.5895 0e+00 0e+00 -0.5797 -0.5596 -0.5596 7 7 AIM2 AIM2
chr14:cg09578794-cg13252700 chr14:104690040-104690882 843 0.5889 0e+00 0e+00 -0.0280 0.0177 0.0177 4 6 LINC02691 LINC02691
chr10:cg04691264-cg04769618 chr10:29697905-29698685 781 0.5886 0e+00 0e+00 -0.0169 -0.0325 -0.0268 3 9 SVIL-AS1 SVIL-AS1

0.10 Top CMEnt DMR Visualization

CMEnt::plotDMR(
    dmrs_cment,
    dmr_index = top_dmr_index,
    beta = beta_handler,
    pheno = pheno,
    sample_group_col = sample_group_col,
    genome = "hg19",
    plot_motif = FALSE,
    array = array_type
)
Aligned panels show the selected DMR, seed and surrounding-site positions, annotated gene context, and a heatmap of sample beta values along the locus.

(#fig:top_dmr)Structure, gene context, and sample methylation of the top-ranked CMEnt DMR in GSE63409.

0.11 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 2: Genome-wide CMEnt DMR scores in GSE63409, colored by primary gene annotation class

0.12 CMEnt DMR Interaction Circos Visualization

CMEnt::plotDMRsCircos(
    dmrs_cment,
    genome = "hg19",
    array = array_type
)
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 3: Circos summary of CMEnt DMR locations, methylation direction, and detected interactions in GSE63409

0.13 Insights and Interpretation

This CMEnt re-analysis of GSE63409 extends the original Hu et al.  study from a gene-level methylation and expression biomarker analysis into a region-level interpretation of AML-associated methylation architecture. Whereas the original publication integrated methylation and expression data to nominate AML-relevant pathways and three prognostic genes, ATP11A, ITGAM, and ZNRF2, this notebook shows that the same methylation dataset contains a broad and spatially coherent regional signal: 72,556 significant DMPs were assembled into 8,654 CMEnt DMRs, covering about 15.1 Mb of the genome, with a slight excess of hypermethylated regions in AML.

A main new observation is that the regional signal is not randomly distributed across gene annotation space. Most CMEnt DMRs overlap annotated genes, especially gene bodies, and many overlap both promoter and gene-body annotations. This suggests that AML-associated methylation changes in this dataset are not limited to isolated promoter CpGs, but frequently form extended regulatory blocks spanning transcriptional units. The effect-size plot further indicates that both promoter-overlapping and gene-body-overlapping DMRs contain clear hypermethylated and hypomethylated subclasses, with gene-body DMRs contributing the larger number of region-level events.

Although significant DMPs are most frequent in OpenSea regions, CpGs retained inside CMEnt DMRs are strongly enriched in CpG islands and shores. This indicates that CMEnt preferentially captures locally coherent methylation changes in CpG-dense regulatory neighborhoods, rather than simply preserving the genome-wide distribution of significant single CpGs. In practical terms, the analysis reframes part of the AML methylation signal as coordinated regional dysregulation around CpG island-associated loci.

The enrichment results further reinforce the biological interpretation. The strongest GO terms are broad developmental, cell-periphery, cell-signaling, and differentiation categories, while AML-context terms include cell differentiation, regulation of differentiation, stem-cell pluripotency signaling, hematopoietic cell lineage, and cytokine-related pathways. Compared with the original study’s pathway-level analysis, CMEnt adds the observation that these biological themes are supported by coherent DMR structures rather than only by individually significant CpGs or methylation-expression gene intersections.

Overall, this analysis supports the original conclusion that AML is associated with dysregulated methylation in genes and pathways linked to hematopoietic differentiation, immune signaling, and development. Its main added value is that it converts the single-CpG signal into interpretable regional methylation structures, revealing CpG island-enriched DMR organization, promoter and gene-body-specific effect patterns, and candidate loci that can be prioritized for biological validation.

0.14 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] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.1 [2] CMEnt_0.99.14
[3] pathview_1.52.0
[4] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0 [5] scales_1.4.0
[6] minfi_1.58.0
[7] bumphunter_1.54.0
[8] locfit_1.5-9.12
[9] iterators_1.0.14
[10] foreach_1.5.2
[11] Biostrings_2.80.2
[12] XVector_0.52.0
[13] DSS_2.60.0
[14] BiocParallel_1.46.0
[15] bsseq_1.48.0
[16] SummarizedExperiment_1.42.0
[17] Biobase_2.72.0
[18] MatrixGenerics_1.24.0
[19] matrixStats_1.5.0
[20] ggplot2_4.0.3
[21] GenomicRanges_1.64.0
[22] Seqinfo_1.2.0
[23] IRanges_2.46.0
[24] S4Vectors_0.50.3
[25] BiocGenerics_0.58.1
[26] generics_0.1.4
[27] dplyr_1.2.1
[28] data.table_1.18.6.1
[29] testthat_3.3.2
[30] 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] xml2_1.6.0
[130] GenomicAlignments_1.48.0
[131] aplot_0.3.2
[132] hummingbird_1.22.0
[133] BSgenome_1.80.0
[134] ape_5.8-1
[135] viridisLite_0.4.3
[136] biovizBase_1.60.0
[137] TxDb.Hsapiens.UCSC.hg19.knownGene_3.22.1 [138] gridBase_0.4-7
[139] xtable_1.8-8
[140] interp_1.1-6
[141] plyr_1.8.9
[142] httr_1.4.9
[143] tools_4.6.1
[144] pkgbuild_1.4.8
[145] htmlTable_2.5.0
[146] checkmate_2.3.4
[147] nlme_3.1-170
[148] futile.logger_1.4.9
[149] lambda.r_1.2.6
[150] dbplyr_2.6.0
[151] ExperimentHub_3.2.2
[152] digest_0.6.39
[153] permute_0.9-10
[154] bookdown_0.48
[155] Matrix_1.7-6
[156] farver_2.1.2
[157] tzdb_0.5.0
[158] AnnotationFilter_1.36.0
[159] reshape2_1.4.5
[160] viridis_0.6.5
[161] yulab.utils_0.2.5
[162] DirichletMultinomial_1.54.1
[163] rpart_4.1.27
[164] glue_1.8.1
[165] cachem_1.1.0
[166] bspm_0.5.8
[167] VennDiagram_1.8.2
[168] BiocFileCache_3.2.0
[169] polyclip_1.10-7
[170] Hmisc_5.3-0
[171] sysfonts_0.8.9
[172] pkgload_1.5.3
[173] statmod_1.5.2
[174] fontBitstreamVera_0.1.1
[175] GEOquery_2.80.0
[176] httr2_1.3.0
[177] showtextdb_3.0
[178] gson_0.2.2
[179] dmrseq_1.32.0
[180] graphlayouts_1.2.5
[181] siggenes_1.86.0
[182] gtools_3.9.5
[183] readxl_1.6.0
[184] preprocessCore_1.74.0
[185] gridExtra_2.3.1
[186] shiny_1.14.0
[187] tidydr_0.0.7
[188] R.utils_2.13.0
[189] rhdf5filters_1.24.1
[190] RCurl_1.98-1.20
[191] memoise_2.0.1
[192] rmarkdown_2.32
[193] R.methodsS3_1.8.2
[194] svglite_2.2.2
[195] reshape_0.8.10
[196] fontLiberation_0.1.0
[197] illuminaio_0.54.0
[198] rstudioapi_0.19.0
[199] cluster_2.1.8.3
[200] hms_1.1.4
[201] DMRcate_3.8.0
[202] colorspace_2.1-4
[203] FNN_1.1.4.1
[204] rlang_1.3.0
[205] quadprog_1.5-8
[206] GenomeInfoDb_1.48.0
[207] DelayedMatrixStats_1.34.0
[208] sparseMatrixStats_1.24.0
[209] ggforce_0.5.0
[210] circlize_0.4.18
[211] ggtangle_0.1.4
[212] xfun_0.62
[213] ggseqlogo_0.2.2
[214] e1071_1.7-17
[215] aisdk_1.4.12
[216] abind_1.4-8
[217] GOSemSim_2.38.3
[218] tibble_3.3.1
[219] treeio_1.36.1
[220] Rhdf5lib_2.0.0
[221] readr_2.2.0
[222] futile.options_1.0.1
[223] bitops_1.1-0
[224] ps_1.9.3
[225] promises_1.5.0
[226] qvalue_2.44.0
[227] scatterpie_0.2.6
[228] inline_0.3.21
[229] RSQLite_3.53.3
[230] DelayedArray_0.38.2
[231] proxy_0.4-29
[232] GO.db_3.23.1
[233] compiler_4.6.1
[234] prettyunits_1.2.0
[235] beachmat_2.28.0
[236] graphite_1.58.0
[237] Rcpp_1.1.2
[238] benchmarkmeData_2.0.0
[239] fontquiver_0.2.1
[240] edgeR_4.10.6
[241] AnnotationHub_4.2.2
[242] usethis_3.2.2
[243] MASS_7.3-66
[244] progress_1.2.3
[245] R6_2.6.1
[246] fastmap_1.2.0
[247] ensembldb_2.36.1
[248] enrichit_0.2.6
[249] nnet_7.3-21
[250] gtable_0.3.6
[251] strex_2.1.0
[252] latticeExtra_0.6-31
[253] deldir_2.0-4
[254] htmltools_0.5.9
[255] bit64_4.8.6
[256] lifecycle_1.0.5
[257] S7_0.2.2
[258] processx_3.9.0
[259] callr_3.8.0
[260] Gviz_1.56.0
[261] restfulr_0.0.17
[262] sass_0.4.10
[263] vctrs_0.7.3
[264] DOSE_4.6.0
[265] ggfun_0.2.2
[266] bslib_0.12.0
[267] pillar_1.11.1
[268] GenomicFeatures_1.64.0
[269] magick_2.9.1
[270] otel_0.2.0
[271] jsonlite_2.0.0
[272] cigarillo_1.2.1
[273] GetoptLong_1.1.1