Fit methylKit DMRs on families, not libraries (methylation.replicate_unit)

[sr320]

What and why

#14 (CIBNOR2021_HEAT_RRSS) is blocked on this. The study's thermal-tolerance phenotype is assigned per family: 12 WGBS libraries from 2 heat-resistant (RR52, RR59) and 2 heat-susceptible (SS05, SS35) families. DMR_METHYLKIT fitted RR vs SS on the 12 libraries as independent, so family-to-family variation would have passed for a phenotype effect. This is the methylation counterpart of --replicate_unit for RNA-seq (#12).

• methylation.replicate_unit names a column in methylation.design_sheet (a TSV with run_accession or sample_id; for CIBNOR, its committed samplesheet.tsv and family).
• What DMR_METHYLKIT does with it:
  • After unite, it sums each unit's methylated / unmethylated counts per region into one pseudo-library.
  • It tests the units with methylKit's McCullagh-Nelder overdispersion correction and an F test (overdispersion = "MN", test = "F"), so the error comes from family-to-family variation.
  • It refuses a unit spanning both arms, fewer than 2 units per arm, a missing column, or a library with no unit.
• Provenance: the run manifest's parameters include replicate_unit and design_sheet. The methylKit QC table records n_replicate_units, overdispersion and test. The DMR table gains n_replicate_units.
• Without replicate_unit, nothing changes. A placeholder design sheet (assets/NO_DESIGN_SHEET) is passed, the same pattern as the RNA-seq viral QC.

Also: four defects found by the first real run (second to fourth commits), plus alignment options (fifth)

• COVERAGE_FILTER used zcat, which on macOS only reads .Z files (can't stat: x.gz.Z), so every library would have failed. Now gzip -dc. On a gzipped CpG report it keeps exactly the rows an independent count does.
• BISMARK_METHYLATION_EXTRACTOR wrote a --CX report (CpG + CHG + CHH). For this genome that's about 193 M rows per library against about 24 M CpGs. The DMR step uses CpG only, and methylKit reads the whole file into memory before filtering by context, so 12 libraries would not fit on a 26 GB machine. --CX is dropped, and the output is now *.CpG_report.txt.gz. Non-CpG methylation would need its own path.
• BISMARK_GENOME_PREPARATION skipped the genome (third commit). Bismark only recognises .fa / .fasta, and NCBI genomes are .fna, so the first real run stopped with "does not contain any sequence files in FastA format". The genome is now copied in as genome.fa. Checked with Bismark 0.24.2: the bisulfite indexes build.
• BISMARK_ALIGN passed --basename with --multicore (fourth commit). Bismark 0.24 rejects that combination (exit 60). --multicore is kept, --basename is dropped, and the read-1-named outputs are renamed. On toy reads, alignment, the rename, deduplicate_bismark and the CpG-only extractor all ran and produced the files each next step consumes.

Also: per-study Bismark options (fifth commit)

The first real alignment mapped only 13.3–14.0% of read pairs uniquely across 6 CIBNOR libraries. A diagnostic on 50,000 trimmed pairs of one library:

Setting	Unique mapping	CHH methylation
default --score_min L,0,-0.2	15.5%	3.0%
--score_min L,0,-0.6	53.9%	2.9%
--non_directional	15.5%	3.0%
-X 1000	15.8%	3.0%
read 1 alone	37.2%	0.8%
read 2 alone (--pbat)	17.2%	5.0%

• The cause is the default alignment-score threshold, which is too strict for a genome as heterozygous as M. gigas. The library is directional, and the insert-size limit is not the problem.
• methylation.bismark_align_args and methylation.extractor_args pass options through as written and are recorded in the manifest. Both are empty by default, so other studies don't change. CIBNOR will set --score_min L,0,-0.6.
• extractor_args is for M-bias trimming: read 2 alone shows 5% CHH methylation, against 0.8% for read 1.

Testing

The R code is not exercised by CI (see below). I extracted the module's methylKit script and ran it with methylKit 1.36.0 (Bioconductor 3.22, R 4.5.1) on synthetic Bismark cytosine reports: 12 libraries in 4 families × 3, 400 tiles × 25 CpGs. The data have real family-to-family variation and no RR-vs-SS difference, except 20 true regions with +30 points.

Fit	False positives (of 380 null)	True regions found (of 20)
12 libraries as independent (old behaviour)	162	20
4 families, MN + F (replicate_unit = family)	0	0

• Low power at 2 vs 2: true regions rank near the top (7 of the top 20 by q; median difference 26 points) but none reaches q < 0.05. That's the honest power of this design, and [study] CIBNOR2021_HEAT_RRSS (PRJNA690951): first reanalysis #14 will say so.
• A known methylKit quirk: when a region's estimated overdispersion falls below 1, methylKit resets it to 1. So 3 of 380 null regions still get q < 0.05, though none passes the 10-point meth_diff_cutoff. Documented in workflows/methylation/README.md.
• Default path unchanged: without replicate_unit, the DMR table is byte-identical to main's module on the same input, apart from the new column.
• Refusals triggered as expected: missing column; a unit spanning both arms; one unit per arm.
• CI checks run locally: the methylation demo (CI's command) passes; pytest -q: 270 passed; ruff check src tests scripts: clean.
• Nextflow wiring: a raw-mode run with replicate_unit set reaches DMR_METHYLKIT with the design sheet attached.
• Not covered by CI: no methylation module has a stub: block, so -stub-run executes real tools and the raw path can't be stub-tested. That's existing scaffolding, out of scope here.

Checklist

• pytest -q and ruff check src tests pass
• New logic has a test. Not in CI: the logic is R inside the Nextflow module, with no harness that executes it (as in Pool condition levels per arm, and fit DESeq2 on tanks with --replicate_unit #12). 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, so they can rebase after it merges

Study branches waiting on this

• study/PRJNA690951/CIBNOR2021_HEAT_RRSS (#14): set methylation.replicate_unit = 'family' and methylation.design_sheet to its samplesheet.tsv.

🤖 Generated with Claude Code