diff --git a/modules/nf-core/modkit/dmr/environment.yml b/modules/nf-core/modkit/dmr/environment.yml new file mode 100644 index 000000000000..62b9786305ae --- /dev/null +++ b/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 diff --git a/modules/nf-core/modkit/dmr/main.nf b/modules/nf-core/modkit/dmr/main.nf new file mode 100644 index 000000000000..0a5ce6799d25 --- /dev/null +++ b/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 + // .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 + """ +} diff --git a/modules/nf-core/modkit/dmr/meta.yml b/modules/nf-core/modkit/dmr/meta.yml new file mode 100644 index 000000000000..63c958daa17f --- /dev/null +++ b/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 + .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" diff --git a/modules/nf-core/modkit/dmr/tests/main.nf.test b/modules/nf-core/modkit/dmr/tests/main.nf.test new file mode 100644 index 000000000000..d7d4e6de8cb1 --- /dev/null +++ b/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() } + ) + } + + } + +} diff --git a/modules/nf-core/modkit/dmr/tests/main.nf.test.snap b/modules/nf-core/modkit/dmr/tests/main.nf.test.snap new file mode 100644 index 000000000000..a7039d7dac2d --- /dev/null +++ b/modules/nf-core/modkit/dmr/tests/main.nf.test.snap @@ -0,0 +1,80 @@ +{ + "[bedmethyl_a, tbi], [bedmethyl_b, tbi], regions_bed, fasta": { + "content": [ + [ + [ + { + "id": "test" + }, + "test.bed:md5,c57c0b35ce975dbafa19f3b4b5604dec" + ] + ], + { + "versions_modkit": [ + [ + "MODKIT_DMR", + "modkit", + "0.6.1" + ] + ] + } + ], + "meta": { + "nf-test": "0.9.0", + "nextflow": "25.10.2" + }, + "timestamp": "2026-09-29T11:00:29.834449326" + }, + "[bedmethyl_a, tbi], [bedmethyl_b, tbi], [], fasta": { + "content": [ + [ + [ + { + "id": "test" + }, + "test.bed:md5,e1ed04436fd4aee1c05fa4a612846417" + ] + ], + { + "versions_modkit": [ + [ + "MODKIT_DMR", + "modkit", + "0.6.1" + ] + ] + } + ], + "meta": { + "nf-test": "0.9.0", + "nextflow": "25.10.2" + }, + "timestamp": "2026-09-29T11:00:38.70209544" + }, + "[bedmethyl_a, tbi], [bedmethyl_b, tbi], regions_bed, fasta - stub": { + "content": [ + [ + [ + { + "id": "test" + }, + "test.bed:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + { + "versions_modkit": [ + [ + "MODKIT_DMR", + "modkit", + "0.6.1" + ] + ] + } + ], + "meta": { + "nf-test": "0.9.0", + "nextflow": "25.10.2" + }, + "timestamp": "2026-09-29T11:00:46.52021922" + } +} \ No newline at end of file diff --git a/modules/nf-core/modkit/dmr/tests/nextflow.config b/modules/nf-core/modkit/dmr/tests/nextflow.config new file mode 100644 index 000000000000..a6b69568628b --- /dev/null +++ b/modules/nf-core/modkit/dmr/tests/nextflow.config @@ -0,0 +1,14 @@ +process { + withName: 'MODKIT_PILEUP' { + // --modified-bases 5mC 5hmC (matching modkit/pileup's own phased test): if this BAM's + // basecalling calls both, requesting only one leaves the other counted as "other" in + // modkit's bedMethyl output, which MODKIT_DMR's internal consistency check rejects. + ext.args = '--phased --modified-bases 5mC 5hmC' + } + withName: 'HTSLIB_BGZIPTABIX' { + ext.args2 = '-p bed' + } + withName: 'MODKIT_DMR' { + ext.args = params.module_args + } +}