Add haplotype DMR (differential methylation region) calling

[YannVRB]

PR checklist

• This comment contains a description of changes (with reason).
• Module-level tests added (modules/nf-core/modkit/dmr, real, non-stub, uses a public nf-core/test-datasets fixture; see below).
• Follows pipeline conventions (new nf-core module vendored ahead of upstream merge; local modules follow the standard main.nf/environment.yml/meta.yml layout).
• nf-core pipelines lint (known 2.14.1 linter bug on this checkout, reproduces identically on the pre-existing modkit/pileup module, not something this PR introduces).
• End-to-end wiring verified locally against real phased dual-modification ONT data (private, not included in this PR, see Testing below).
• Usage/Output docs, CHANGELOG, README not yet updated — will follow up before merge.

What

Adds haplotype-vs-haplotype differential methylation region (DMR) calling, discussed with @ljwharbers: when --modkit_phased is set, lrsomatic now automatically runs a DMR comparison between a sample's two haplotypes at CpG islands, annotated with the nearest gene. Tumor/normal DMR is explicitly out of scope for this PR (deferred; cross-sample DMR stays in a downstream consumer pipeline).

New subworkflow subworkflows/local/dmr.nf:

MODKIT_PILEUP (hp1/hp2 bedMethyl)
-> HTSLIB_BGZIPTABIX (tabix-index each haplotype)
-> DMR_HAPLOTYPE_REGIONS (restrict to CpG islands covered in both haplotypes; new local module, ports the NCP downstream pipeline's dmr_calling.sh pattern)
-> MODKIT_DMR (nf-core modkit/dmr, vendored ahead of upstream merge -- nf-core/modules#13021)
-> DMR_NEAREST_GENE (nearest-gene annotation via bedtools closest; new local module)

Runs automatically whenever --modkit_phased is set; a new --skip_dmr (default false) turns it off independently.

New reference data

dmr_cpg_islands_bed / dmr_gencode_gene_bed added to conf/igenomes.config for both GRCh38 and CHM13. Currently point at the IntGenomicsLab/test-datasets fork branch (YannVRB/test-datasets@add-dmr-reference-fixtures) pending that PR's merge -- TODO comments mark where to re-point once it lands.

Testing

• Module level (modules/nf-core/modkit/dmr/tests/): real, non-stub nf-test against a public bedMethyl fixture in nf-core/test-datasets (Add modkit dmr bedmethyl nf-core/test-datasets#2289), itself derived from an existing, already-merged nf-core fixture BAM. Passes.
• Pipeline level: none of lrsomatic's existing CI BAMs carry any MM/ML methylation tags, so --modkit_phased` has never been exercisable in CI. This PR was verified end-to-end (all four subworkflow steps, both tumor and normal, real gene annotations and per-haplotype methylation statistics) against real ONT SUP-basecalled dual-modification data locally, but that data is private and is not included in this PR, so there is currently no committed pipeline-level nf-test exercising the full real-data path -- CI coverage for this feature is presently limited to the module-level test above. Flagging this explicitly for reviewers rather than skipping the disclosure; happy to discuss options (e.g. sourcing a suitable public phased-methylation fixture) as a follow-up.

Bugs found and fixed via real end-to-end testing (not stubs)

• HTSLIB_BGZIPTABIX output filename collision between a sample's hp1/hp2 outputs --ext.prefixnow includesmeta.haplotype`.
• DMR_NEAREST_GENE's sort -k1,1V isn't supported by BusyBox sort (bundled in the bedtools:2.31.1 biocontainer, which ships a coreutils-minimal sort) -- switched to plain -k1,1 -k2,2n.

Follow-ups (not blocking, tracked separately)

• Report a separate finding: the pipeline's default modkit_args (--cpg --modified-bases 5mC) is incompatible with modkit dmr pair against real dual-modification (5mC+5hmC) SUP-basecalled data -- see maintainer report.
• Once Add modkit/dmr module nf-core/modules#13021 and its test-datasets fixture PR merge, re-point modules.json and conf/igenomes.config at the canonical locations.
• Source a public phased-methylation pipeline-level test fixture, or accept module-level coverage as sufficient for now.
• Usage/output docs, CHANGELOG, README updates.

[ljwharbers]

Thanks for this, @YannVRB. Haplotype DMR is a good fit for --modkit_phased, and the channel wiring holds up. I checked that the _hp1/_hp2.bed.gz selection matches MODKIT_PILEUP's phased output names, that the haplotype tag is added to and then stripped from meta, and that the join/multiMap afterwards is correct. I also downloaded all four reference BEDs. They use chr naming, and the GENCODE files are sorted lexicographically, which matches the sort -k1,1 -k2,2n in DMR_NEAREST_GENE. Thanks also for re-pointing the reference BEDs at IntGenomicsLab/test-datasets/main in a39317d. All four URLs resolve.

I'm requesting changes for the items below. The main problem is that the step runs by default whenever --modkit_phased is set, yet CI never runs it and it can crash on inputs that currently work.

Blocking

1. CI never ran any of the new code, and the vendored modkit/dmr test can't pass.
nf-test-changes found the modkit/dmr, htslib/bgziptabix, bedtools/closest and bedtools/intersect tests, but the PR run filters on --tag small, which drops all four. That's why CI is green: none of them ran. You closed nf-core/test-datasets#2289 in favour of generating the fixture in-test, which works upstream (nf-core/modules#13021 at d886465 now runs a setup { run("MODKIT_PILEUP") } on test.sorted.phased.bam). But the copy vendored here is still pinned at 9d93dfa, and it reads genomics/homo_sapiens/nanopore/bedmethyl/test_hp1.bed.gz, which returns 404. Please re-vendor from the current #13021 head and update the git_sha in modules.json. Since manascripts has just asked for changes upstream, the best time is after you've addressed that. Their point, that the log output isn't deterministic, matches the test.log md5 I was going to raise. Asserting process.out.log.size() == 1 and snapshotting the other outputs, as they suggest, fixes it in both places.

2. We can have a pipeline-level test now: a public phased, methylation-tagged fixture already exists.
genomics/homo_sapiens/nanopore/bam/test.sorted.phased.bam is on nf-core/test-datasets, and this repo already uses it in the MODKIT_PILEUP --phased --modified-bases 5mC 5hmC module test. Your updated upstream module test already shows the pattern: a setup { run("MODKIT_PILEUP") } on that BAM. A subworkflows/local/dmr.nf nf-test that does the same, with a small chr22 CpG-island BED and a gene BED, tagged small, would cover the whole chain in PR CI. I'd rather have that in this PR than as a follow-up.

3. With --modkit_phased, any genome other than GRCh38/CHM13 hits a hard failure.
getGenomeAttribute('dmr_cpg_islands_bed') returns null for a custom --fasta/--genome, and file(null, checkIfExists: true) then fails. Users who never asked for DMR get this, because the step is on by default. The two params also aren't declared in nextflow.config or nextflow_schema.json, so users can't point to their own BEDs. Suggested fix: declare dmr_cpg_islands_bed / dmr_gencode_gene_bed as params with schema entries, let an explicit value override the genome default, and skip DMR with a log.warn (or give a clear error) when neither is available.

4. Default modkit_args vs DMR on dual-modification data.
Your follow-up notes say the default --cpg --modified-bases 5mC is incompatible with modkit dmr pair on real 5mC+5hmC SUP data. If that's true, most --modkit_phased runs on standard ONT data would now fail, or silently produce something wrong, unless the user knows to pass --skip_dmr. What exactly goes wrong, and which modkit_args did your end-to-end run use? I think this should be fixed in this PR, for example by deriving the MODKIT_DMR args from the pileup args. The alternative is to make DMR opt-in until it's fixed.

Should fix

• modules.json still points at YannVRB/modules@add-modkit-dmr. Vendoring ahead of nf-core/modules#13021 is fine. Please just keep the git_sha in sync with that PR (see 1) and switch it to nf-core once it merges.
• The vendored bedtools/closest and bedtools/intersect are never used. They add about 790 lines, including tests and snapshots. Either use BEDTOOLS_CLOSEST in place of the local DMR_NEAREST_GENE, which is only sort plus bedtools closest, or remove both modules and their modules.json entries.
• DMR_HAPLOTYPE_REGIONS unpacks and sorts whole-genome bedMethyl twice (see the inline comment).
• Docs: docs/usage.md (--skip_dmr and the reference BEDs), docs/output.md (methylation/<type>/dmr/*.nearestGene.bed and its columns: 17 modkit dmr columns, then the gene BED6, then the distance), and a CHANGELOG entry.

Nits

• -D a in DMR_NEAREST_GENE (see inline comment).
• MODKIT_DMR uses the 0.6.1 biocontainer, while MODKIT_PILEUP uses the patched 0.6.4 image. That's harmless, since DMR only reads bedMethyl, but a one-line comment would stop someone from "aligning" them later.
• MODKIT_DMR gets --ref without a .fai. Please check whether modkit re-indexes the 3 GB FASTA in every task. If it does, pass fai through.
• The GENCODE BEDs contain duplicate records at identical coordinates (e.g. hsa-mir-1253 / MIR1253) and some bare ENSG… names. With -t first, which of the duplicates is reported depends on file order.
• The checklist mentions a "2.14.1 linter bug", but this repo is on template 4.1.0 and the CI nf-core lint job passed. You can probably tick that box.
• Consider publishing the MODKIT_DMR log, for debugging.

[ljwharbers]

For any genome that isn't GRCh38/CHM13 (custom --fasta/--genome), params.dmr_cpg_islands_bed is null here, so file(null, checkIfExists: true) fails the run for every --modkit_phased user. Can you declare both params in nextflow.config and the schema, let an explicit value override the genome default, and skip DMR with a log.warn (or raise a clear error) when neither is set?

[ljwharbers]

On a whole-genome run this zcat | awk | sort -us tens of millions of CpG rows per haplotype. As the PR notes, the container ships BusyBox sort, which sorts in memory, and the process is process_single (6 GB). bedtools intersect -u doesn't need sorted input, and the staged .tbi files are never used. tabix -R ${cpg_islands_bed} hp1.bed.gz (islands sorted), or bedtools intersect -sorted -g genome.txt -u -a islands -b hp1.bed.gz, would stream the data instead of loading it all into memory.

Also, "covered" currently means at least one bedMethyl row in each haplotype, at any depth. Should there be a minimum valid coverage (column 10)?

[ljwharbers]

The DMR regions come from unstranded CpG islands, so modkit's column 6 (strand) is ., and -D a doesn't do anything meaningful. -D b, which reports distance relative to the gene's strand so negative means upstream of the gene, is probably what you want. Please document the sign convention in output.md.

[ljwharbers]

bedtools/closest and bedtools/intersect are installed but nothing includes them. Either use BEDTOOLS_CLOSEST in place of DMR_NEAREST_GENE, or drop both modules and these entries.

[ljwharbers]

This fixture returns 404: nf-core/test-datasets#2289 was closed, and upstream #13021 now generates the bedMethyl in-test with setup { run("MODKIT_PILEUP") } on test.sorted.phased.bam. This vendored copy (9d93dfa) predates that change. Please re-vendor from the current #13021 head. When you do, also pick up the fix for the non-deterministic log snapshot that manascripts asked for upstream. PR CI didn't catch this because the --tag small filter skips this test.