Methylation gene-body units, poolable SEs, and three first-run fixes (the part of #67 that missed its merge)

modules/methylation/annotate_regions.nf

@@ -102,6 +102,17 @@ process ANNOTATE_REGIONS {
         gene_ids[i] <- ifelse(is.na(gid), "", gid)
     }
 
+    # A DMR step that classifies its own units (DMR_POOLED_TILES with
+    # dmr_features = 'genes') keeps that: a gene-body unit stays with its gene,
+    # and an intergenic tile stays gene-less even where it borders a gene
+    # window (its CpGs are only the ones outside it). Re-assigning by overlap
+    # would hand both the same gene and recreate an identifier collision.
+    preset <- if ("annotation_context" %in% names(dmr)) !is.na(dmr\$annotation_context) & nzchar(as.character(dmr\$annotation_context)) else rep(FALSE, nrow(dmr))
+    if (any(preset)) {
+        contexts[preset] <- as.character(dmr\$annotation_context[preset])
+        gid <- as.character(dmr\$gene_id[preset])
+        gene_ids[preset] <- ifelse(is.na(gid), "", gid)
+    }
     dmr\$annotation_context <- contexts
     dmr\$gene_id <- gene_ids
 

modules/methylation/bismark_align.nf

@@ -105,7 +105,7 @@ process BISMARK_DEDUPLICATE {
 
     cat <<-END_VERSIONS > versions.yml
     ${task.process}:
-        bismark: \$(deduplicate_bismark --version | grep 'Bismark Version' | sed 's/.*Bismark Version: v//')
+        bismark: \$(deduplicate_bismark --version | grep -oE 'Version: v[0-9.]+' | head -1 | sed 's/Version: v//')
     END_VERSIONS
     """
 }

modules/methylation/bismark_methylation_extractor.nf

@@ -43,7 +43,7 @@ process BISMARK_METHYLATION_EXTRACTOR {
 
     cat <<-END_VERSIONS > versions.yml
     ${task.process}:
-        bismark: \$(bismark_methylation_extractor --version | grep 'Bismark Version' | sed 's/.*Bismark Version: v//')
+        bismark: \$(bismark_methylation_extractor --version | grep -oE 'Version: v[0-9.]+' | head -1 | sed 's/Version: v//')
     END_VERSIONS
     """
 }

modules/methylation/dmr_pooled_tiles.nf

@@ -52,6 +52,9 @@ process DMR_POOLED_TILES {
     val replicate_unit
     path design_sheet
     val dmr_test              // 'limma' or 'dss'
+    val feature_units         // 'tiles' or 'genes' (gene bodies + promoters; intergenic tiles)
+    path annotation_gtf
+    val promoter_upstream_bp
 
     output:
     path "${study_id}_${comparison_id}_dmr_raw.tsv", emit: dmr_table
@@ -75,6 +78,10 @@ process DMR_POOLED_TILES {
     meth_diff_cut <- as.numeric(${meth_diff_cutoff})
     tile_size     <- as.integer(${tile_size})
     tile_min_cpgs <- as.integer(${tile_min_cpgs})
+    feature_units <- "${feature_units}"
+    stopifnot(feature_units %in% c("tiles", "genes"))
+    annotation_gtf <- "${annotation_gtf}"
+    promoter_bp   <- as.integer(${promoter_upstream_bp})
     unit_col      <- trimws("${replicate_unit}")
     design_path   <- "${design_sheet}"
 
@@ -128,19 +135,55 @@ process DMR_POOLED_TILES {
     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")
+    sum_by <- function(d, keys) d[, c(list(n_cpg = .N), lapply(.SD, sum)), by = keys, .SDcols = c(mcols, ucols)]
 
-    M   <- as.matrix(tiles[, ..mcols])
-    Cov <- M + as.matrix(tiles[, ..ucols])
+    if (feature_units == "genes") {
+        # Gene-body units: every CpG in a gene body or its promoter
+        # (promoter_bp upstream of the strand-aware start) counts toward that
+        # gene, so a gene contributes one estimate however many tiles it spans.
+        # CpGs outside every gene window fall back to tile_size tiles.
+        gtf <- fread(cmd = paste("grep -v '^#'", shQuote(annotation_gtf)), header = FALSE, sep = "\\t", quote = "",
+                     select = c(1, 3, 4, 5, 7, 9),
+                     col.names = c("chr", "type", "g_start", "g_end", "strand", "attr"))
+        gtf <- gtf[type == "gene"]
+        gtf[, gene_id := sub('.*gene_id "([^"]+)".*', "\\\\1", attr)]
+        gtf[, w_start := ifelse(strand == "-", g_start, pmax(1L, g_start - promoter_bp))]
+        gtf[, w_end   := ifelse(strand == "-", g_end + promoter_bp, g_end)]
+        win <- gtf[, .(chr, w_start, w_end, gene_id, g_start, g_end)]
+        setkey(win, chr, w_start, w_end)
+        cpg <- wide[, .(chr, s = pos, e = pos, row = .I)]
+        hits <- foverlaps(cpg, win, by.x = c("chr", "s", "e"), by.y = c("chr", "w_start", "w_end"),
+                          type = "within", nomatch = 0L)
+        genic_rows <- unique(hits\$row)
+        g <- sum_by(cbind(wide[hits\$row], gene_id = hits\$gene_id), "gene_id")
+        g <- merge(g, unique(win[, .(gene_id, chr, start = w_start, end = w_end)]), by = "gene_id")
+        g <- g[!duplicated(gene_id)]
+        g[, `:=`(annotation_context = "gene_body", unit_kind = "gene")]
+        inter <- wide[!seq_len(nrow(wide)) %in% genic_rows]
+        inter[, tile := (pos - 1L) %/% tile_size]
+        t <- sum_by(inter, c("chr", "tile"))
+        t[, `:=`(start = as.integer(tile * tile_size + 1L), end = as.integer((tile + 1L) * tile_size),
+                 gene_id = "", annotation_context = "intergenic", unit_kind = "tile")]
+        t[, tile := NULL]
+        units <- rbindlist(list(g, t), use.names = TRUE, fill = TRUE)
+        message(sprintf("gene-body units: %d genes; intergenic CpGs in %d tiles", nrow(g), nrow(t)))
+    } else {
+        wide[, tile := (pos - 1L) %/% tile_size]
+        units <- sum_by(wide, c("chr", "tile"))
+        units[, `:=`(start = as.integer(tile * tile_size + 1L), end = as.integer((tile + 1L) * tile_size),
+                     gene_id = NA_character_, annotation_context = NA_character_, unit_kind = "tile")]
+        units[, tile := NULL]
+    }
+    units <- units[n_cpg >= tile_min_cpgs]
+    message(sprintf("%s units with >= %d CpGs: %d", feature_units, tile_min_cpgs, nrow(units)))
+    if (nrow(units) == 0) stop("no units pass the coverage and CpG-count filters")
+
+    M   <- as.matrix(units[, ..mcols])
+    Cov <- M + as.matrix(units[, ..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") {
@@ -151,27 +194,41 @@ process DMR_POOLED_TILES {
         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
-        )
+        diff <- rowMeans(beta[, treat_units, drop = FALSE]) - rowMeans(beta[, control_units, drop = FALSE])
+        # treatment minus control, as for methylKit's meth.diff. diff.se is
+        # chosen so that a normal (Wald) test of diff / se - which is what
+        # inverse-variance meta-analysis applies - reproduces limma's own
+        # p-value. Using |diff / t| instead would drop limma's small-sample df:
+        # at 2 vs 2 units (df ~4), a moderated t of 10 is p ~3e-4 but a z of 10
+        # is p ~1e-23 (CIBNOR2021_HEAT_RRSS, first rebuild).
+        zq <- qnorm(1 - tt\$P.Value / 2)
+        res <- data.table(diff = diff, diff.se = ifelse(zq > 0, abs(diff) / zq, NA_real_),
+                          pval = tt\$P.Value, 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.
+        # Index positions: overlapping genes can share a start, and DSS needs
+        # unique loci; results are mapped back by index.
+        idx <- seq_len(nrow(units))
+        bs <- BSseq(chr = rep("units", nrow(units)), pos = idx, M = M, Cov = Cov, sampleNames = unit_levels)
         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]
+        res <- dml[match(idx, dml\$pos), .(diff, diff.se, pval, fdr)]
     }
 
+    label <- ifelse(units\$unit_kind == "gene", paste0("GENE_", units\$gene_id),
+                    paste0("TILE_", units\$chr, "_", units\$start))
     out <- data.frame(
-        region_id         = paste0("DMR_", seq_len(nrow(res))),
-        chrom             = res\$chr,
-        start             = res\$pos,
-        end               = res\$pos + tile_size - 1L,
+        region_id         = label,
+        chrom             = units\$chr,
+        start             = units\$start,
+        end               = units\$end,
+        gene_id           = ifelse(is.na(units\$gene_id), "", units\$gene_id),
+        annotation_context = ifelse(is.na(units\$annotation_context), "", units\$annotation_context),
+        n_cpg             = units\$n_cpg,
         meth_diff_percent = round(100 * res\$diff, 3),
         qvalue            = signif(res\$fdr, 4),
+        pvalue            = signif(res\$pval, 4),
+        se_meth_diff_percent = signif(100 * res\$diff.se, 4),
+        analysis_method   = paste0(ifelse(dmr_test == "limma", "limma", "DSS"), "_pooled_", feature_units),
         n_samples         = length(sample_ids),
         n_replicate_units = length(unit_levels),
         passes_significance_filter = !is.na(res\$fdr) & res\$fdr <= qvalue_cut &
@@ -187,7 +244,7 @@ process DMR_POOLED_TILES {
         n_samples             = length(sample_ids),
         n_candidate_regions   = nrow(out),
         n_significant_regions = sum(out\$passes_significance_filter),
-        dmr_mode              = "tile",
+        dmr_mode              = feature_units,
         tile_size             = tile_size,
         min_coverage          = unit_min_cov,
         qvalue_cutoff         = qvalue_cut,

modules/methylation/standardize_output.nf

@@ -40,6 +40,11 @@ process STANDARDIZE_OUTPUT {
         "region_id", "chrom", "start", "end", "gene_id",
         "annotation_context", "meth_diff_percent", "qvalue", "direction",
     ]
+    # Written only when the DMR table supplies them (DMR_POOLED_TILES does;
+    # methylKit does not), so a methylKit run's output is unchanged. Without a
+    # standard error and an unadjusted p-value, a methylation result cannot
+    # be pooled.
+    OPTIONAL_OUT_COLS = ["pvalue", "se_meth_diff_percent", "analysis_method"]
 
     def sniff_delimiter(path):
         with open(path, newline="") as fh:
@@ -92,7 +97,9 @@ process STANDARDIZE_OUTPUT {
         if direction not in ("hyper", "hypo"):
             direction = direction_from_diff(meth_diff)
 
+        extra = {c: r.get(c, "") for c in OPTIONAL_OUT_COLS if c in in_cols}
         out_rows.append({
+            **extra,
             "region_id": r.get("region_id", ""),
             "chrom": chrom,
             "start": start,
@@ -105,7 +112,8 @@ process STANDARDIZE_OUTPUT {
         })
 
     with open(out_path, "w", newline="") as fh:
-        writer = csv.DictWriter(fh, fieldnames=REQUIRED_OUT_COLS, delimiter="\\t")
+        out_cols = REQUIRED_OUT_COLS + [c for c in OPTIONAL_OUT_COLS if c in in_cols]
+        writer = csv.DictWriter(fh, fieldnames=out_cols, delimiter="\\t")
         writer.writeheader()
         writer.writerows(out_rows)
 

src/harmonize/methylation.py

@@ -52,9 +52,15 @@ def harmonize_methylation(
             feature_id_original = r["region_id"]
 
         meth_diff = float(r["meth_diff_percent"]) if pd.notna(r.get("meth_diff_percent")) else None
+        # Present when the DMR step reports them (pooled-tile engines), absent
+        # for methylKit tables, which then stay unpoolable as before.
+        se = r.get("se_meth_diff_percent")
+        pval = r.get("pvalue")
+        method = r.get("analysis_method")
+        method = method if isinstance(method, str) and method.strip() else analysis_method
 
         rows.append({
-            "evidence_id": make_evidence_id(study["study_id"], comparison["comparison_id"], feature_id_original, analysis_method),
+            "evidence_id": make_evidence_id(study["study_id"], comparison["comparison_id"], feature_id_original, method),
             "study_id": study["study_id"],
             "comparison_id": comparison["comparison_id"],
             "simulated": bool(study.get("simulated", False)),
@@ -67,10 +73,10 @@ def harmonize_methylation(
             "molecular_direction": DIRECTION_MAP.get(r.get("direction"), "ambiguous"),
             "effect_size": meth_diff,
             "effect_size_type": "methylation_diff_percent",
-            "standard_error": None,
+            "standard_error": float(se) if pd.notna(se) else None,
             "ci_lower": None,
             "ci_upper": None,
-            "p_value": None,
+            "p_value": float(pval) if pd.notna(pval) else None,
             "adjusted_p_value": float(r["qvalue"]) if pd.notna(r.get("qvalue")) else None,
             "sample_size": comparison["sample_size"],
             "tissue": comparison["tissue"],
@@ -79,7 +85,7 @@ def harmonize_methylation(
             "phenotype": comparison["phenotype"],
             "phenotype_direction": comparison["phenotype_direction"],
             "resilience_classification": comparison.get("resilience_classification"),
-            "analysis_method": analysis_method,
+            "analysis_method": method,
             "mapping_confidence": resolved.mapping_confidence,
             "quality_flags": compute_quality_flags(study, comparison, resolved.mapping_confidence),
             "source_file": source_ref,

tests/test_methylation_harmonize.py

@@ -0,0 +1,48 @@
+"""The methylation harmonizer carries a DMR table's standard error, p-value and
+method when the DMR step supplies them, and leaves older tables unchanged.
+
+Without a standard error and an unadjusted p-value a methylation result
+cannot enter inverse-variance pooling. methylKit tables carry neither;
+DMR_POOLED_TILES (dmr_engine limma or dss) writes both.
+"""
+from __future__ import annotations
+
+import pandas as pd
+
+from common import REPO_ROOT, STUDIES_DIR, load_yaml
+
+STUDY_ID = "GIGAS_PATH03"
+
+
+def _demo_table():
+    study = load_yaml(STUDIES_DIR / f"{STUDY_ID}.yaml")
+    (comp,) = study["comparisons"]
+    return comp["comparison_id"], pd.read_csv(REPO_ROOT / comp["results_file"], sep="\t")
+
+
+def test_methylkit_table_stays_unpoolable_and_labelled_methylkit(isolated_reports, tmp_path):
+    from harmonize.core import harmonize_processed_table
+
+    cid, table = _demo_table()
+    path = tmp_path / "plain_dmr_standardized.tsv"
+    table.to_csv(path, sep="\t", index=False)
+    ev = harmonize_processed_table(STUDY_ID, path, date_generated="2026-01-01", comparison_id=cid)
+    assert ev["standard_error"].isna().all()
+    assert ev["p_value"].isna().all()
+    assert set(ev["analysis_method"]) == {"methylKit"}
+
+
+def test_pooled_tile_columns_reach_the_evidence_table(isolated_reports, tmp_path):
+    from harmonize.core import harmonize_processed_table
+
+    cid, table = _demo_table()
+    table["pvalue"] = 0.01
+    table["se_meth_diff_percent"] = 2.5
+    table["analysis_method"] = "limma_pooled_tiles"
+    path = tmp_path / "pooled_dmr_standardized.tsv"
+    table.to_csv(path, sep="\t", index=False)
+    ev = harmonize_processed_table(STUDY_ID, path, date_generated="2026-01-01", comparison_id=cid)
+    assert (ev["standard_error"] == 2.5).all()
+    assert (ev["p_value"] == 0.01).all()
+    assert set(ev["analysis_method"]) == {"limma_pooled_tiles"}
+    assert ev["evidence_id"].str.contains("limma_pooled_tiles").all()

workflows/methylation/README.md

@@ -130,6 +130,16 @@ filter is the second guard.
   - `dss`: `DSS::DMLtest`, a beta-binomial model with shrunken dispersion
     and a Wald test.
 
+Units (`methylation.dmr_features`), for `limma` and `dss`:
+- `tiles` (default): `tile_size` windows genome-wide.
+- `genes`: one unit per gene, pooling every CpG in the gene body and its
+  `promoter_upstream_bp` promoter (strand-aware). CpGs outside every gene
+  window fall into intergenic tiles. Each gene then gives one record, so
+  harmonized methylation evidence links one-to-one to RNA-seq genes; with
+  `tiles`, many tiles share a gene and meta-analysis refuses to pool them
+  (identifier collision). `ANNOTATE_REGIONS` keeps the gene and context these
+  units arrive with.
+
 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:
@@ -176,7 +186,7 @@ which is clearly labeled simulated data.
 ## Standardized output schema
 
 `STANDARDIZE_OUTPUT` emits `<study_id>_<comparison_id>_dmr_standardized.tsv`
-with exactly these tab-separated columns, matching
+with these tab-separated columns (the last three optional), matching
 `data/demo/methylation/GIGAS_PATH03_OsHV1_challenge_vs_sham_dmr_demo.tsv` and
 what `src/harmonize/methylation.py` expects:
 
@@ -191,6 +201,12 @@ what `src/harmonize/methylation.py` expects:
 | `meth_diff_percent` | signed methylation difference, treatment − control |
 | `qvalue` | multiple-testing-adjusted p-value from methylKit (or as supplied in processed mode) |
 | `direction` | `hyper` \| `hypo`, derived from the sign of `meth_diff_percent` when not directly supplied |
+| `pvalue` | *optional*: unadjusted p-value. Written when the DMR step supplies it (`limma` / `dss` engines); methylKit does not |
+| `se_meth_diff_percent` | *optional*: standard error of `meth_diff_percent`. For `limma`, abs(diff) / qnorm(1 - p/2), so a normal (Wald) test, which is what meta-analysis applies, reproduces limma's own p-value and keeps its small-sample df. Using abs(diff / t) instead gave p ~1e-23 where limma gave 3e-4 at 2 vs 2 units |
+| `analysis_method` | *optional*: e.g. `limma_pooled_tiles`; the harmonizer uses it as the evidence `analysis_method` (default `methylKit`) |
+
+Without `pvalue` and `se_meth_diff_percent`, harmonized methylation
+evidence has no standard error and cannot enter inverse-variance pooling.
 
 The only rows genuinely dropped by `STANDARDIZE_OUTPUT` are ones missing
 `chrom`/`start`/`end` entirely (no genomic coordinates at all) — that count

workflows/methylation/assets/report_template.qmd

@@ -108,10 +108,22 @@ n_no_gene <- sum(is.na(dmr$gene_id) | trimws(as.character(dmr$gene_id)) == "")
 cat(n_no_gene, "of", nrow(dmr), "regions have no gene_id (kept as intergenic/unresolved).\n")
 ```
 
-## Full standardized region table
-
-```{r full-table}
-knitr::kable(dmr)
+## Standardized region table
+
+```{r full-table, results='asis'}
+# A whole-genome tile table (319,088 rows for CIBNOR2021_HEAT_RRSS) embedded
+# in a self-contained HTML page exhausts Quarto's 8 GB JavaScript heap, so
+# large tables are cut to the regions with the smallest q-values. The full
+# table is the standardized TSV in the run outputs.
+max_rows <- 500
+if (nrow(dmr) > max_rows) {
+    cat(sprintf("Showing the %d regions with the smallest q-values of %s. The full table is `%s`.\n\n",
+                max_rows, format(nrow(dmr), big.mark = ","), basename(params$standardized_table)))
+    shown <- head(dmr[order(dmr$qvalue), , drop = FALSE], max_rows)
+} else {
+    shown <- dmr
+}
+print(knitr::kable(shown, row.names = FALSE))
 ```
 
 ## Provenance