Benchmarking CMEnt on Repeated Simulated ONT Data
CMEnt vs DSS, DMRcate, DMRseq, and hummingbird across repeated simulateDMRs runs
2026-10-19
Contents
0.1 Introduction
This vignette repeats the simulated ONT benchmark over multiple independently
simulated datasets. Each dataset is generated from the same augmented ONT BSseq
input with CMEnt::simulateDMRs(), using simulation parameters drawn under a
master seed of 42. Outputs are cached per simulation and method, and summary
plots show medians with observed interval ranges across simulations.
Methods:
- CMEnt: seed-guided, correlation-based DMR assembly
- DSS: smoothed two-group beta-binomial testing and native DMR calling
- DMRcate: Gaussian kernel smoothing-based DMR detection
- DMRseq: BSseq-native region-level inference
- hummingbird: HMM-based DMR detection
DSS uses DMLtest(smoothing = TRUE, smoothing.span = 500) followed by
callDMR() per chromosome with its default site-level p-value threshold of 1e-5 and
minCG = 50, matching the benchmark’s minimum-site setting. This p-value
threshold is not a region-level FDR cutoff. Other DMR-calling settings retain
their DSS defaults. The DSS runtime includes both testing and region calling.
0.2 Setup and Dependencies
.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(bsseq)
library(data.table)
library(dplyr)
library(GenomicRanges)
library(ggplot2)
library(SummarizedExperiment)
library(S4Vectors)
})
getBenchmarkFile <- function(filename) {
folder <- "simulated_ont_multisim"
benchmark_file <- if (getwd() == "notebooks") {
file.path("cached_results", folder, filename)
} else if ("NAMESPACE" %in% dir()) {
file.path("notebooks", "cached_results", folder, filename)
} else {
file.path("cached_results", folder, filename)
}
dir.create(dirname(benchmark_file), showWarnings = FALSE, recursive = TRUE)
benchmark_file
}
getSimulationFile <- function(sim_id, filename) {
getBenchmarkFile(file.path(sprintf("sim_%02d", sim_id), filename))
}
read_cache <- function(path, params = NULL) {
if (!file.exists(path)) return(NULL)
out <- readRDS(path)
if (is.null(params) || identical(out$params, params)) out else NULL
}
save_cache <- function(path, value, params = NULL) {
saveRDS(list(value = value, params = params), path)
}
nsim <- suppressWarnings(as.integer(Sys.getenv("CMENT_BENCHMARK_NSIM", "10")))
nsim <- ifelse(is.na(nsim) || nsim < 1L, 10L, nsim)
plan_file <- getBenchmarkFile(paste0("simulation_parameters", nsim, ".tsv"))
if (file.exists(plan_file)) {
simulation_plan <- data.table::fread(plan_file)
} else {
set.seed(42)
simulation_plan <- data.frame(
simulation_id = seq_len(nsim),
seed = sample(0:1000, nsim),
delta_beta_max = round(runif(nsim, 0.2, 0.5), 2),
expected_corr = round(runif(nsim, 0.3, 0.7), 2)
)
data.table::fwrite(simulation_plan, plan_file, sep = "\t")
}
knitr::kable(simulation_plan, caption = "Simulation parameters generated under master seed 42")| simulation_id | seed | delta_beta_max | expected_corr |
|---|---|---|---|
| 1 | 560 | 0.34 | 0.66 |
| 2 | 996 | 0.42 | 0.36 |
| 3 | 320 | 0.48 | 0.70 |
| 4 | 152 | 0.28 | 0.68 |
| 5 | 73 | 0.34 | 0.33 |
| 6 | 227 | 0.48 | 0.51 |
| 7 | 145 | 0.49 | 0.46 |
| 8 | 633 | 0.24 | 0.66 |
| 9 | 48 | 0.34 | 0.48 |
| 10 | 127 | 0.37 | 0.63 |
0.2.1 Install External Dependencies
if (!requireNamespace("benchmarkme", quietly = TRUE)) install.packages("benchmarkme")
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
if (!requireNamespace("DSS", quietly = TRUE)) BiocManager::install("DSS", update = FALSE, ask = FALSE)
if (!requireNamespace("DMRcate", quietly = TRUE)) BiocManager::install("DMRcate", update = FALSE, ask = FALSE)
if (!requireNamespace("dmrseq", quietly = TRUE)) BiocManager::install("dmrseq", update = FALSE, ask = FALSE)
if (!requireNamespace("hummingbird", quietly = TRUE)) BiocManager::install("hummingbird", update = FALSE, ask = FALSE)
if (!requireNamespace("tidyr", quietly = TRUE)) install.packages("tidyr")
suppressPackageStartupMessages(library(DSS))0.3 Load Input Data and Simulate
input_bsseq_path <- getBenchmarkFile("ont_bsseq_augmented.rds")
zenodo_url <- "https://zenodo.org/records/20592944/files/ont_bsseq_augmented.rds?download=1"
input_bsseq <- NULL
base_simulation <- list(
delta_jitter = 0.1,
max_gap = 500L,
min_dmr_size = 500L,
max_dmr_size = 5000L,
min_sites = 50L,
max_sites = 2000L,
num_dmrs = 1000L,
genome = "hg38"
)
min_seeds <- 2L
min_sites <- 50L
max_pval <- 0.05
testing_mode <- "parametric"
ext_site_delta_beta <- 0.1
max_bridge_seeds_gaps <- 1L
max_bridge_extension_gaps <- 2L
max_lookup_dist <- 1000L
simulate_one <- function(plan_row) {
sim_file <- getSimulationFile(plan_row$simulation_id, "simulation.rds")
params <- c(base_simulation, list(
seed = plan_row$seed,
delta_beta_max = plan_row$delta_beta_max,
expected_corr = plan_row$expected_corr
))
cached <- read_cache(sim_file, params)
if (!is.null(cached)) return(cached$value)
cat(
sprintf(
"Simulating dataset %d with seed %d, delta_beta_max %.3f, expected_corr %.3f",
plan_row$simulation_id, plan_row$seed, plan_row$delta_beta_max, plan_row$expected_corr
)
)
if (is.null(input_bsseq)) {
if (!file.exists(input_bsseq_path)) {
curl::curl_download(zenodo_url, input_bsseq_path)
}
input_bsseq <<- sort(readRDS(input_bsseq_path))
stopifnot(inherits(input_bsseq, "BSseq"))
}
set.seed(plan_row$seed)
sim <- CMEnt::simulateDMRs(
beta = input_bsseq,
sample_group_col = "Sample_Group",
num_dmrs = base_simulation$num_dmrs,
delta_max0 = plan_row$delta_beta_max,
delta_jitter = base_simulation$delta_jitter,
expected_correlation = plan_row$expected_corr,
max_gap = base_simulation$max_gap,
min_dmr_size = base_simulation$min_dmr_size,
max_dmr_size = base_simulation$max_dmr_size,
min_sites = base_simulation$min_sites,
max_sites = base_simulation$max_sites,
verbose = 1
)
sim$pheno$Sample_ID <- rownames(sim$pheno)
sim$pheno$casecontrol <- sim$pheno$Sample_Group == "case"
SummarizedExperiment::colData(sim$simulated) <- S4Vectors::DataFrame(sim$pheno, row.names = sim$pheno$Sample_ID)
data.table::fwrite(sim$truth, getSimulationFile(plan_row$simulation_id, "truth.tsv"), sep = "\t")
data.table::fwrite(sim$pheno, getSimulationFile(plan_row$simulation_id, "samplesheet.tsv"), sep = "\t")
save_cache(sim_file, sim, params)
sim
}
NJOBS <- as.integer(
Sys.getenv(
"NJOBS",
unset = as.character(max(1L, min(8L, parallel::detectCores(logical = TRUE) - 1L)))
)
)
doParallel::registerDoParallel(cores = NJOBS)
benchmark_specs_file <- getBenchmarkFile("benchmark_specs.rds")
if (!file.exists(benchmark_specs_file)) {
benchmark_specs <- list(
cpu = benchmarkme::get_cpu(),
ram = benchmarkme::get_ram(),
platform = benchmarkme::get_platform_info(),
rversion = R.version.string,
njobs = NJOBS
)
saveRDS(benchmark_specs, benchmark_specs_file)
} else {
benchmark_specs <- readRDS(benchmark_specs_file)
}
cpu <- benchmark_specs$cpu
ram <- benchmark_specs$ram
platform <- benchmark_specs$platform
rversion <- benchmark_specs$rversion
njobs <- benchmark_specs$njobs
cat("CPU:", cpu$vendor, cpu$model, "(", cpu$no_of_cores, "cores )\n")CPU: GenuineIntel 13th Gen Intel(R) Core(TM) i7-13700H ( 20 cores )
cat("RAM:", capture.output(print(ram)), "\n")RAM: 18.4 GB
cat("OS:", platform$OS.type, "\n")OS: unix
cat("R Version:", rversion, "\n")R Version: R version 4.6.1 (2026-06-24)
cat("Using", njobs, "parallel jobs\n")Using 8 parallel jobs
0.4 Method Runners
as_truth_gr <- function(sim) {
makeGRangesFromDataFrame(
sim$truth,
seqnames.field = "seqnames",
start.field = "start",
end.field = "end",
keep.extra.columns = TRUE,
seqinfo = GenomeInfoDb::Seqinfo(genome = base_simulation$genome)
)
}
prepare_bsseq_data <- function(sim, plan_row) {
bs <- sort(sim$simulated)
pheno <- as.data.frame(SummarizedExperiment::colData(bs))
pheno$Sample_ID <- rownames(pheno)
pheno$Sample_Group <- factor(pheno$Sample_Group, levels = c("control", "case"))
pheno$casecontrol <- pheno$Sample_Group == "case"
SummarizedExperiment::colData(bs) <- S4Vectors::DataFrame(pheno, row.names = pheno$Sample_ID)
seeds_file <- getSimulationFile(plan_row$simulation_id, sprintf("seeds.p%s.tsv.gz", gsub("[^0-9]+", "", sprintf("%.3f", max_pval))))
if (!file.exists(seeds_file)) {
seeds <- CMEnt::findDMPsBSSeq(
bsseq = bs,
samplesheet = pheno,
sample_group_col = "Sample_Group",
id_col = "Sample_ID",
case_group = "case",
output_file = seeds_file,
njobs = NJOBS
)
} else {
seeds <- data.table::fread(seeds_file)
}
seeds <- as.data.frame(seeds, stringsAsFactors = FALSE)
seeds <- seeds[seeds$qval < max_pval, , drop = FALSE]
rownames(seeds) <- seeds$site_id
list(bs = bs, pheno = pheno, seeds = seeds)
}
run_dss_bsseq <- function(bs_obj, params) {
sample_groups <- SummarizedExperiment::colData(bs_obj)$Sample_Group
dml <- DSS::DMLtest(
BSobj = bs_obj,
group1 = which(sample_groups == "case"),
group2 = which(sample_groups == "control"),
smoothing = params$smoothing,
smoothing.span = params$smoothing_span,
ncores = NJOBS
)
# DSS 2.60 callDMR detects region boundaries from position gaps alone.
dmrs <- do.call(rbind, lapply(split(dml, dml$chr), function(dml_chr) {
DSS::callDMR(dml_chr, p.threshold = params$p_threshold, minCG = params$min_sites)
}))
if (is.null(dmrs) || nrow(dmrs) == 0L) return(GRanges())
sort(makeGRangesFromDataFrame(
dmrs,
seqnames.field = "chr",
keep.extra.columns = TRUE,
seqinfo = GenomeInfoDb::Seqinfo(genome = base_simulation$genome)
))
}
run_dmrcate_bsseq <- function(bs_obj, design, params) {
if (!requireNamespace("DMRcate", quietly = TRUE)) return(GRanges())
dmrcate_design <- as.matrix(design)
colnames(dmrcate_design) <- gsub("Sample_Group", "", colnames(dmrcate_design))
colnames(dmrcate_design)[1] <- "Intercept"
dmrcate_design <- edgeR::modelMatrixMeth(dmrcate_design)
coef_names <- colnames(dmrcate_design)
sample_cols <- grepl("^Sample[0-9]+$", coef_names)
contrast_candidates <- coef_names[!(sample_cols | coef_names %in% c("Intercept", "(Intercept)"))]
if (length(contrast_candidates) == 0L) return(GRanges())
selected_coef <- match(ifelse("case" %in% contrast_candidates, "case", tail(contrast_candidates, 1L)), coef_names)
ann <- DMRcate::sequencing.annotate(bs_obj, all.cov = TRUE, methdesign = dmrcate_design, coef = selected_coef, fdr = params$max_pval)
out <- DMRcate::extractRanges(
DMRcate::dmrcate(ann, lambda = params$lambda, C = params$C, min.cpgs = params$min_sites),
genome = base_simulation$genome
)
out
}
run_dmrseq_bsseq <- function(bs_obj, params) {
if (!requireNamespace("dmrseq", quietly = TRUE)) return(GRanges())
materialize <- function(x) {
bsseq::BSseq(
M = as.matrix(bsseq::getCoverage(x, type = "M")),
Cov = as.matrix(bsseq::getCoverage(x, type = "Cov")),
gr = bsseq::granges(x),
pData = as.data.frame(SummarizedExperiment::colData(x)),
sampleNames = colnames(x)
)
}
select_calibration_sites <- function(bs, params) {
window_sites <- as.integer(params$calibration_window_sites)
windows_per_chr <- as.integer(params$calibration_windows_per_chr)
window_sites <- max(window_sites, as.integer(params$min_sites))
windows_per_chr <- max(windows_per_chr, 1L)
gr <- bsseq::granges(bs)
seqnames_vec <- as.character(seqnames(gr))
chr_vec <- unique(seqnames_vec[!is.na(seqnames_vec) & nzchar(seqnames_vec)])
idx <- lapply(chr_vec, function(chr_name) {
chr_idx <- which(seqnames_vec == chr_name)
n_chr <- length(chr_idx)
if (n_chr <= window_sites * windows_per_chr) return(chr_idx)
max_start <- n_chr - window_sites + 1L
starts <- unique(as.integer(round(seq.int(1L, max_start, length.out = windows_per_chr))))
unlist(lapply(starts, function(start_i) {
chr_idx[seq.int(start_i, min(start_i + window_sites - 1L, n_chr))]
}), use.names = FALSE)
})
sort(unique(unlist(idx, use.names = FALSE)))
}
calibrate_cutoff <- function(bs, params, bpparam) {
cutoffs <- unique(as.numeric(c(params$cutoff, params$calibration_cutoffs)))
cutoffs <- sort(cutoffs[is.finite(cutoffs) & cutoffs > 0 & cutoffs <= 1], decreasing = TRUE)
if (length(cutoffs) == 0L) stop("No valid dmrseq cutoff values supplied")
calibration_idx <- select_calibration_sites(bs, params)
if (length(calibration_idx) < as.integer(params$min_sites)) return(params$cutoff)
old_seed_exists <- exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)
old_seed <- if (old_seed_exists) get(".Random.seed", envir = .GlobalEnv) else NULL
on.exit({
if (old_seed_exists) {
assign(".Random.seed", old_seed, envir = .GlobalEnv)
} else if (exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) {
rm(".Random.seed", envir = .GlobalEnv)
}
}, add = TRUE)
set.seed(as.integer(params$calibration_seed))
calibration_bs <- bs[calibration_idx, , drop = FALSE]
for (cutoff in cutoffs) {
cat(sprintf(
"dmrseq calibration: testing cutoff %.3f on %s sites with %d permutation(s)\n",
cutoff,
format(length(calibration_idx), big.mark = ","),
as.integer(params$calibration_max_perms)
))
fit <- tryCatch(
suppressMessages(withCallingHandlers(
dmrseq::dmrseq(
bs = calibration_bs,
testCovariate = "Sample_Group",
cutoff = cutoff,
minNumRegion = params$min_sites,
maxPerms = as.integer(params$calibration_max_perms),
verbose = FALSE,
BPPARAM = bpparam
),
warning = function(w) invokeRestart("muffleWarning")
)),
error = function(e) e
)
if (inherits(fit, "error") || is.null(fit) || length(fit) == 0L) next
qvals <- mcols(fit)$qval
finite_qvals <- sum(is.finite(qvals))
cat(sprintf("dmrseq calibration: cutoff %.3f produced %d finite q-value(s)\n", cutoff, finite_qvals))
if (finite_qvals > 0L) return(cutoff)
}
warning(
"dmrseq cutoff calibration did not produce finite q-values on the calibration subset; ",
"using the lowest tested cutoff for the full run."
)
tail(cutoffs, 1L)
}
bpparam <- CMEnt:::.makeBiocParallelParam(max(1L, min(8L, parallel::detectCores(logical = TRUE) - 1L)), progressbar = TRUE)
bs_obj <- materialize(bs_obj)
cutoff <- calibrate_cutoff(bs_obj, params, bpparam)
cat(sprintf("dmrseq: using calibrated cutoff %.3f for full run\n", cutoff))
dmrs_by_chr <- dmrseq::dmrseq(
bs = bs_obj,
testCovariate = "Sample_Group",
cutoff = cutoff,
minNumRegion = params$min_sites,
maxPerms = as.integer(params$max_perms),
BPPARAM = bpparam
)
if (is.null(dmrs_by_chr) || length(dmrs_by_chr) == 0L) stop("dmrseq returned no DMRs")
qvals <- mcols(dmrs_by_chr)$qval
if (any(is.finite(qvals) & qvals < params$qval)) {
dmrs_by_chr <- dmrs_by_chr[is.finite(qvals) & qvals < params$qval]
}
sort(dmrs_by_chr)
}
run_hummingbird_bsseq <- function(bs_obj, max_gap = 300L, min_sites = 50L, min_length = 50L) {
if (!requireNamespace("hummingbird", quietly = TRUE)) return(GRanges())
sample_df <- as.data.frame(SummarizedExperiment::colData(bs_obj))
ctrl_idx <- which(sample_df$Sample_Group == "control")
case_idx <- which(sample_df$Sample_Group == "case")
if (length(ctrl_idx) == 0L || length(case_idx) == 0L) return(GRanges())
row_gr <- bsseq::granges(bs_obj)
row_gr <- GPos(seqnames = seqnames(row_gr), pos = start(row_gr))
to_count_matrix <- function(x) {
x <- as.matrix(x)
x[!is.finite(x) | x < 0] <- 0
storage.mode(x) <- "integer"
x
}
all_seqnames <- as.character(seqnames(row_gr))
chr_vec <- unique(all_seqnames[!is.na(all_seqnames) & nzchar(all_seqnames)])
dmrs_by_chr <- lapply(chr_vec, function(chr_name) {
chr_idx <- which(all_seqnames == chr_name)
if (length(chr_idx) < as.integer(min_sites)) return(NULL)
chr_bs <- bs_obj[chr_idx, , drop = FALSE]
meth_mat <- bsseq::getCoverage(chr_bs, type = "M")
cov_mat <- bsseq::getCoverage(chr_bs, type = "Cov")
unmeth_mat <- cov_mat - meth_mat
control_se <- SummarizedExperiment::SummarizedExperiment(
rowRanges = row_gr[chr_idx],
assays = list(normM = to_count_matrix(meth_mat[, ctrl_idx, drop = FALSE]), normUM = to_count_matrix(unmeth_mat[, ctrl_idx, drop = FALSE]))
)
case_se <- SummarizedExperiment::SummarizedExperiment(
rowRanges = row_gr[chr_idx],
assays = list(
abnormM = to_count_matrix(meth_mat[, case_idx, drop = FALSE]),
abnormUM = to_count_matrix(unmeth_mat[, case_idx, drop = FALSE])
)
)
tryCatch({
fit <- hummingbird::hummingbirdEM(experimentInfoControl = control_se, experimentInfoCase = case_se, binSize = 40)
fit <- hummingbird::hummingbirdPostAdjustment(
emInfo = fit,
experimentInfoControl = control_se,
experimentInfoCase = case_se,
maxGap = as.integer(max_gap),
minCpGs = as.integer(min_sites),
minLength = as.integer(min_length)
)
dmrs <- fit$DMRs
if (!inherits(dmrs, "GenomicRanges") || length(dmrs) == 0L) return(NULL)
as(dmrs, "GRanges")
}, error = function(e) {
warning(sprintf("hummingbird failed (%s): %s", chr_name, e$message))
NULL
})
})
dmrs_by_chr <- Filter(function(x) !is.null(x) && length(x) > 0L, dmrs_by_chr)
if (length(dmrs_by_chr) == 0L) {
GRanges()
} else {
sort(unlist(GRangesList(dmrs_by_chr), use.names = FALSE))
}
}
time_method <- function(expr) {
invisible(gc(reset = TRUE, verbose = FALSE))
start <- Sys.time()
value <- force(expr)
elapsed <- Sys.time() - start
gc_stats <- gc(verbose = FALSE)
max_heap_mb <- if (ncol(gc_stats) >= 6L) sum(gc_stats[, 6L], na.rm = TRUE) else sum(gc_stats[, 2L], na.rm = TRUE)
list(dmrs = value, time = elapsed, max_r_heap_mb = max_heap_mb)
}0.5 Run Benchmark
run_one_simulation <- function(plan_row) {
sim_dir <- sprintf("sim_%02d", plan_row$simulation_id)
common <- list(delta_beta_max = plan_row$delta_beta_max, simulation_seed = plan_row$seed)
method_specs <- list(
CMEnt = list(
params = c(common, list(
min_seeds = min_seeds, min_sites = min_sites / 2,
max_pval = max_pval, testing_mode = testing_mode,
ext_site_delta_beta = ext_site_delta_beta,
max_bridge_seeds_gaps = max_bridge_seeds_gaps,
max_bridge_extension_gaps = max_bridge_extension_gaps,
max_lookup_dist = max_lookup_dist
)),
run = function(params) {
time_method({
beta_handler <- CMEnt::getBetaHandler(dat$bs)
CMEnt::buildDMRs(
beta = beta_handler,
seeds = dat$seeds,
pheno = dat$pheno,
seeds_id_col = "site_id",
sample_group_col = "Sample_Group",
casecontrol_col = "casecontrol",
genome = base_simulation$genome,
array = "NULL",
entanglement = "weak",
testing_mode = params$testing_mode,
max_bridge_seeds_gaps = params$max_bridge_seeds_gaps,
max_bridge_extension_gaps = params$max_bridge_extension_gaps,
min_seeds = params$min_seeds,
min_sites = params$min_sites,
max_pval = params$max_pval,
ext_site_delta_beta = params$ext_site_delta_beta,
max_lookup_dist = params$max_lookup_dist,
annotate_with_genes = FALSE,
extract_motifs = FALSE,
njobs = NJOBS
)
})
}
),
DSS = list(
params = c(common, list(min_sites = min_sites, p_threshold = 1e-5, smoothing = TRUE, smoothing_span = 500L)),
run = function(params) time_method(run_dss_bsseq(dat$bs, params))
),
DMRcate = list(
params = c(common, list(min_sites = min_sites, max_pval = max_pval, lambda = 1000, C = 2)),
run = function(params) time_method(run_dmrcate_bsseq(dat$bs, model.matrix(~ Sample_Group, data = dat$pheno), params))
),
DMRseq = list(
params = c(common, list(
min_sites = min_sites,
cutoff = 0.05,
calibration_cutoffs = c(0.05, 0.025, 0.01),
calibration_window_sites = 50000L,
calibration_windows_per_chr = 2L,
calibration_max_perms = 1L,
calibration_seed = 1L,
max_perms = 10L,
qval = max_pval
)),
run = function(params) time_method(run_dmrseq_bsseq(dat$bs, params))
),
hummingbird = list(
params = c(common, list(max_gap = 300L, min_sites = min_sites, min_length = 50L)),
run = function(params) {
time_method(
run_hummingbird_bsseq(dat$bs, max_gap = params$max_gap, min_sites = params$min_sites, min_length = params$min_length)
)
}
)
)
cached_methods <- lapply(names(method_specs), function(method) {
spec <- method_specs[[method]]
cache_file <- getBenchmarkFile(file.path(sim_dir, sprintf("dmrs.%s.rds", gsub("[^A-Za-z0-9]+", "_", method))))
cached <- read_cache(cache_file, spec$params)
if (
!is.null(cached) &&
inherits(cached$value$dmrs, "GenomicRanges") &&
length(cached$value$dmrs) > 0L &&
!is.null(cached$value$max_r_heap_mb)
) {
cached$value
} else {
NULL
}
})
names(cached_methods) <- names(method_specs)
summary_file <- getSimulationFile(plan_row$simulation_id, "simulation_summary.rds")
summary_params <- c(base_simulation, list(
seed = plan_row$seed,
delta_beta_max = plan_row$delta_beta_max,
expected_corr = plan_row$expected_corr,
max_pval = max_pval
))
cached_summary <- read_cache(summary_file, summary_params)
if (!is.null(cached_summary) && all(!vapply(cached_methods, is.null, logical(1)))) {
truth_dmrs <- cached_summary$value$truth
seeds_n <- cached_summary$value$seeds_n
} else {
sim <- simulate_one(plan_row)
truth_dmrs <- as_truth_gr(sim)
dat <- prepare_bsseq_data(sim, plan_row)
seeds_n <- nrow(dat$seeds)
save_cache(summary_file, list(truth = truth_dmrs, seeds_n = seeds_n), summary_params)
}
method_outputs <- lapply(names(method_specs), function(method) {
spec <- method_specs[[method]]
cached <- cached_methods[[method]]
out <- if (!is.null(cached)) {
cached
} else {
cache_file <- getBenchmarkFile(file.path(sim_dir, sprintf("dmrs.%s.rds", gsub("[^A-Za-z0-9]+", "_", method))))
cat(sprintf("Running method %s for simulation %d", method, plan_row$simulation_id))
value <- spec$run(spec$params)
if (length(value$dmrs) > 0L) {
save_cache(cache_file, value, spec$params)
}
value
}
list(method = method, dmrs = out$dmrs, time = out$time, max_r_heap_mb = out$max_r_heap_mb)
})
names(method_outputs) <- names(method_specs)
list(
simulation_id = plan_row$simulation_id,
params = plan_row,
truth = truth_dmrs,
seeds_n = seeds_n,
methods = method_outputs
)
}
benchmark_results <- lapply(seq_len(nrow(simulation_plan)), function(i) {
run_one_simulation(simulation_plan[i, , drop = FALSE])
})0.6 Comparative Analysis
get_dmr_stats <- function(dmrs, method, simulation_id) {
data.frame(
simulation_id = simulation_id,
Method = method,
Number_of_DMRs = length(dmrs),
Median_DMR_Width = ifelse(length(dmrs) == 0L, NA_real_, median(width(dmrs))),
Total_Coverage_bp = ifelse(length(dmrs) == 0L, 0, sum(width(reduce(dmrs)))),
Min_Width = ifelse(length(dmrs) == 0L, NA_real_, min(width(dmrs))),
Max_Width = ifelse(length(dmrs) == 0L, NA_real_, max(width(dmrs)))
)
}
evaluate_against_truth <- function(pred, truth) {
n_pred <- length(pred)
n_truth <- length(truth)
pred_hit <- sum(overlapsAny(pred, truth))
truth_hit <- sum(overlapsAny(truth, pred))
precision <- ifelse(n_pred == 0L, 0, pred_hit / n_pred)
recall <- ifelse(n_truth == 0L, 0, truth_hit / n_truth)
pred_red <- reduce(pred)
truth_red <- reduce(truth)
inter_bp <- sum(width(GenomicRanges::intersect(pred_red, truth_red)))
pred_bp <- sum(width(pred_red))
truth_bp <- sum(width(truth_red))
union_bp <- pred_bp + truth_bp - inter_bp
precision_bp <- ifelse(pred_bp == 0L, 0, inter_bp / pred_bp)
recall_bp <- ifelse(truth_bp == 0L, 0, inter_bp / truth_bp)
data.frame(
TP_regions = truth_hit,
FP_regions = n_pred - pred_hit,
FN_regions = n_truth - truth_hit,
Precision_regions = precision,
Recall_regions = recall,
F1_regions = ifelse(precision + recall == 0, 0, 2 * precision * recall / (precision + recall)),
Precision_bp = precision_bp,
Recall_bp = recall_bp,
F1_bp = ifelse(precision_bp + recall_bp == 0, 0, 2 * precision_bp * recall_bp / (precision_bp + recall_bp)),
IoU_bp = ifelse(union_bp == 0L, 0, inter_bp / union_bp)
)
}
format_mean_pm <- function(x, digits = 2L, suffix = "") {
x <- x[is.finite(x)]
if (length(x) == 0L) return(NA_character_)
sd_x <- if (length(x) > 1L) stats::sd(x) else 0
fmt <- paste0("%.", digits, "f")
values <- sprintf(fmt, c(mean(x), sd_x))
paste0(values[1L], suffix, " (+/- ", values[2L], suffix, ")")
}
summarise_mean_pm <- function(df, group_cols, value_cols, digits = 2L) {
metric_digits <- if (length(digits) == 1L) {
stats::setNames(rep(digits, length(value_cols)), value_cols)
} else {
stats::setNames(as.integer(digits), names(digits))
}
df |>
tidyr::pivot_longer(dplyr::all_of(value_cols), names_to = "Metric", values_to = "Value") |>
dplyr::group_by(dplyr::across(dplyr::all_of(c(group_cols, "Metric")))) |>
dplyr::summarise(
Value = format_mean_pm(Value, digits = metric_digits[dplyr::first(Metric)]),
.groups = "drop"
) |>
tidyr::pivot_wider(names_from = Metric, values_from = Value)
}
summarise_benchmark_table <- function(stats_all, runtime_all) {
stats_all |>
dplyr::select(
simulation_id,
Method,
dmrs_number = Number_of_DMRs,
median_width_bp = Median_DMR_Width,
min_width_bp = Min_Width,
max_width_bp = Max_Width,
Total_Coverage_bp
) |>
dplyr::mutate(total_coverage_kb = Total_Coverage_bp / 1000) |>
dplyr::select(-Total_Coverage_bp) |>
dplyr::left_join(
runtime_all |>
dplyr::rename(runtime_sec = Seconds, max_r_heap_mb = Max_R_Heap_MB),
by = c("simulation_id", "Method")
) |>
dplyr::rename(method = Method) |>
dplyr::mutate(
relative_runtime_pct = {
cment_mean_runtime <- mean(runtime_sec[method == "CMEnt" & is.finite(runtime_sec)])
if (is.finite(cment_mean_runtime) && cment_mean_runtime > 0) {
runtime_sec / cment_mean_runtime * 100
} else {
NA_real_
}
}
) |>
dplyr::group_by(method) |>
dplyr::summarise(
dmrs_number = format_mean_pm(dmrs_number, digits = 0L),
median_width_bp = format_mean_pm(median_width_bp, digits = 0L),
min_width_bp = format_mean_pm(min_width_bp, digits = 0L),
max_width_bp = format_mean_pm(max_width_bp, digits = 0L),
total_coverage_kb = format_mean_pm(total_coverage_kb, digits = 1L),
runtime_sec = format_mean_pm(runtime_sec, digits = 1L),
relative_runtime_pct = format_mean_pm(relative_runtime_pct, digits = 1L, suffix = "%"),
max_r_heap_mb = format_mean_pm(max_r_heap_mb, digits = 1L),
.groups = "drop"
)
}
stats_all <- do.call(rbind, lapply(benchmark_results, function(res) {
rbind(
do.call(rbind, lapply(res$methods, function(x) get_dmr_stats(x$dmrs, x$method, res$simulation_id))),
get_dmr_stats(res$truth, "Truth", res$simulation_id)
)
}))
truth_eval_all <- do.call(rbind, lapply(benchmark_results, function(res) {
do.call(rbind, lapply(res$methods, function(x) {
out <- evaluate_against_truth(x$dmrs, res$truth)
out$simulation_id <- res$simulation_id
out$Method <- x$method
out
}))
}))
runtime_all <- do.call(rbind, lapply(benchmark_results, function(res) {
do.call(rbind, lapply(res$methods, function(x) {
data.frame(
simulation_id = res$simulation_id,
Method = x$method,
Seconds = as.numeric(x$time, units = "secs"),
Max_R_Heap_MB = x$max_r_heap_mb
)
}))
}))0.6.1 Basic Statistics
stats_summary <- summarise_benchmark_table(stats_all, runtime_all)
knitr::kable(stats_summary, caption = "Benchmark statistics across simulations, mean (+/- SD)")| method | dmrs_number | median_width_bp | min_width_bp | max_width_bp | total_coverage_kb | runtime_sec | relative_runtime_pct | max_r_heap_mb |
|---|---|---|---|---|---|---|---|---|
| CMEnt | 968 (+/- 70) | 3764 (+/- 381) | 264 (+/- 70) | 11907 (+/- 2338) | 3660.5 (+/- 559.1) | 106.0 (+/- 31.3) | 100.0% (+/- 29.5%) | 3645.1 (+/- 126.0) |
| DMRcate | 337 (+/- 197) | 2453 (+/- 438) | 832 (+/- 141) | 6886 (+/- 1738) | 928.5 (+/- 640.7) | 164.2 (+/- 4.1) | 154.9% (+/- 3.8%) | 14400.9 (+/- 1189.3) |
| DMRseq | 886 (+/- 51) | 3024 (+/- 49) | 562 (+/- 150) | 11801 (+/- 2278) | 2754.6 (+/- 188.6) | 6651.6 (+/- 186.2) | 6275.2% (+/- 175.6%) | 8616.6 (+/- 1011.0) |
| DSS | 668 (+/- 174) | 3044 (+/- 413) | 266 (+/- 0) | 9288 (+/- 2222) | 2103.5 (+/- 725.8) | 532.1 (+/- 2.4) | 502.0% (+/- 2.2%) | 9010.9 (+/- 484.5) |
| Truth | 1000 (+/- 0) | 3509 (+/- 23) | 665 (+/- 109) | 4998 (+/- 1) | 3461.5 (+/- 20.8) | NA | NA | NA |
| hummingbird | 357 (+/- 19) | 2134 (+/- 53) | 432 (+/- 96) | 4724 (+/- 182) | 786.2 (+/- 45.4) | 29.1 (+/- 2.1) | 27.5% (+/- 2.0%) | 7325.9 (+/- 932.4) |
stats_plot <- stats_all |>
dplyr::select(Method, simulation_id, Number_of_DMRs, Total_Coverage_bp) |>
tidyr::pivot_longer(c(Number_of_DMRs, Total_Coverage_bp), names_to = "Metric", values_to = "Value")
ggplot(stats_plot, aes(Method, Value, colour = Method)) +
geom_point(alpha = 0.45, position = position_jitter(width = 0.08, height = 0)) +
stat_summary(fun = median, geom = "point", size = 2.6, colour = "black") +
stat_summary(fun.min = min, fun.max = max, geom = "errorbar", width = 0.25, colour = "black") +
facet_wrap(~ Metric, scales = "free_y") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none") +
labs(title = "DMR counts and coverage across simulations", x = "", y = "Value")
Figure 1: DMR counts and total genomic coverage by method across simulated ONT datasets
0.7 Ground-Truth Evaluation
regions_truth_eval <- summarise_mean_pm(
truth_eval_all,
"Method",
c("TP_regions", "FP_regions", "FN_regions", "Precision_regions", "Recall_regions", "F1_regions"),
digits = c(TP_regions = 0L, FP_regions = 0L, FN_regions = 0L, Precision_regions = 2L, Recall_regions = 2L, F1_regions = 2L)
)
bp_truth_eval <- summarise_mean_pm(
truth_eval_all,
"Method",
c("Precision_bp", "Recall_bp", "F1_bp", "IoU_bp"),
digits = 2
)
knitr::kable(regions_truth_eval, caption = "Macro region-level performance against ground truth, mean (+/- SD)")| Method | F1_regions | FN_regions | FP_regions | Precision_regions | Recall_regions | TP_regions |
|---|---|---|---|---|---|---|
| CMEnt | 0.97 (+/- 0.05) | 44 (+/- 86) | 10 (+/- 7) | 0.99 (+/- 0.01) | 0.96 (+/- 0.09) | 956 (+/- 86) |
| DMRcate | 0.48 (+/- 0.24) | 657 (+/- 203) | 0 (+/- 0) | 1.00 (+/- 0.00) | 0.34 (+/- 0.20) | 343 (+/- 203) |
| DMRseq | 0.96 (+/- 0.03) | 81 (+/- 58) | 0 (+/- 0) | 1.00 (+/- 0.00) | 0.92 (+/- 0.06) | 919 (+/- 58) |
| DSS | 0.79 (+/- 0.15) | 325 (+/- 184) | 5 (+/- 1) | 0.99 (+/- 0.00) | 0.67 (+/- 0.18) | 675 (+/- 184) |
| hummingbird | 0.52 (+/- 0.02) | 648 (+/- 18) | 0 (+/- 0) | 1.00 (+/- 0.00) | 0.35 (+/- 0.02) | 352 (+/- 18) |
knitr::kable(bp_truth_eval, caption = "Micro base-pair-level performance against ground truth, mean (+/- SD)")| Method | F1_bp | IoU_bp | Precision_bp | Recall_bp |
|---|---|---|---|---|
| CMEnt | 0.88 (+/- 0.06) | 0.78 (+/- 0.09) | 0.86 (+/- 0.02) | 0.91 (+/- 0.13) |
| DMRcate | 0.39 (+/- 0.23) | 0.27 (+/- 0.18) | 0.99 (+/- 0.00) | 0.27 (+/- 0.18) |
| DMRseq | 0.87 (+/- 0.03) | 0.77 (+/- 0.05) | 0.98 (+/- 0.00) | 0.78 (+/- 0.05) |
| DSS | 0.73 (+/- 0.19) | 0.60 (+/- 0.21) | 0.99 (+/- 0.00) | 0.60 (+/- 0.21) |
| hummingbird | 0.37 (+/- 0.02) | 0.23 (+/- 0.01) | 1.00 (+/- 0.00) | 0.23 (+/- 0.01) |
metric_plot <- truth_eval_all |>
dplyr::select(Method, simulation_id, F1_regions, F1_bp, IoU_bp) |>
tidyr::pivot_longer(c(F1_regions, F1_bp, IoU_bp), names_to = "Metric", values_to = "Value")
metric_plot$Metric <- factor(
metric_plot$Metric,
levels = c("F1_regions", "F1_bp", "IoU_bp"),
labels = c("F1 (Region)", "F1 (Base-pair)", "IoU (Base-pair)")
)
ggplot(metric_plot, aes(Method, Value, colour = Method)) +
geom_point(alpha = 0.5, position = position_jitter(width = 0.08, height = 0)) +
stat_summary(fun = median, geom = "point", size = 2.6, colour = "black") +
stat_summary(fun.min = min, fun.max = max, geom = "errorbar", width = 0.25, colour = "black") +
facet_wrap(~ Metric, nrow = 1) +
coord_cartesian(ylim = c(0, 1)) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none") +
labs(title = "Method performance vs ground truth across simulations", x = "", y = "Score")
Figure 2: Region-level F1, base-pair-level F1, and base-pair intersection over union against ground truth across simulated ONT datasets
0.7.1 Seed Recovery
seed_summary <- do.call(rbind, lapply(benchmark_results, function(res) {
data.frame(simulation_id = res$simulation_id, Seeds = res$seeds_n, Truth_DMRs = length(res$truth))
}))
knitr::kable(summarise_mean_pm(seed_summary, character(0), c("Seeds", "Truth_DMRs"), digits = 2), caption = "Seed and truth-set sizes, mean (+/- SD)")| Seeds | Truth_DMRs |
|---|---|
| 43332.30 (+/- 19953.97) | 1000.00 (+/- 0.00) |
ggplot(seed_summary, aes(simulation_id, Seeds)) +
geom_line() +
geom_point() +
theme_minimal() +
labs(title = "CMEnt seed counts by simulation", x = "Simulation", y = "Seeds passing FDR threshold")
Figure 3: CMEnt seed counts passing the FDR threshold in each simulated ONT dataset
0.7.2 Runtime and Memory
runtime_summary <- summarise_mean_pm(runtime_all, "Method", c("Seconds", "Max_R_Heap_MB"), digits = 2)
knitr::kable(runtime_summary, caption = "Runtime and maximum R heap across simulations, mean (+/- SD)")| Method | Max_R_Heap_MB | Seconds |
|---|---|---|
| CMEnt | 3645.09 (+/- 125.96) | 106.00 (+/- 31.27) |
| DMRcate | 14400.87 (+/- 1189.30) | 164.18 (+/- 4.06) |
| DMRseq | 8616.63 (+/- 1010.97) | 6651.56 (+/- 186.16) |
| DSS | 9010.91 (+/- 484.53) | 532.06 (+/- 2.38) |
| hummingbird | 7325.94 (+/- 932.43) | 29.13 (+/- 2.07) |
runtime_plot <- runtime_all |>
tidyr::pivot_longer(c(Seconds, Max_R_Heap_MB), names_to = "Metric", values_to = "Value")
ggplot(runtime_plot, aes(Method, Value, colour = Method)) +
geom_point(alpha = 0.45, position = position_jitter(width = 0.08, height = 0)) +
stat_summary(fun = median, geom = "point", size = 2.6, colour = "black") +
stat_summary(fun.min = min, fun.max = max, geom = "errorbar", width = 0.25, colour = "black") +
facet_wrap(~ Metric, scales = "free_y") +
scale_y_continuous(labels = scales::comma) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none") +
labs(title = "Runtime and maximum R heap across simulations", x = "", y = "Value")
(#fig:runtime_plot)Runtime and maximum R heap usage by method across simulated ONT datasets.
0.7.3 Width Distributions
width_data <- do.call(rbind, lapply(benchmark_results, function(res) {
sets <- c(lapply(res$methods, `[[`, "dmrs"), list(Truth = res$truth))
do.call(rbind, lapply(names(sets), function(method) {
if (length(sets[[method]]) == 0L) return(NULL)
data.frame(simulation_id = res$simulation_id, Method = method, Width = width(sets[[method]]))
}))
}))
ggplot(width_data, aes(Method, Width, fill = Method)) +
geom_violin(alpha = 0.55) +
geom_boxplot(width = 0.18, alpha = 0.55, outlier.size = 0.4) +
scale_y_log10(labels = scales::comma) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none") +
labs(title = "DMR width distributions pooled across simulations", x = "", y = "Width (bp, log scale)")
(#fig:dmr_width_distributions_plot)DMR width distributions by method and ground truth, pooled across simulated ONT datasets.
0.8 Summary
best_f1 <- truth_eval_all |>
dplyr::group_by(Method) |>
dplyr::summarise(Median_F1 = median(F1_regions), Median_IoU = median(IoU_bp), .groups = "drop")
cat(sprintf("- Best median region-level F1: **%s** (%.3f)\n", best_f1$Method[which.max(best_f1$Median_F1)], max(best_f1$Median_F1)))- Best median region-level F1: CMEnt (0.991)
cat(sprintf("- Best median base-pair IoU: **%s** (%.3f)\n", best_f1$Method[which.max(best_f1$Median_IoU)], max(best_f1$Median_IoU)))- Best median base-pair IoU: CMEnt (0.822)
0.9 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] DSS_2.60.0 BiocParallel_1.46.0
[3] bsseq_1.48.0 SummarizedExperiment_1.42.0
[5] Biobase_2.72.0 MatrixGenerics_1.24.0
[7] matrixStats_1.5.0 CMEnt_0.99.14
[9] ggplot2_4.0.3 GenomicRanges_1.64.0
[11] Seqinfo_1.2.0 IRanges_2.46.0
[13] S4Vectors_0.50.3 BiocGenerics_0.58.1
[15] generics_0.1.4 dplyr_1.2.1
[17] data.table_1.18.6.1 testthat_3.3.2
[19] BiocStyle_2.40.0
loaded via a namespace (and not attached):
[1] R.methodsS3_1.8.2
[2] dichromat_2.0-1
[3] ggseqlogo_0.2.2
[4] strex_2.1.0
[5] progress_1.2.3
[6] nnet_7.3-21
[7] DT_0.34.0
[8] Biostrings_2.80.2
[9] HDF5Array_1.40.0
[10] vctrs_0.7.3
[11] digest_0.6.39
[12] png_0.1-9
[13] shape_1.4.6.1
[14] proxy_0.4-29
[15] deldir_2.0-4
[16] permute_0.9-10
[17] magick_2.9.1
[18] MASS_7.3-66
[19] pkgdown_2.2.1
[20] reshape_0.8.10
[21] reshape2_1.4.5
[22] httpuv_1.6.17
[23] foreach_1.5.2
[24] bumphunter_1.54.0
[25] withr_3.0.3
[26] showtextdb_3.0
[27] xfun_0.62
[28] survival_3.8-12
[29] bedr_1.1.5
[30] doRNG_1.8.6.3
[31] memoise_2.0.1
[32] benchmarkme_1.0.8
[33] sysfonts_0.8.9
[34] GlobalOptions_0.1.4
[35] gtools_3.9.5
[36] R.oo_1.27.1
[37] Formula_1.2-6
[38] prettyunits_1.2.0
[39] KEGGREST_1.52.2
[40] promises_1.5.0
[41] otel_0.2.0
[42] httr_1.4.9
[43] restfulr_0.0.17
[44] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.1
[45] rhdf5filters_1.24.1
[46] ps_1.9.3
[47] rhdf5_2.56.1
[48] rstudioapi_0.19.0
[49] UCSC.utils_1.8.0
[50] base64enc_0.1-6
[51] processx_3.9.0
[52] curl_8.1.0
[53] h5mread_1.4.1
[54] quadprog_1.5-8
[55] ExperimentHub_3.2.2
[56] SparseArray_1.12.3
[57] xtable_1.8-8
[58] stringr_1.6.0
[59] desc_1.4.3
[60] doParallel_1.0.17
[61] evaluate_1.0.5
[62] S4Arrays_1.12.1
[63] BiocFileCache_3.2.0
[64] preprocessCore_1.74.0
[65] hms_1.1.4
[66] bookdown_0.48
[67] colorspace_2.1-4
[68] filelock_1.0.3
[69] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
[70] VennDiagram_1.8.2
[71] readr_2.2.0
[72] magrittr_2.0.5
[73] later_1.4.8
[74] lattice_0.23-1
[75] genefilter_1.94.0
[76] bspm_0.5.8
[77] XML_3.99-0.25
[78] class_7.3-24
[79] Hmisc_5.3-0
[80] pillar_1.11.1
[81] nlme_3.1-170
[82] iterators_1.0.14
[83] pwalign_1.8.0
[84] hummingbird_1.22.0
[85] gridBase_0.4-7
[86] caTools_1.18.4
[87] compiler_4.6.1
[88] beachmat_2.28.0
[89] stringi_1.8.9
[90] shinycssloaders_1.1.0
[91] showtext_0.9-8
[92] GenomicAlignments_1.48.0
[93] plyr_1.8.9
[94] crayon_1.5.3
[95] abind_1.4-8
[96] BiocIO_1.22.0
[97] locfit_1.5-9.12
[98] org.Hs.eg.db_3.23.1
[99] bit_4.6.0
[100] codetools_0.2-20
[101] openssl_2.4.2
[102] bslib_0.12.0
[103] e1071_1.7-17
[104] biovizBase_1.60.0
[105] GetoptLong_1.1.1
[106] plotly_4.12.1
[107] multtest_2.68.0
[108] mime_0.13
[109] splines_4.6.1
[110] circlize_0.4.18
[111] Rcpp_1.1.2
[112] dbplyr_2.6.0
[113] sparseMatrixStats_1.24.0
[114] interp_1.1-6
[115] knitr_1.52
[116] blob_1.3.0
[117] clue_0.3-68
[118] BiocVersion_3.23.1
[119] seqLogo_1.78.0
[120] AnnotationFilter_1.36.0
[121] fs_2.1.0
[122] checkmate_2.3.4
[123] DelayedMatrixStats_1.34.0
[124] pkgbuild_1.4.8
[125] Gviz_1.56.0
[126] tibble_3.3.1
[127] Matrix_1.7-6
[128] callr_3.8.0
[129] statmod_1.5.2
[130] tzdb_0.5.0
[131] pkgconfig_2.0.3
[132] tools_4.6.1
[133] cachem_1.1.0
[134] cigarillo_1.2.1
[135] RSQLite_3.53.3
[136] viridisLite_0.4.3
[137] DBI_1.3.0
[138] fastmap_1.2.0
[139] rmarkdown_2.32
[140] scales_1.4.0
[141] grid_4.6.1
[142] outliers_0.15
[143] Rsamtools_2.28.0
[144] AnnotationHub_4.2.2
[145] sass_0.4.10
[146] FNN_1.1.4.1
[147] BiocManager_1.30.27
[148] VariantAnnotation_1.58.0
[149] scrime_1.3.7
[150] rpart_4.1.27
[151] minfi_1.58.0
[152] farver_2.1.2
[153] yaml_2.3.12
[154] latticeExtra_0.6-31
[155] foreign_0.8-91
[156] rtracklayer_1.72.0
[157] illuminaio_0.54.0
[158] cli_3.6.6
[159] purrr_1.2.2
[160] siggenes_1.86.0
[161] GEOquery_2.80.0
[162] lifecycle_1.0.5
[163] askpass_1.2.1
[164] lambda.r_1.2.6
[165] backports_1.5.1
[166] annotate_1.90.0
[167] gtable_0.3.6
[168] rjson_0.2.23
[169] missMethyl_1.46.0
[170] limma_3.68.5
[171] jsonlite_2.0.0
[172] edgeR_4.10.6
[173] TFBSTools_1.50.0
[174] bitops_1.1-0
[175] benchmarkmeData_2.0.0
[176] bit64_4.8.6
[177] brio_1.1.5
[178] base64_2.0.2
[179] annotatr_1.38.3
[180] futile.options_1.0.1
[181] jquerylib_0.1.4
[182] R.utils_2.13.0
[183] DMRcate_3.8.0
[184] lazyeval_0.2.3
[185] shiny_1.14.0
[186] htmltools_0.5.9
[187] rappdirs_0.3.4
[188] formatR_1.14
[189] ensembldb_2.36.1
[190] tinytex_0.61
[191] glue_1.8.1
[192] TFMPvalue_1.0.0
[193] httr2_1.3.0
[194] XVector_0.52.0
[195] RCurl_1.98-1.20
[196] rprojroot_2.1.1
[197] mclust_6.1.3
[198] BSgenome_1.80.0
[199] jpeg_0.1-11
[200] futile.logger_1.4.9
[201] gridExtra_2.3.1
[202] igraph_2.3.4
[203] R6_2.6.1
[204] tidyr_1.3.2
[205] DelayedDataFrame_1.28.0
[206] labeling_0.4.3
[207] GenomicFeatures_1.64.0
[208] cluster_2.1.8.3
[209] rngtools_1.5.2
[210] regioneR_1.44.0
[211] pkgload_1.5.3
[212] Rhdf5lib_2.0.0
[213] beanplot_1.3.1
[214] GenomeInfoDb_1.48.0
[215] DirichletMultinomial_1.54.1
[216] DelayedArray_0.38.2
[217] tidyselect_1.2.1
[218] ProtGenerics_1.44.0
[219] dmrseq_1.32.0
[220] htmlTable_2.5.0
[221] xml2_1.6.0
[222] inline_0.3.21
[223] AnnotationDbi_1.74.0
[224] S7_0.2.2
[225] nor1mix_1.3-3
[226] htmlwidgets_1.6.4
[227] ComplexHeatmap_2.28.0
[228] RColorBrewer_1.1-3
[229] biomaRt_2.68.0
[230] rlang_1.3.0
[231] rentrez_1.2.4