diff --git a/modules/methylation/annotate_regions.nf b/modules/methylation/annotate_regions.nf index f88f968..5e330f8 100644 --- a/modules/methylation/annotate_regions.nf +++ b/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 diff --git a/modules/methylation/bismark_align.nf b/modules/methylation/bismark_align.nf index 4c7b1e4..951f734 100644 --- a/modules/methylation/bismark_align.nf +++ b/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 """ } diff --git a/modules/methylation/bismark_methylation_extractor.nf b/modules/methylation/bismark_methylation_extractor.nf index 4bdc221..76b9c9a 100644 --- a/modules/methylation/bismark_methylation_extractor.nf +++ b/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 """ } diff --git a/modules/methylation/dmr_pooled_tiles.nf b/modules/methylation/dmr_pooled_tiles.nf index 4bbc1ed..06a63f0 100644 --- a/modules/methylation/dmr_pooled_tiles.nf +++ b/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, diff --git a/modules/methylation/standardize_output.nf b/modules/methylation/standardize_output.nf index 873b1b7..9d9656c 100644 --- a/modules/methylation/standardize_output.nf +++ b/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) diff --git a/src/harmonize/methylation.py b/src/harmonize/methylation.py index 3f0844f..c4ea72f 100644 --- a/src/harmonize/methylation.py +++ b/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, diff --git a/tests/test_methylation_harmonize.py b/tests/test_methylation_harmonize.py new file mode 100644 index 0000000..fc72374 --- /dev/null +++ b/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() diff --git a/workflows/methylation/README.md b/workflows/methylation/README.md index ee56de0..7784378 100644 --- a/workflows/methylation/README.md +++ b/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 `__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 diff --git a/workflows/methylation/assets/report_template.qmd b/workflows/methylation/assets/report_template.qmd index 4a0ef62..65f7a19 100644 --- a/workflows/methylation/assets/report_template.qmd +++ b/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 diff --git a/workflows/methylation/main.nf b/workflows/methylation/main.nf index 5b973a6..6132dcc 100644 --- a/workflows/methylation/main.nf +++ b/workflows/methylation/main.nf @@ -48,6 +48,7 @@ def paramsAsJson() { tile_size : params.methylation?.tile_size, tile_step : params.methylation?.tile_step, dmr_engine : params.methylation?.dmr_engine ?: 'methylkit', + dmr_features : params.methylation?.dmr_features ?: 'tiles', unit_min_coverage : params.methylation?.unit_min_coverage, tile_min_cpgs : params.methylation?.tile_min_cpgs, replicate_unit : params.methylation?.replicate_unit ?: null, @@ -153,7 +154,10 @@ workflow { params.methylation.tile_min_cpgs, replicate_unit, design_sheet_ch, - dmr_engine + dmr_engine, + params.methylation.dmr_features ?: 'tiles', + annotation_gtf_ch, + params.methylation.promoter_upstream_bp ) : DMR_METHYLKIT( dmr_ch.map { it[0] }, diff --git a/workflows/methylation/nextflow.config b/workflows/methylation/nextflow.config index a1d9e7b..eacf5a1 100644 --- a/workflows/methylation/nextflow.config +++ b/workflows/methylation/nextflow.config @@ -60,6 +60,12 @@ params { // either, set min_coverage low (e.g. 1): coverage is then applied // after pooling, not per library. dmr_engine = 'methylkit' + // With 'limma' / 'dss': 'tiles' tests tile_size windows genome-wide; + // 'genes' tests one unit per gene (all CpGs in its body and + // promoter_upstream_bp promoter) plus tiles for intergenic CpGs. Use + // 'genes' when results must be one record per gene, e.g. to pool or + // link with RNA-seq evidence. + dmr_features = 'tiles' unit_min_coverage = 10 tile_min_cpgs = 10