Add modkit/dmr module

modules/nf-core/modkit/dmr/environment.yml

@@ -0,0 +1,7 @@
+---
+# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/environment-schema.json
+channels:
+  - conda-forge
+  - bioconda
+dependencies:
+  - ont-modkit=0.6.1

modules/nf-core/modkit/dmr/main.nf

@@ -0,0 +1,53 @@
+process MODKIT_DMR {
+    tag "${meta.id}"
+    label 'process_medium'
+
+    conda "${moduleDir}/environment.yml"
+    container "${workflow.containerEngine in ['singularity', 'apptainer'] && !task.ext.singularity_pull_docker_container
+        ? 'https://depot.galaxyproject.org/singularity/ont-modkit:0.6.1--hcdda2d0_0'
+        : 'quay.io/biocontainers/ont-modkit:0.6.1--hcdda2d0_0'}"
+
+    input:
+    tuple val(meta), path(bedmethyl_a), path(bedmethyl_a_tbi)
+    tuple val(meta2), path(bedmethyl_b), path(bedmethyl_b_tbi)
+    tuple val(meta3), path(regions_bed)
+    tuple val(meta4), path(fasta), path(fai)
+
+    output:
+    tuple val(meta), path("*.bed"), emit: bed
+    tuple val(meta), path("*.log"), emit: log
+    tuple val("${task.process}"), val('modkit'), eval("modkit --version | sed 's/modkit //'"), emit: versions_modkit, topic: versions
+
+    when:
+    task.ext.when == null || task.ext.when
+
+    script:
+    def args = task.ext.args ?: ''
+    def prefix = task.ext.prefix ?: "${meta.id}"
+    def regions = regions_bed ? "-r ${regions_bed}" : ''
+    // fai is never referenced directly -- modkit dmr pair has no --fai flag and looks for
+    // <fasta>.fai next to --ref itself. Declaring it as an input only stages it alongside
+    // fasta so modkit finds and reuses it instead of building one from scratch, which for a
+    // multi-GB genome is a real, repeated cost across every task invocation. Confirmed by
+    // direct testing: with no .fai present, modkit dmr pair creates one; with one already
+    // staged next to fasta, it's read (unmodified mtime) rather than rebuilt.
+    """
+    modkit \\
+        dmr pair \\
+        -a ${bedmethyl_a} \\
+        -b ${bedmethyl_b} \\
+        ${regions} \\
+        --ref ${fasta} \\
+        -o ${prefix}.bed \\
+        --log-filepath ${prefix}.log \\
+        -t ${task.cpus} \\
+        ${args}
+    """
+
+    stub:
+    def prefix = task.ext.prefix ?: "${meta.id}"
+    """
+    touch ${prefix}.bed
+    touch ${prefix}.log
+    """
+}

modules/nf-core/modkit/dmr/meta.yml

@@ -0,0 +1,140 @@
+name: modkit_dmr
+description: Compare regions between a pair of bedMethyl samples (e.g. two haplotypes,
+  or tumor and normal) and report differentially methylated regions
+keywords:
+  - methylation
+  - ont
+  - long-read
+  - differential methylation
+tools:
+  - "modkit":
+      description: A bioinformatics tool for working with modified bases in Oxford
+        Nanopore sequencing data
+      homepage: https://github.com/nanoporetech/modkit
+      documentation: https://github.com/nanoporetech/modkit
+      tool_dev_url: https://github.com/nanoporetech/modkit
+      licence: ["Oxford Nanopore Technologies PLC. Public License Version 1.0"]
+      identifier: ""
+input:
+  - - meta:
+        type: map
+        description: |
+          Groovy Map containing sample information for the first (control) input
+          e.g. `[ id:'test' ]`
+    - bedmethyl_a:
+        type: file
+        description: Bgzipped bedMethyl file for the first (control) sample, as
+          produced by modkit pileup with --bgzf
+        pattern: "*.bed.gz"
+        ontologies:
+          - edam: "http://edamontology.org/format_3003" # BED
+          - edam: "http://edamontology.org/format_3989" # GZIP format
+    - bedmethyl_a_tbi:
+        type: file
+        description: Tabix index for bedmethyl_a. Must have the same basename with
+          a .tbi extension next to bedmethyl_a
+        pattern: "*.bed.gz.tbi"
+        ontologies:
+          - edam: "http://edamontology.org/format_3700" # Tabix index format
+  - - meta2:
+        type: map
+        description: |
+          Groovy Map containing sample information for the second (experimental)
+          input e.g. `[ id:'test' ]`
+    - bedmethyl_b:
+        type: file
+        description: Bgzipped bedMethyl file for the second (experimental) sample,
+          as produced by modkit pileup with --bgzf
+        pattern: "*.bed.gz"
+        ontologies:
+          - edam: "http://edamontology.org/format_3003" # BED
+          - edam: "http://edamontology.org/format_3989" # GZIP format
+    - bedmethyl_b_tbi:
+        type: file
+        description: Tabix index for bedmethyl_b. Must have the same basename with
+          a .tbi extension next to bedmethyl_b
+        pattern: "*.bed.gz.tbi"
+        ontologies:
+          - edam: "http://edamontology.org/format_3700" # Tabix index format
+  - - meta3:
+        type: map
+        description: |
+          Groovy Map containing regions information
+          e.g. `[ id:'regions' ]`
+    - regions_bed:
+        type: file
+        description: Optional BED file of regions over which to compare methylation
+          levels. When omitted, methylation levels are compared at each site instead
+        pattern: "*.bed"
+        ontologies:
+          - edam: "http://edamontology.org/format_3003" # BED
+  - - meta4:
+        type: map
+        description: |
+          Groovy Map containing reference information
+          e.g. `[ id:'hg38' ]`
+    - fasta:
+        type: file
+        description: Reference sequence in FASTA format, used during pileup/alignment
+        pattern: "*.{fa,fasta}"
+        ontologies:
+          - edam: "http://edamontology.org/format_1929" # FASTA
+    - fai:
+        type: file
+        description: |
+          Optional FASTA index for fasta. modkit dmr pair has no --fai flag and looks for
+          <fasta>.fai next to --ref itself; passing one here only stages it alongside fasta
+          so modkit reuses it instead of building one from scratch on every invocation, which
+          for a multi-GB genome is a real, repeated cost. Pass `[]` to omit.
+        pattern: "*.fai"
+        ontologies: []
+output:
+  bed:
+    - - meta:
+          type: map
+          description: |
+            Groovy Map containing sample information
+            e.g. `[ id:'test' ]`
+      - "*.bed":
+          type: file
+          description: BED file with the score column indicating the magnitude of
+            the difference in methylation between the two samples
+          pattern: "*.bed"
+          ontologies:
+            - edam: "http://edamontology.org/format_3003" # BED
+  log:
+    - - meta:
+          type: map
+          description: |
+            Groovy Map containing sample information
+            e.g. `[ id:'test' ]`
+      - "*.log":
+          type: file
+          description: File for debug logs to be written to
+          pattern: "*.log"
+          ontologies: []
+  versions_modkit:
+    - - ${task.process}:
+          type: string
+          description: The name of the process
+      - modkit:
+          type: string
+          description: The name of the tool
+      - modkit --version | sed 's/modkit //':
+          type: eval
+          description: The expression to obtain the version of the tool
+topics:
+  versions:
+    - - ${task.process}:
+          type: string
+          description: The name of the process
+      - modkit:
+          type: string
+          description: The name of the tool
+      - modkit --version | sed 's/modkit //':
+          type: eval
+          description: The expression to obtain the version of the tool
+authors:
+  - "@YannVRB"
+maintainers:
+  - "@YannVRB"

modules/nf-core/modkit/dmr/tests/main.nf.test

@@ -0,0 +1,194 @@
+nextflow_process {
+
+    name "Test Process MODKIT_DMR"
+    script "../main.nf"
+    tag "modules"
+    tag "modules_nfcore"
+    tag "modkit"
+    tag "modkit/dmr"
+    tag "modkit/pileup"
+    tag "htslib/bgziptabix"
+    process "MODKIT_DMR"
+    config "./nextflow.config"
+
+    setup {
+        run("MODKIT_PILEUP") {
+            script "../../pileup/main.nf"
+            process {
+                """
+                input[0] = [
+                    [ id: 'test' ],
+                    file(params.modules_testdata_base_path + 'genomics/homo_sapiens/nanopore/bam/test.sorted.phased.bam', checkIfExists: true),
+                    file(params.modules_testdata_base_path + 'genomics/homo_sapiens/nanopore/bam/test.sorted.phased.bam.bai', checkIfExists: true)
+                ]
+                input[1] = [
+                    [ id: 'test_ref' ],
+                    file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta', checkIfExists: true),
+                    file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta.fai', checkIfExists: true)
+                ]
+                input[2] = [[],[]]
+                """
+            }
+        }
+
+        run("HTSLIB_BGZIPTABIX") {
+            script "../../../htslib/bgziptabix/main.nf"
+            process {
+                """
+                // MODKIT_PILEUP --phased emits one [meta, [combined, hp1, hp2]] tuple per sample
+                // ("bedgz_files" is a single Path rather than a List when only one file matches
+                // the module's output glob -- true for -stub, since modkit pileup's stub script
+                // never splits by haplotype -- so normalize via [x].flatten() first). Split into
+                // one [meta, bedgz] tuple per haplotype so each can be tabix-indexed (modkit
+                // pileup --bgzf already bgzips, but never indexes) and fed to MODKIT_DMR as its
+                // own bedmethyl_a/bedmethyl_b input, exactly like real --modkit_phased usage.
+                input[0] = MODKIT_PILEUP.out.bedgz
+                    .flatMap { meta, bedgz_files ->
+                        def files = [ bedgz_files ].flatten()
+                        def phased = files.findAll { f -> f.name =~ /_hp[12]\\.bed\\.gz\$/ }
+                        if (phased) {
+                            phased.collect { f ->
+                                def haplotype = (f.name =~ /_(hp[12])\\.bed\\.gz\$/)[0][1]
+                                [ meta + [ id: "\${meta.id}_\${haplotype}" ], f, [], [] ]
+                            }
+                        } else {
+                            // -stub: modkit pileup's stub always emits one generic, unsplit file
+                            // (it never reads --phased in stub mode) -- reuse it for both
+                            // haplotypes, since MODKIT_DMR's own stub never reads its inputs.
+                            files.collectMany { f ->
+                                ['hp1', 'hp2'].collect { haplotype -> [ meta + [ id: "\${meta.id}_\${haplotype}" ], f, [], [] ] }
+                            }
+                        }
+                    }
+                input[1] = "compress"
+                input[2] = true
+                input[3] = []
+                """
+            }
+        }
+    }
+
+    test("[bedmethyl_a, tbi], [bedmethyl_b, tbi], regions_bed, fasta") {
+
+        when {
+            params {
+                module_args = '--base C'
+            }
+            process {
+                """
+                input[0] = HTSLIB_BGZIPTABIX.out.output
+                    .join(HTSLIB_BGZIPTABIX.out.index)
+                    .filter { meta, bedgz, tbi -> meta.id.endsWith('_hp1') }
+                    .map { meta, bedgz, tbi -> [ [ id: 'test' ], bedgz, tbi ] }
+                input[1] = HTSLIB_BGZIPTABIX.out.output
+                    .join(HTSLIB_BGZIPTABIX.out.index)
+                    .filter { meta, bedgz, tbi -> meta.id.endsWith('_hp2') }
+                    .map { meta, bedgz, tbi -> [ [ id: 'test' ], bedgz, tbi ] }
+                input[2] = Channel.of('chr22\t0\t40001')
+                    .collectFile(name: 'chr22.bed', newLine: true)
+                    .map { file -> [ [ id:'chr22' ], file ] }
+                input[3] = [
+                    [ id: 'test_ref' ],
+                    file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta', checkIfExists: true),
+                    file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta.fai', checkIfExists: true)
+                ]
+                """
+            }
+        }
+
+        then {
+            assertAll (
+                { assert process.success },
+                // log-file content includes a timestamp/runtime that varies between runs --
+                // assert existence only, and snapshot everything else.
+                { assert process.out.log.size() == 1 },
+                { assert snapshot(process.out.bed, process.out.findAll { key, val -> key.startsWith('versions') }
+                ).match() }
+            )
+        }
+
+    }
+
+    test("[bedmethyl_a, tbi], [bedmethyl_b, tbi], [], fasta") {
+
+        when {
+            params {
+                module_args = '--base C'
+            }
+            process {
+                """
+                input[0] = HTSLIB_BGZIPTABIX.out.output
+                    .join(HTSLIB_BGZIPTABIX.out.index)
+                    .filter { meta, bedgz, tbi -> meta.id.endsWith('_hp1') }
+                    .map { meta, bedgz, tbi -> [ [ id: 'test' ], bedgz, tbi ] }
+                input[1] = HTSLIB_BGZIPTABIX.out.output
+                    .join(HTSLIB_BGZIPTABIX.out.index)
+                    .filter { meta, bedgz, tbi -> meta.id.endsWith('_hp2') }
+                    .map { meta, bedgz, tbi -> [ [ id: 'test' ], bedgz, tbi ] }
+                input[2] = [ [ id:'regions' ], [] ]
+                input[3] = [
+                    [ id: 'test_ref' ],
+                    file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta', checkIfExists: true),
+                    file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta.fai', checkIfExists: true)
+                ]
+                """
+            }
+        }
+
+        then {
+            assertAll (
+                { assert process.success },
+                // log-file content includes a timestamp/runtime that varies between runs --
+                // assert existence only, and snapshot everything else.
+                { assert process.out.log.size() == 1 },
+                { assert snapshot(process.out.bed, process.out.findAll { key, val -> key.startsWith('versions') }
+                ).match() }
+            )
+        }
+
+    }
+
+    test("[bedmethyl_a, tbi], [bedmethyl_b, tbi], regions_bed, fasta - stub") {
+
+        options "-stub"
+
+        when {
+            params {
+                module_args = '--base C'
+            }
+            process {
+                """
+                input[0] = HTSLIB_BGZIPTABIX.out.output
+                    .join(HTSLIB_BGZIPTABIX.out.index)
+                    .filter { meta, bedgz, tbi -> meta.id.endsWith('_hp1') }
+                    .map { meta, bedgz, tbi -> [ [ id: 'test' ], bedgz, tbi ] }
+                input[1] = HTSLIB_BGZIPTABIX.out.output
+                    .join(HTSLIB_BGZIPTABIX.out.index)
+                    .filter { meta, bedgz, tbi -> meta.id.endsWith('_hp2') }
+                    .map { meta, bedgz, tbi -> [ [ id: 'test' ], bedgz, tbi ] }
+                input[2] = Channel.of('chr22\t0\t40001')
+                    .collectFile(name: 'chr22.bed', newLine: true)
+                    .map { file -> [ [ id:'chr22' ], file ] }
+                input[3] = [
+                    [ id: 'test_ref' ],
+                    file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta', checkIfExists: true),
+                    file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta.fai', checkIfExists: true)
+                ]
+                """
+            }
+        }
+
+        then {
+            assertAll (
+                { assert process.success },
+                // log-file content includes a timestamp/runtime that varies between runs --
+                // assert existence only, and snapshot everything else.
+                { assert process.out.log.size() == 1 },
+                { assert snapshot(process.out.bed, process.out.findAll { key, val -> key.startsWith('versions') }
+).match() }
+            )
+        }
+
+    }
+
+}