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

[sr320]

What and why

Follow-up to #15, merged on 2026-10-12. The first complete real methylation run, CIBNOR2021_HEAT_RRSS (12 WGBS libraries in 2 heat-resistant and 2 heat-susceptible families; draft #14), showed two problems with the family-level methylKit test from #15:

1. The q-values were anti-conservative. 397 tiles had q < 0.05, but 296 of them differed by 5 points or less, and q reached 1.7e-127. With 2 families per arm, methylKit's per-region McCullagh-Nelder dispersion falls below 1 wherever the families happen to agree. methylKit resets it to 1, and the test becomes a binomial test on summed reads.
2. Most of the data was discarded. Coverage was required per library and per strand (10× in all 12 libraries) before pooling, so only 23,350 1-kb tiles, about 4% of the genome, were testable. Each library has a mean of about 5× per strand.

DMR_POOLED_TILES (methylation.dmr_engine = 'limma' | 'dss'):

• merges the two strands of each CpG;
• pools by replicate_unit, and keeps CpGs with at least unit_min_coverage reads in every unit;
• sums them into tile_size tiles with at least tile_min_cpgs CpGs;
• tests with either:
  • limma: a moderated t on unit M-values, with empirical-Bayes variance and small-sample df; or
  • dss: DSS::DMLtest, a Wald test.

methylkit stays the default engine. The DMR output columns are unchanged, and the QC table records the engine, test and df_prior.

Testing

Synthetic Bismark CpG reports, both strands at 5–19× each: 12 libraries in 4 families × 3, real family-to-family variation (SD 12 points), 20 true tiles (+30 points) among 400.

Test (4 families)	True found	False positives (of 380)
methylKit, MN + F (#15)	0 / 20	0
DSS DMLtest	17 / 20	57
limma moderated t	4 / 20	1

• limma is calibrated: 5.0% of null tiles have p < 0.05. DSS's Wald test has no small-sample correction, so it's overconfident at 2 vs 2. That's why limma is the recommended engine, documented in workflows/methylation/README.md.
• The module reproduces the prototype exactly when extracted and run with DSS 2.58.0, bsseq 1.46.0, limma (R 4.5.1) and data.table.
• Versions reach the manifest: DMR_POOLED_TILES writes versions.yml (limma, DSS, bsseq, data.table), and it's wired to EMIT_MANIFEST. A ternary reads only the DMR process that ran, which satisfies tests/test_workflow_versions.py.
• Checks: pytest -q: 594 passed. The methylation demo (CI's command) passes.
• Real data: CIBNOR is re-running its DMR step with dmr_engine = 'limma' now. The 12 alignments, deduplications and extractions are reused from cache, because those modules are identical on main.

Later commits, all found by running CIBNOR end to end

• Report table capped at 500 rows (9c74a9f). Embedding all 319,088 tiles exhausted Quarto's 8 GB JavaScript heap (exit 133). The report now shows the 500 rows with the smallest q-values, and the full table is the standardized TSV.
• Bismark version captured in deduplication and extraction (faecd31). Those tools print "Deduplicator Version" and "Bismark Extractor Version", so versions.yml got an empty field.
• A standard error and p-value carried through to evidence (0b5621a). Methylation evidence had standard_error and p_value hard-coded to None, so no methylation study could pool. DMR_POOLED_TILES now writes pvalue, se_meth_diff_percent and analysis_method, standardizing passes them through when present, and the harmonizer reads them. There are 2 new tests in tests/test_methylation_harmonize.py. methylKit tables are unchanged.
• Gene-body units, dmr_features = 'genes' (f2a383f). Genome-wide tiles give many records per gene, so they collided in meta-analysis, and build-evidence-cards then failed for every study. Gene mode pools each gene's body and promoter CpGs into one unit, with intergenic tiles otherwise, and ANNOTATE_REGIONS keeps the gene and context those units carry.
• The limma SE now reproduces limma's p under a Wald test (7c8664e). Using |diff / t| dropped limma's small-sample df. One CIBNOR gene went from limma p 2.7e-4 to pooled p 4e-27, and 383 features became "significant", 53 of them promoted to multi-omics convergence. The SE is now |diff| / qnorm(1 − p/2).

End-to-end check: a full Pages-sequence rebuild with CIBNOR's evidence (register → demo → harmonize_committed_results.py → gene meta-analysis → cards → top candidates → dashboard → all-feature meta-analysis):

• it completes with no errors;
• CIBNOR contributes 30,732 emerging, resilience_associated candidates, none significant, as in limma;
• 0 multi-omics candidates;
• 889 high-priority candidates, unchanged from before CIBNOR's evidence;
• pytest -q on that tree: 601 passed.

Checklist

• pytest -q and ruff check src tests pass
• New logic has a test. Not in CI, as for Pool condition levels per arm, and fit DESeq2 on tanks with --replicate_unit #12 and Fit methylKit DMRs on families, not libraries (methylation.replicate_unit) #15: the logic is R inside a Nextflow module. It was tested on synthetic data as above.
• Relevant docs/ page updated (workflows/methylation/README.md)
• Study branches that need this change are listed below

Study branches waiting on this

• study/PRJNA690951/CIBNOR2021_HEAT_RRSS (#14): set dmr_engine = 'limma', min_coverage = 1, unit_min_coverage = 10, tile_min_cpgs = 10 in its config.

🤖 Generated with Claude Code

[sr320]

Note: this PR merged at its first commit (3926b34). The five later commits listed in the description (report cap, Bismark version capture, SE/p carry-through, gene-body units, Wald-matched SE) were pushed after the merge and are not on main. They're rebased, unchanged, in #74, which CIBNOR2021_HEAT_RRSS (#14) needs before it can merge.

🤖 Generated with Claude Code