[YannVRB]

PR checklist

• This comment contains a description of changes (with reason).
• Module-level tests added (modules/nf-core/modkit/dmr, real, non-stub, chains MODKIT_PILEUP on a public nf-core/test-datasets BAM fixture; see below).
• Pipeline/subworkflow-level test added (subworkflows/local/tests/dmr.nf.test, tag small, real non-stub chain; see below).
• Follows pipeline conventions (new nf-core module vendored ahead of upstream merge, now re-pointed at nf-core/modules proper; local modules follow the standard main.nf/environment.yml/meta.yml layout).
• nf-core pipelines lint and pre-commit both pass in CI (the module-provenance/lint issue below was fixed, not a persistent linter bug).
• 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 and CHANGELOG updated.

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, queries via tabix -R rather than unpacking the whole genome, ports the NCP downstream pipeline's dmr_calling.sh pattern)
-> MODKIT_DMR (nf-core modkit/dmr, now pointed at the merged nf-core/modules#13021)
-> DMR_NEAREST_GENE (nearest-gene annotation via bedtools closest -D b; new local module)

Runs automatically whenever --modkit_phased is set; --skip_dmr (default false) turns it off independently. On any genome other than GRCh38/CHM13, DMR is skipped with a log.warn unless --dmr_cpg_islands_bed/--dmr_gencode_gene_bed are supplied explicitly.

Real bug fixed, not just worked around: the pipeline's default --modkit_args (--cpg --modified-bases 5mC) crashes modkit dmr pair on real dual-modification (5mC+5hmC) SUP-basecalled ONT data, since reads carrying an unrequested modification code land in N_other and fail dmr pair's valid_coverage == N_canonical + N_modified consistency check. Fixed by deriving the actual pileup args in conf/modules.config: when --modkit_phased is set and DMR will run, 5hmC is appended automatically if --modified-bases names 5mC without it, regardless of what --modkit_args the user supplied otherwise.
Full writeup in the PR comments.

New reference data

dmr_cpg_islands_bed / dmr_gencode_gene_bed are declared in conf/igenomes.config for both GRCh38 and CHM13, pointed at IntGenomicsLab/test-datasets/main (merged via IntGenomicsLab/test-datasets#4).

Testing

• Module level (modules/nf-core/modkit/dmr/tests/): real, non-stub nf-test that chains MODKIT_PILEUP (--phased --modified-bases 5mC 5hmC) on the public test.sorted.phased.bam nf-core/test-datasets fixture into MODKIT_DMR, instead of a static fixture. Passes.
• Subworkflow level (subworkflows/local/tests/dmr.nf.test, tag small so it actually runs in PR CI): real, non-stub nf-test chaining the same public phased BAM fixture through the full MODKIT_PILEUP -> HTSLIB_BGZIPTABIX -> DMR_HAPLOTYPE_REGIONS -> MODKIT_DMR -> DMR_NEAREST_GENE path against a small synthetic chr22 CpG-island/gene pair. Verifies real (not fabricated) per-haplotype methylation counts/fractions/deltas, correctly annotated with the test's synthetic gene.
• End-to-end: all four subworkflow steps, both tumor and normal, real gene annotations and per-haplotype methylation statistics, verified locally against real ONT SUP-basecalled dual-modification data.
  That data is private and not included in this PR (see Testing disclosure note below), but the subworkflow-level test above now gives public, CI-visible coverage of the same chain.

Private-data disclosure: the full real-data end-to-end run above was done against confidential data that is not, and will not be, committed to this repo or referenced by path/name anywhere in it; the public subworkflow-level test is a separate, independently-verified real (non-fabricated) run against public fixtures only.

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.prefix now includes meta.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.
• Any --genome other than GRCh38/CHM13 hard-failed under --modkit_phased, because getGenomeAttribute() unconditionally overwrote an explicit --dmr_cpg_islands_bed override with null. Fixed: both BEDs declared in schema + config, explicit value wins over the genome default, and DMR is skipped with a warning (not a crash) when neither is available.
• MODKIT_DMR got --ref with no .fai, so modkit dmr pair rebuilt one from scratch on every task for a multi-GB genome. Fixed upstream in Add modkit/dmr module nf-core/modules#13021 (now merged) and re-vendored.
• bedtools closest -D a -> -D b: distance sign should be relative to the gene's strand, not the unstranded DMR region's.
• Default --modkit_args incompatible with real dual-modification data (see "What" above) -- fixed by deriving DMR's pileup args instead of making DMR opt-in.
• Removed vendored-but-unused bedtools/closest/bedtools/intersect modules.
• MODKIT_DMR now uses the same patched modkit 0.6.4 image as MODKIT_PILEUP, for one consistent modkit binary across the chain.

Follow-ups (not blocking, tracked separately)

• IntGenomicsLab/test-datasets's dmr_gencode_gene_bed has a known data-quality issue (~46% bare ENSG... IDs instead of gene symbols, plus some duplicate-coordinate records). Not yet filed as a separate issue against that repo; deferred, out of scope for this PR.
• Source a public phased dual-modification ONT fixture for an even closer end-to-end CI test, if one becomes available; the current subworkflow-level test already gives real public-fixture coverage of the full chain.

[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.

[YannVRB]

What goes wrong: modkit dmr pair has its own internal consistency check on each bedMethyl row: valid_coverage must equal N_canonical + N_modified. It doesn't, whenever a read's modification call isn't one of the codes named in --modified-bases, that read gets bucketed into N_other instead, and dmr pair rejects the row outright:

encountered invalid bedMethyl record(s), invalid data, valid coverage (2) is not equal to the sum of canonical and modified counts (1), [BedMethylLine { ... raw_mod_code: Code('m') ... }]

Which modkit_args the end-to-end run used: the pipeline default, --cpg --modified-bases 5mC. Real ONT SUP-basecalled data calls both 5mC and 5hmC at the same CpG sites, so with only 5mC requested, every read with a 5hmC call lands in N_other, and MODKIT_DMR crashes on the first such position. This isn't PacBio-specific (the schema's existing N_other note was written with PacBio's conflicting-probability case in mind), it reproduces identically on dual-modification ONT data, since the mechanism (uncounted mod code -> N_other) is the same.

Fix, pushed in 515b522: rather than make DMR opt-in, conf/modules.config's MODKIT_PILEUP block now derives the actual pileup args instead of passing modkit_args straight through: when --modkit_phased is set and DMR isn't skipped, it appends 5hmC to --modified-bases if that list names 5mC without it, regardless of what modkit_args the user supplied otherwise. An unphased run, or a --skip_dmr run, is untouched; the user's modkit_args passes through exactly as given. Also updated modkit_args's help_text in nextflow_schema.json to document this.

This matches what subworkflows/local/tests/dmr.nf.test already does explicitly in its own local config (ext.args = '--phased --modified-bases 5mC 5hmC') — that test's real (non-fabricated) DMR output is itself evidence this is the right pair of codes for the production path, not just the test fixture.

[ljwharbers]

Thanks @YannVRB, this addresses all four points from the first review: modkit/dmr now points at the merged #13021, dmr.nf.test runs in PR CI, the params are declared and skip cleanly on other genomes, and the 5hmC/N_other problem is handled.

While checking it on Mindwell (Apptainer, real test-profile run with --modkit_phased), I found two more bugs. I pushed fixes for them directly to the branch, plus a few smaller fixes:

1. The genome-default BEDs were never applied (843ec83). nextflow.config declares dmr_cpg_islands_bed/dmr_gencode_gene_bed as null, and Nextflow silently ignores a later params.x = ... in the workflow for a param a config already defines (checked on 26.04.3). So --genome GRCh38 --modkit_phased printed "no default for --genome 'GRCh38'" and skipped DMR unless both BEDs were passed by hand. The other getGenomeAttribute params only work because they aren't declared. The BEDs are now resolved into local variables.
2. An empty regions BED killed the run (ed5e7e2). A BAM without MM/ML tags gives empty pileups, so no island is covered in both haplotypes. modkit dmr pair then exits with "zero non-comment lines in regions", and the whole pipeline stops. Those samples are now skipped with a warning, and there is a second dmr.nf.test case for this.
3. Docs about the 5hmC append (1e2bbda). usage.md said the published modkit_pileup files keep --modkit_args exactly as given. They don't: DMR reads the same MODKIT_PILEUP task that publishes them, so they gain a 5hmC row per CpG. I corrected usage.md and the schema help, and moved the CHANGELOG entry to Changed. The config append is now also skipped when DMR itself will be skipped for lack of BEDs. I also noted in output.md that --base C scores 5mC and 5hmC together.
4. The vendored modkit/dmr module test couldn't run here (8ebea39). Its setup passed upstream's 3-input MODKIT_PILEUP call to our 4-input patched pileup, and the snapshot recorded modkit 0.6.1 while the patch runs 0.6.4. I fixed the setup, re-recorded the snapshot from a real run (3/3 pass) and regenerated the patch.
5. Smaller fixes: the DMR_HAPLOTYPE_REGIONS coverage threshold now reads ext.min_valid_coverage instead of overloading ext.args; the skip warning no longer prints --genome 'null'; and current dev is merged in.

Checks run locally:

• default pipeline snapshot: unchanged after the dev merge;
• dmr.nf.test: 2/2 pass;
• modkit/dmr module test: 3/3 pass;
• phased test-profile run: completes;
• pre-commit and nf-core pipelines lint: clean.

The only modules lint failure is oras_singularity_tag, which dev's modkit/pileup shares. I've left the linting.yml prek→pre-commit change as it is.

From my side this is good to go once CI is green on the new head.