ljwharbers · Oct 27, 2026 · Sep 23, 2026 · Sep 23, 2026 · Sep 23, 2026 · Sep 23, 2026
diff --git a/.github/workflows/linting.yml b/.github/workflows/linting.yml
@@ -16,8 +16,19 @@ jobs:
       - name: Install Nextflow
         uses: nf-core/setup-nextflow@893c28b667aedeba26e37f296d260ccc5bc4d914 # v3
 
-      - name: Run prek
-        uses: j178/prek-action@5337cb91e0fa35a7ff31b9ca345126d8bbbcdf16 # v2
+      - uses: actions/setup-python@5fda3b95a4ea91299a34e894583c3862153e4b97 # v7
+        with:
+          python-version: "3.13"
+          architecture: "x64"
+
+      - name: Setup uv
+        uses: astral-sh/setup-uv@c771a70e6277c0a99b617c7a806ffedaca235ff9 # v9.0.0
+
+      - name: Install pre-commit
+        run: uv tool install pre-commit
+
+      - name: Run pre-commit
+        run: pre-commit run --all-files --show-diff-on-failure
 
   nf-core:
     runs-on: ubuntu-latest

diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -7,6 +7,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### `Added`
 
+- [#204](https://github.com/IntGenomicsLab/lrsomatic/pull/204) - Added haplotype differential methylation region (DMR) calling: with `--modkit_phased` set, `modkit dmr pair` compares methylation between the `_hp1`/`_hp2` bedMethyl pileups over CpG islands covered in both haplotypes, and each DMR call is annotated with its nearest gene via `bedtools closest`. Runs automatically whenever `--modkit_phased` is set; turn it off with `--skip_dmr`. `--dmr_cpg_islands_bed`/`--dmr_gencode_gene_bed` have built-in defaults for `--genome GRCh38`/`CHM13` and must be supplied explicitly for any other genome, or DMR is skipped with a warning rather than failing the run. Tumour/normal (cross-sample) DMR is explicitly out of scope for this PR (@YannVRB).
+- [#204](https://github.com/IntGenomicsLab/lrsomatic/pull/204) - Added a real, public subworkflow-level nf-test for the `DMR` subworkflow (tag `small`), chaining `MODKIT_PILEUP` on the public phased, methylation-tagged `test.sorted.phased.bam` fixture through the full DMR chain against a small synthetic CpG-island/gene region (@YannVRB).
 - [#205](https://github.com/IntGenomicsLab/lrsomatic/pull/205) - Added ecDNA and focal amplification analysis: CoRAL reconstructs amplicon structures from each tumour BAM, seeded by ASCAT's copy-number calls, and AmpliconClassifier labels each amplicon as ecDNA, BFB or linear. Runs by default on tumour samples with ASCAT calls; turn it off with `--skip_coral`, or keep reconstruction without classification with `--skip_ampliconclassifier`. CoRAL's optimisation defaults to the open-source SCIP solver rather than Gurobi: across the 1102 solver logs from the previous standalone cohort the largest model was 620 rows by 438 columns and the slowest solve took 0.49 s, so a licensed solver buys nothing at this scale, though `--coral_solver gurobi_direct` with `--gurobi_license` remains available. AmpliconClassifier's AmpliconArchitect data repository is downloaded automatically for GRCh38, checked against a pinned `--aa_data_repo_md5`, and can be supplied with `--aa_data_repo` (the reference directory or its parent); CHM13 has no published repository, so without `--aa_data_repo` the classifier is skipped with a warning (@robert-a-forsyth).
 - [#205](https://github.com/IntGenomicsLab/lrsomatic/pull/205) - Added `ascat_to_coral_bed.py` and the `ASCAT_TO_CORAL_BED` module, converting ASCAT's `cnvs.txt` into the headerless BED CoRAL seeds from. Contigs are respelled to match the reference, since ASCAT writes `1` where the BAM may say `chr1` and CoRAL builds its chromosome sizes from the BAM header; zero-length and inverted segments are dropped, which CoRAL's own parser would otherwise reject (@robert-a-forsyth).
 - [#205](https://github.com/IntGenomicsLab/lrsomatic/pull/205) - Added stub nf-tests for the `ECDNA` subworkflow and `ASCAT_TO_CORAL_BED` (tag `small`), covering the seeded path, the empty-seed branch, the opt-in cycle re-extraction with plotting, the `--skip_ampliconclassifier` path and a run with no data repository, plus a stub test for `AADATAREPO_DOWNLOAD` (@robert-a-forsyth).
@@ -28,6 +30,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### `Changed`
 
+- [#204](https://github.com/IntGenomicsLab/lrsomatic/pull/204) - With `--modkit_phased` set and haplotype DMR running, `5hmC` is appended to `--modified-bases` whenever `--modkit_args` names `5mC` without it, so the published `modkit_pileup/*_{hp1,hp2,combined}.bed.gz` now carry a 5hmC row next to each 5mC row. `modkit dmr pair` rejects bedMethyl rows whose valid coverage includes uncounted (`N_other`) calls, which ONT SUP data calling both 5mC and 5hmC produces under the `--cpg --modified-bases 5mC` default. Runs with `--skip_dmr`, without `--modkit_phased`, or without DMR references keep `--modkit_args` exactly as given (@YannVRB).
 - [#205](https://github.com/IntGenomicsLab/lrsomatic/pull/205) - Moved the fiber-seq steps into a `FIBERSEQ` subworkflow and the AmpliconArchitect MD5 default into `aaDataRepoMd5()`, because with dev merged the `LRSOMATIC` workflow body exceeded the 65,535-character limit on a compiled string and failed to compile; behaviour is unchanged (@robert-a-forsyth).
 - [#205](https://github.com/IntGenomicsLab/lrsomatic/pull/205) - `--download_sigprofiler_genome` now defaults to `true`, so a run without `--sigprofiler_genome_dir` installs the SigProfiler genome instead of stopping at launch; a given directory still wins (@robert-a-forsyth).
 - [#205](https://github.com/IntGenomicsLab/lrsomatic/pull/205) - Every small-variant path now rejoins each arm to one record per position (`bcftools norm -m +any`) before Longphase, which cannot take two records at one position; single-caller runs previously passed an SNV and an indel at the same base (@robert-a-forsyth).

diff --git a/conf/igenomes.config b/conf/igenomes.config
@@ -24,6 +24,8 @@ params.genomes = [
         vep_species     : "homo_sapiens",
         savana_contigs   : "https://raw.githubusercontent.com/cortes-ciriano-lab/savana/main/example/contigs.chr.hg38.txt",
         savana_g1000_vcf : "1000g_hg38",
+        dmr_cpg_islands_bed  : "https://raw.githubusercontent.com/IntGenomicsLab/test-datasets/main/references/dmr/GRCh38.cpg_islands.bed.gz",
+        dmr_gencode_gene_bed : "https://raw.githubusercontent.com/IntGenomicsLab/test-datasets/main/references/dmr/GRCh38.gencode_gene.bed.gz",
         coral_ref        : "hg38",
         ac_ref           : "GRCh38",
         // Plain build, not GRCh38_indexed: the extra BWA index is for AmpliconArchitect's
@@ -64,6 +66,8 @@ params.genomes = [
         vep_species     : "homo_sapiens_gca009914755v4",
         savana_contigs   : "https://raw.githubusercontent.com/IntGenomicsLab/test-datasets/main/references/savana/contigs.chr.chm13.txt",
         savana_g1000_vcf : "1000g_t2t",
+        dmr_cpg_islands_bed  : "https://raw.githubusercontent.com/IntGenomicsLab/test-datasets/main/references/dmr/CHM13.cpg_islands.bed.gz",
+        dmr_gencode_gene_bed : "https://raw.githubusercontent.com/IntGenomicsLab/test-datasets/main/references/dmr/CHM13.gencode_gene.bed.gz",
         coral_ref        : "t2t",
         ac_ref           : "CHM13",
         // No published AA data repo for CHM13; supply one with --aa_data_repo

diff --git a/conf/modules.config b/conf/modules.config
@@ -397,10 +397,28 @@ process {
 
     withName: '.*:MODKIT_PILEUP' {
         ext.args = {
+                // A bare or empty `--modkit_args` reaches here as true (Nextflow 25) or 'true' (26): treat both as no arguments.
+                // The documented form is `--modkit_args=''`.
+                def modkit_args = (params.modkit_args == null || params.modkit_args in [true, 'true']) ? '' : params.modkit_args
+
+                // modkit dmr pair rejects rows with N_other > 0 and SUP models also call 5hmC, so add
+                // 5hmC whenever DMR will run (same gate as workflows/lrsomatic.nf).
+                def genome_attrs = params.genomes?.get(params.genome) ?: [:]
+                def dmr_beds = (params.dmr_cpg_islands_bed ?: genome_attrs.dmr_cpg_islands_bed) &&
+                    (params.dmr_gencode_gene_bed ?: genome_attrs.dmr_gencode_gene_bed)
+                if (params.modkit_phased && !params.skip_dmr && dmr_beds) {
+                    def match = (modkit_args =~ /--modified-bases(?:\s+[A-Za-z0-9]+)+/)
+                    if (match.find()) {
+                        def modified_bases_arg = match.group()
+                        def codes = modified_bases_arg.replaceFirst('--modified-bases', '').trim().split(/\s+/) as List
+                        if (codes.contains('5mC') && !codes.contains('5hmC')) {
+                            modkit_args = modkit_args.replace(modified_bases_arg, "${modified_bases_arg} 5hmC")
+                        }
+                    }
+                }
+
                 [
-                    // A bare or empty `--modkit_args` reaches here as true (Nextflow 25) or 'true' (26): treat both as no arguments.
-                    // The documented form is `--modkit_args=''`.
-                    (params.modkit_args == null || params.modkit_args in [true, 'true']) ? '' : params.modkit_args,
+                    modkit_args,
                     params.modkit_phased ? '--phased' : ''
                 ].join(' ').trim()
             }
@@ -411,6 +429,38 @@ process {
         ]
     }
 
+    //
+    // SUBWORKFLOW: DMR -- only the nearest-gene BED and the modkit dmr log are published
+    //
+    withName: '.*:DMR:HTSLIB_BGZIPTABIX' {
+        // Haplotype in the filename, or hp1/hp2 collide once meta.haplotype is dropped
+        ext.prefix = { "${meta.id}_${meta.haplotype}" }
+        ext.args2 = '-p bed'
+        publishDir = [
+            enabled: false
+        ]
+    }
+    withName: '.*:DMR:DMR_HAPLOTYPE_REGIONS' {
+        publishDir = [
+            enabled: false
+        ]
+    }
+    withName: '.*:DMR:MODKIT_DMR' {
+        ext.args = '--base C'
+        publishDir = [
+            path: { "${params.outdir}/${meta.id}/methylation/${meta.type}/dmr" },
+            mode: params.publish_dir_mode,
+            saveAs: { filename -> filename.endsWith('.log') ? filename : null }
+        ]
+    }
+    withName: '.*:DMR:DMR_NEAREST_GENE' {
+        publishDir = [
+            path: { "${params.outdir}/${meta.id}/methylation/${meta.type}/dmr" },
+            mode: params.publish_dir_mode,
+            saveAs: { filename -> filename.equals('versions.yml') ? null : filename }
+        ]
+    }
+
     withName: '.*:FIBERTOOLSRS_PREDICTM6A' {
         ext.args = {
             [

diff --git a/docs/output.md b/docs/output.md
@@ -264,6 +264,27 @@ QC outputs are placed under `tumor/` for all samples, and additionally under `no
 | `{tumor,normal}/modkit_pileup/sample_{hp1,hp2}.bed.gz` | With `--modkit_phased`: bedMethyl tables restricted to reads carrying haplotype tag `HP:1` or `HP:2` from the Longphase-haplotagged BAM.                                                                        |
 | `{tumor,normal}/modkit_pileup/sample_combined.bed.gz`  | With `--modkit_phased`: bedMethyl table over all reads, including untagged ones (equivalent to the unphased default output). There is no separate file for untagged reads; they only contribute to `_combined`. |
 
+With `--modkit_phased` set (and `--skip_dmr` not set), haplotype differential methylation region (DMR) calling also runs, comparing the `_hp1`/`_hp2` bedMethyl pileups over CpG islands covered in both haplotypes and annotating each result with its nearest gene:
+
+```
+├── methylation
+│   ├── tumor
+│   │   ├── modkit_pileup
+│   │   │   └── ...
+│   │   └── dmr
+│   │       ├── sample.log                   # modkit dmr pair's debug log
+│   │       └── sample.nearestGene.bed       # final DMR deliverable
+│   ├── normal                               # paired samples only
+│   │   └── ...                              # same layout as tumor
+```
+
+| File                                        | Description                                                                                                                                                                                                                                                                                                                                                                                                                                                              |
+| ------------------------------------------- | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ |
+| `{tumor,normal}/dmr/sample.log`             | `modkit dmr pair`'s own debug log, published for troubleshooting (e.g. the `invalid-bedmethyl-data` class of failure) without needing `-resume`. The intermediate raw `modkit dmr pair` BED is not published, since this file supersedes it.                                                                                                                                                                                                                             |
+| `{tumor,normal}/dmr/sample.nearestGene.bed` | Each `modkit dmr pair` DMR call (17 columns — chrom/start/end, haplotype-vs-haplotype methylation counts/fractions and the difference score; see `modkit dmr pair --help` for the full column list), joined via `bedtools closest -D b -t first` to its nearest gene from `--dmr_gencode_gene_bed` (that gene's own BED6: chrom/start/end/name/score/strand), followed by the signed distance to that gene (relative to the gene's own strand; negative means upstream). |
+
+DMR runs `modkit dmr pair --base C`, which counts every modification on C as modified. The pileup feeding it requests both 5mC and 5hmC (see `--modkit_args` in the usage docs), so the methylation fractions and the difference score describe 5mC and 5hmC combined, not 5mC alone.
+
 </details>
 
 ### `variants`

diff --git a/docs/usage.md b/docs/usage.md
@@ -193,6 +193,7 @@ Padfoot bundles gene and repeat annotations for `hg38` and `mm10` only. With `--
 | `--skip_normalfiber`   | A boolean to skip fibertools processing for the normal sample. Default = `false`                                                                                                                                                                                        |
 | `--skip_modcall`       | A boolean to skip Longphase `modcall`, the 5mC base-modification calling whose VCF is used as extra evidence during phasing. Unrelated to the modkit pileup (see `--skip_modkit`). Default = `false`                                                                    |
 | `--skip_modkit`        | A boolean to skip the modkit pileup step. Default = `false`                                                                                                                                                                                                             |
+| `--skip_dmr`           | A boolean to skip haplotype differential methylation region (DMR) calling. Only takes effect when `--modkit_phased` is also set; has no effect otherwise. Default = `false`                                                                                             |
 | `--skip_whatshapstats` | A boolean to skip WhatsHap phasing statistics. Default = `false`                                                                                                                                                                                                        |
 | `--skip_signatures`    | A boolean to skip mutational signature analysis (SigProfilerMatrixGenerator + SigProfilerAssignment). Default = `false`                                                                                                                                                 |
 | `--skip_report`        | A boolean to skip the final per-sample HTML report. Default = `false`                                                                                                                                                                                                   |
@@ -206,6 +207,18 @@ Padfoot bundles gene and repeat annotations for `hg38` and `mm10` only. With `--
 
 The pileup runs a patched modkit 0.6.4 image (`ghcr.io/ljwharbers/modkit`, `linux/amd64` only, built from [ljwharbers/modkit@pacbio-conflict-fix](https://github.com/ljwharbers/modkit/tree/pacbio-conflict-fix)): stock modkit drops PacBio HiFi reads whose 5mC and 5hmC probabilities sum above 1 and returns empty `--cpg` pileups on them ([nanoporetech/modkit#612](https://github.com/nanoporetech/modkit/issues/612)). Conda and `arm64` are not supported for this step; use `--skip_modkit` there.
 
+#### DMR options:
+
+With `--modkit_phased` set, haplotype differential methylation region (DMR) calling runs automatically unless `--skip_dmr` is given: `modkit dmr pair` compares methylation between the `_hp1`/`_hp2` bedMethyl pileups over CpG islands covered in both haplotypes, and the result is annotated with its nearest gene.
+
+| Parameter                | Description                                                                                                                                                                                                                            |
+| ------------------------ | -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
+| `--skip_dmr`             | Skip DMR calling. Only takes effect when `--modkit_phased` is also set. Default = `false`                                                                                                                                              |
+| `--dmr_cpg_islands_bed`  | CpG-islands BED to restrict DMR calling to. Has a built-in default for `--genome GRCh38`/`CHM13`; required as an explicit override for any other genome. Default = `null` (built-in default for GRCh38/CHM13)                          |
+| `--dmr_gencode_gene_bed` | Gene-model BED used to annotate each DMR call with its nearest gene. Has a built-in default for `--genome GRCh38`/`CHM13`; required as an explicit override for any other genome. Default = `null` (built-in default for GRCh38/CHM13) |
+
+On any other genome, DMR is skipped with a `log.warn` (not a hard failure) unless both BEDs are supplied explicitly. `--modkit_args`'s default only requests 5mC calls, but ONT SUP-basecalled data also calls 5hmC at the same sites, and `modkit dmr pair` rejects bedMethyl rows whose valid coverage includes calls of an unrequested modification code. So whenever DMR will run (`--modkit_phased` set, `--skip_dmr` not set, both BEDs available), `5hmC` is appended to `--modified-bases` if it names `5mC` without it. DMR reads the same pileup that is published, so the published `modkit_pileup/*_{hp1,hp2,combined}.bed.gz` then contain a 5hmC row next to each 5mC row. Use `--skip_dmr` to keep `--modkit_args` exactly as given.
+
 #### LONGPHASE options:
 
 | Parameter                       | Description                                                                        |