diff --git a/modules/methylation/dmr_pooled_tiles.nf b/modules/methylation/dmr_pooled_tiles.nf new file mode 100644 index 0000000..4bbc1ed --- /dev/null +++ b/modules/methylation/dmr_pooled_tiles.nf @@ -0,0 +1,220 @@ +// Tile-level differential methylation on pooled replicate units +// (raw_reanalysis mode only), selected by params.methylation.dmr_engine = +// 'limma' (recommended) or 'dss'. +// +// Written for designs with very few replicate units per arm, where +// methylKit's per-region McCullagh-Nelder overdispersion is unstable: with 2 +// families per arm, a region whose families happen to agree gets an estimated +// dispersion below 1, methylKit resets it to 1, and the test falls back to a +// binomial test on summed read counts (CIBNOR2021_HEAT_RRSS, first pass: +// 296 of 397 tiles at q < 0.05 differed by 5 points or less). +// +// Two tests on the same pooled tiles: +// limma (default): moderated t on per-unit M-values (logit of methylation), +// variance shrunk across tiles by empirical Bayes, with small-sample df. +// dss: DSS::DMLtest, beta-binomial with shrunken dispersion, Wald test. +// On synthetic data with 2 vs 2 families, family-level variation and 20 true +// regions among 400 tiles: limma found 4/20 with 1/380 false positives, and +// 5.0% of null tiles had p < 0.05 (calibrated); DSS found 17/20 but with +// 57/380 false positives, because its Wald test has no small-sample +// correction. Prefer limma unless there are several units per arm. +// +// Steps, per library CpG report (Bismark cytosine report, coverage-filtered +// by COVERAGE_FILTER at params.methylation.min_coverage): +// 1. destrand: a minus-strand CpG at position p is the same CpG as the +// plus-strand C at p - 1; their counts are summed. +// 2. pool by replicate unit (e.g. family), when replicate_unit is set; +// otherwise each library is its own unit. +// 3. keep CpGs covered by >= unit_min_coverage reads in EVERY unit. +// 4. sum kept CpGs into non-overlapping tiles of tile_size bp; keep tiles +// with >= tile_min_cpgs CpGs. +// 5. the chosen test on the tiles, treatment units vs control units. +// meth_diff_percent is treatment minus control mean methylation (%); +// qvalue is the Benjamini-Hochberg FDR (limma) or DSS's FDR. + +process DMR_POOLED_TILES { + tag "${study_id}:${comparison_id}" + label 'process_high' + container 'bioconductor/bioconductor_docker:RELEASE_3_18' + publishDir "${params.outdir}/methylation/dmr", mode: params.publish_mode + + input: + val study_id + val comparison_id + path cx_reports // one coverage-filtered cytosine report per library + val sample_ids + val treatment_labels // 0 = control, 1 = treatment, matched to sample_ids + val unit_min_coverage + val qvalue_cutoff + val meth_diff_cutoff + val tile_size + val tile_min_cpgs + val replicate_unit + path design_sheet + val dmr_test // 'limma' or 'dss' + + output: + path "${study_id}_${comparison_id}_dmr_raw.tsv", emit: dmr_table + path "${study_id}_${comparison_id}_methylkit_qc.tsv", emit: qc_metrics + path "versions.yml", emit: versions + + script: + """ + cat <<-'EOF_R' > run_pooled_tiles.R + suppressPackageStartupMessages({ + library(data.table) + }) + dmr_test <- "${dmr_test}" + stopifnot(dmr_test %in% c("limma", "dss")) + + sample_ids <- strsplit("${sample_ids.join(',')}", ",")[[1]] + treatments <- as.integer(strsplit("${treatment_labels.join(',')}", ",")[[1]]) + cx_files <- strsplit("${cx_reports.join(',')}", ",")[[1]] + unit_min_cov <- as.integer(${unit_min_coverage}) + qvalue_cut <- as.numeric(${qvalue_cutoff}) + meth_diff_cut <- as.numeric(${meth_diff_cutoff}) + tile_size <- as.integer(${tile_size}) + tile_min_cpgs <- as.integer(${tile_min_cpgs}) + unit_col <- trimws("${replicate_unit}") + design_path <- "${design_sheet}" + + stopifnot(length(sample_ids) == length(treatments), length(sample_ids) == length(cx_files)) + # Nextflow stages files under their own names; match each report to its + # sample by the sample id it starts with. + cx_files <- vapply(sample_ids, function(s) { + hit <- cx_files[startsWith(basename(cx_files), paste0(s, "."))] + if (length(hit) != 1) stop("expected one report for ", s, ", found ", length(hit)) + hit + }, character(1)) + + units <- sample_ids + if (nzchar(unit_col)) { + design <- read.delim(design_path, colClasses = "character", check.names = FALSE) + if (!unit_col %in% names(design)) stop("replicate_unit column '", unit_col, "' is not in the design sheet") + key <- if ("run_accession" %in% names(design)) "run_accession" else "sample_id" + units <- design[[unit_col]][match(sample_ids, design[[key]])] + if (any(is.na(units) | !nzchar(units))) { + stop("no '", unit_col, "' for: ", paste(sample_ids[is.na(units) | !nzchar(units)], collapse = ", ")) + } + } + arms_per_unit <- tapply(treatments, units, function(x) length(unique(x))) + if (any(arms_per_unit > 1)) { + stop("replicate_unit(s) spanning both arms: ", paste(names(arms_per_unit)[arms_per_unit > 1], collapse = ", ")) + } + unit_levels <- sort(unique(units)) + unit_trt <- as.integer(tapply(treatments, units, function(x) x[1])[unit_levels]) + if (any(table(factor(unit_trt, levels = c(0, 1))) < 2)) { + stop("each arm needs at least 2 replicate units; got ", + paste(sprintf("%s=%d", unit_levels, unit_trt), collapse = ", ")) + } + + read_destranded <- function(f) { + d <- fread(cmd = paste("gzip -dc", shQuote(f)), header = FALSE, select = 1:5, + col.names = c("chr", "pos", "strand", "m", "u")) + d[strand == "-", pos := pos - 1L] + d[, .(m = sum(m), u = sum(u)), by = .(chr, pos)] + } + + pooled <- list() + for (u in unit_levels) { + libs <- which(units == u) + d <- rbindlist(lapply(cx_files[libs], read_destranded)) + d <- d[, .(m = sum(m), u = sum(u)), by = .(chr, pos)] + d <- d[(m + u) >= unit_min_cov] + setnames(d, c("m", "u"), paste0(c("m_", "u_"), u)) + pooled[[u]] <- d + message(sprintf("unit %s: %d libraries, %d destranded CpGs at >= %dx", u, length(libs), nrow(d), unit_min_cov)) + } + wide <- Reduce(function(a, b) merge(a, b, by = c("chr", "pos")), pooled) + message(sprintf("CpGs at >= %dx in all %d units: %d", unit_min_cov, length(unit_levels), nrow(wide))) + + wide[, tile := (pos - 1L) %/% tile_size] + mcols <- paste0("m_", unit_levels); ucols <- paste0("u_", unit_levels) + tiles <- wide[, c(list(n_cpg = .N), lapply(.SD, sum)), by = .(chr, tile), .SDcols = c(mcols, ucols)] + tiles <- tiles[n_cpg >= tile_min_cpgs] + message(sprintf("tiles of %d bp with >= %d CpGs: %d", tile_size, tile_min_cpgs, nrow(tiles))) + if (nrow(tiles) == 0) stop("no tiles pass the coverage and CpG-count filters") + + M <- as.matrix(tiles[, ..mcols]) + Cov <- M + as.matrix(tiles[, ..ucols]) + colnames(M) <- colnames(Cov) <- unit_levels + treat_units <- unit_levels[unit_trt == 1] + control_units <- unit_levels[unit_trt == 0] + tiles[, start := as.integer(tile * tile_size + 1L)] + df_prior <- NA_real_ + + if (dmr_test == "limma") { + suppressPackageStartupMessages(library(limma)) + beta <- (M + 0.5) / (Cov + 1) + mval <- log2(beta / (1 - beta)) + X <- model.matrix(~ factor(unit_trt, levels = c(0, 1))) + fit <- eBayes(lmFit(mval, X)) + tt <- topTable(fit, coef = 2, number = Inf, sort.by = "none") + df_prior <- fit\$df.prior + res <- data.table( + chr = tiles\$chr, pos = tiles\$start, + # treatment minus control, as for methylKit's meth.diff + diff = rowMeans(beta[, treat_units, drop = FALSE]) - rowMeans(beta[, control_units, drop = FALSE]), + fdr = tt\$adj.P.Val + ) + } else { + suppressPackageStartupMessages({ library(DSS); library(bsseq) }) + bs <- BSseq(chr = tiles\$chr, pos = tiles\$start, M = M, Cov = Cov, sampleNames = unit_levels) + # group1 = treatment, so diff = treatment - control. + dml <- as.data.table(DMLtest(bs, group1 = treat_units, group2 = control_units, smoothing = FALSE)) + res <- dml[tiles[, .(chr, start)], on = .(chr, pos = start), nomatch = 0L] + } + + out <- data.frame( + region_id = paste0("DMR_", seq_len(nrow(res))), + chrom = res\$chr, + start = res\$pos, + end = res\$pos + tile_size - 1L, + meth_diff_percent = round(100 * res\$diff, 3), + qvalue = signif(res\$fdr, 4), + n_samples = length(sample_ids), + n_replicate_units = length(unit_levels), + passes_significance_filter = !is.na(res\$fdr) & res\$fdr <= qvalue_cut & + abs(100 * res\$diff) >= meth_diff_cut, + stringsAsFactors = FALSE + ) + write.table(out, file = "${study_id}_${comparison_id}_dmr_raw.tsv", + sep = "\\t", quote = FALSE, row.names = FALSE) + + qc <- data.frame( + study_id = "${study_id}", + comparison_id = "${comparison_id}", + n_samples = length(sample_ids), + n_candidate_regions = nrow(out), + n_significant_regions = sum(out\$passes_significance_filter), + dmr_mode = "tile", + tile_size = tile_size, + min_coverage = unit_min_cov, + qvalue_cutoff = qvalue_cut, + meth_diff_cutoff = meth_diff_cut, + replicate_unit = ifelse(nzchar(unit_col), unit_col, NA), + n_replicate_units = length(unit_levels), + overdispersion = ifelse(dmr_test == "limma", "empirical-Bayes moderated variance", "DSS shrunken beta-binomial"), + test = ifelse(dmr_test == "limma", "limma moderated t on unit M-values, BH", "DSS::DMLtest (Wald)"), + engine = dmr_test, + df_prior = df_prior, + destranded = TRUE, + tile_min_cpgs = tile_min_cpgs, + n_cpgs_all_units = nrow(wide) + ) + write.table(qc, file = "${study_id}_${comparison_id}_methylkit_qc.tsv", + sep = "\\t", quote = FALSE, row.names = FALSE) + EOF_R + + Rscript run_pooled_tiles.R + + cat <<-END_VERSIONS > versions.yml + ${task.process}: + r-base: \$(R --version | head -n1 | sed 's/R version //; s/ .*//') + bioconductor-limma: \$(Rscript -e 'cat(as.character(packageVersion("limma")))') + bioconductor-dss: \$(Rscript -e 'cat(as.character(packageVersion("DSS")))') + bioconductor-bsseq: \$(Rscript -e 'cat(as.character(packageVersion("bsseq")))') + r-data.table: \$(Rscript -e 'cat(as.character(packageVersion("data.table")))') + END_VERSIONS + """ +} diff --git a/workflows/methylation/README.md b/workflows/methylation/README.md index 9b7d3aa..ee56de0 100644 --- a/workflows/methylation/README.md +++ b/workflows/methylation/README.md @@ -115,6 +115,34 @@ a region's estimated overdispersion falls below 1, methylKit resets it to 1, so a few null regions can still get small q-values; the `meth_diff_cutoff` filter is the second guard. +#### DMR engines (`methylation.dmr_engine`) + +- `methylkit` (default): `DMR_METHYLKIT`, as above. +- `limma` and `dss`: `DMR_POOLED_TILES` (`modules/methylation/dmr_pooled_tiles.nf`). + It merges the two strands of each CpG, pools libraries by `replicate_unit`, + and keeps CpGs with at least `unit_min_coverage` reads in **every unit** + (coverage is applied after pooling, so set `min_coverage` low, e.g. 1). It + sums them into `tile_size` tiles with at least `tile_min_cpgs` CpGs, then + tests treatment units against control units. + - `limma`: a moderated t on per-unit M-values (the logit of methylation), + with variance shrunk across tiles by empirical Bayes. Recommended when + there are only a few units per arm. + - `dss`: `DSS::DMLtest`, a beta-binomial model with shrunken dispersion + and a Wald test. + +Why `limma` is recommended for few units: on synthetic data with 2 vs 2 +families, real family-to-family variation and 20 true regions among 400 +tiles: + +| Test | True regions found | False positives (of 380) | +|---|---|---| +| methylKit (families, MN + F) | 0 | 0 | +| DSS (families) | 17 | 57 | +| limma (families) | 4 | 1 | + +limma's null p-values were calibrated (5.0% below 0.05). DSS's Wald test +has no small-sample correction, so it is overconfident at 2 vs 2. + ### `processed_results_harmonization` For studies where only a processed DMR/DML table is publicly available (the diff --git a/workflows/methylation/main.nf b/workflows/methylation/main.nf index 8357e7f..5b973a6 100644 --- a/workflows/methylation/main.nf +++ b/workflows/methylation/main.nf @@ -28,6 +28,7 @@ include { BISMARK_DEDUPLICATE } from '../../modules/methylation/bism include { BISMARK_METHYLATION_EXTRACTOR } from '../../modules/methylation/bismark_methylation_extractor.nf' include { COVERAGE_FILTER } from '../../modules/methylation/coverage_filter.nf' include { DMR_METHYLKIT } from '../../modules/methylation/dmr_methylkit.nf' +include { DMR_POOLED_TILES } from '../../modules/methylation/dmr_pooled_tiles.nf' include { ANNOTATE_REGIONS } from '../../modules/methylation/annotate_regions.nf' include { STANDARDIZE_OUTPUT } from '../../modules/methylation/standardize_output.nf' include { EMIT_MANIFEST } from '../../modules/methylation/emit_manifest.nf' @@ -46,6 +47,9 @@ def paramsAsJson() { dmr_mode : params.methylation?.dmr_mode, tile_size : params.methylation?.tile_size, tile_step : params.methylation?.tile_step, + dmr_engine : params.methylation?.dmr_engine ?: 'methylkit', + unit_min_coverage : params.methylation?.unit_min_coverage, + tile_min_cpgs : params.methylation?.tile_min_cpgs, replicate_unit : params.methylation?.replicate_unit ?: null, design_sheet : params.methylation?.design_sheet ?: null, bismark_align_args : params.methylation?.bismark_align_args ?: '', @@ -131,7 +135,27 @@ workflow { ? Channel.fromPath(params.methylation.design_sheet, checkIfExists: true).first() : Channel.value(file("${projectDir}/assets/NO_DESIGN_SHEET")) - dmr_result_ch = DMR_METHYLKIT( + def dmr_engine = params.methylation.dmr_engine ?: 'methylkit' + if (!(dmr_engine in ['methylkit', 'limma', 'dss'])) { + error "methylation.dmr_engine must be 'methylkit', 'limma' or 'dss', got '${dmr_engine}'" + } + dmr_result_ch = dmr_engine in ['limma', 'dss'] + ? DMR_POOLED_TILES( + dmr_ch.map { it[0] }, + dmr_ch.map { it[1] }, + dmr_ch.map { it[2] }, + dmr_ch.map { it[3] }, + dmr_ch.map { it[4] }, + params.methylation.unit_min_coverage, + params.methylation.qvalue_cutoff, + params.methylation.meth_diff_cutoff, + params.methylation.tile_size, + params.methylation.tile_min_cpgs, + replicate_unit, + design_sheet_ch, + dmr_engine + ) + : DMR_METHYLKIT( dmr_ch.map { it[0] }, dmr_ch.map { it[1] }, dmr_ch.map { it[2] }, @@ -179,7 +203,8 @@ workflow { .mix(BISMARK_DEDUPLICATE.out.versions) .mix(BISMARK_METHYLATION_EXTRACTOR.out.versions) .mix(COVERAGE_FILTER.out.versions) - .mix(DMR_METHYLKIT.out.versions) + // Only the DMR process that ran has outputs; the ternary reads that one. + .mix(dmr_engine in ['limma', 'dss'] ? DMR_POOLED_TILES.out.versions : DMR_METHYLKIT.out.versions) .mix(ANNOTATE_REGIONS.out.versions) .mix(STANDARDIZE_OUTPUT.out.versions) .collect() diff --git a/workflows/methylation/nextflow.config b/workflows/methylation/nextflow.config index 5e0bd27..a1d9e7b 100644 --- a/workflows/methylation/nextflow.config +++ b/workflows/methylation/nextflow.config @@ -50,6 +50,19 @@ params { bismark_align_args = '' extractor_args = '' + // DMR engine. 'methylkit' (default) tests libraries, or pooled units, + // with methylKit. 'limma' and 'dss' (modules/methylation/ + // dmr_pooled_tiles.nf) destrand CpGs, pool by replicate_unit, keep + // CpGs at >= unit_min_coverage in every unit, and sum them into + // tile_size tiles with >= tile_min_cpgs CpGs. 'limma' then fits a + // moderated t on unit M-values, which stayed calibrated at 2 vs 2 units + // in simulation; 'dss' uses DSS's Wald test, which did not. With + // either, set min_coverage low (e.g. 1): coverage is then applied + // after pooling, not per library. + dmr_engine = 'methylkit' + unit_min_coverage = 10 + tile_min_cpgs = 10 + // Coverage filtering min_coverage = 10