Test pooled family tiles with limma (or DSS), not per-region methylKit MN

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
+    """
+}

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

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

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