diff --git a/.github/workflows/linting.yml b/.github/workflows/linting.yml index 277f12fb..7ca0ebd9 100644 --- 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 index 82e0e8fd..250d345f 100644 --- 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 index 2c589afc..27f70a29 100644 --- 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 index d56f3ffe..41fada0d 100644 --- 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 index 9487ba34..bdf69b86 100644 --- 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. + ### `variants` diff --git a/docs/usage.md b/docs/usage.md index 3a64a4ac..044e9e02 100644 --- 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 | diff --git a/modules.json b/modules.json index 0db390bb..0a33c58a 100644 --- a/modules.json +++ b/modules.json @@ -88,6 +88,11 @@ "installed_by": ["modules"], "patch": "modules/nf-core/ensemblvep/vep/ensemblvep-vep.diff" }, + "htslib/bgziptabix": { + "branch": "master", + "git_sha": "cbe6025183336010a9b716e463f897b0c0f70f8a", + "installed_by": ["modules"] + }, "longphase/haplotag": { "branch": "master", "git_sha": "6d46786420b4d7bc88eba026eb389c0c5535d120", @@ -110,6 +115,12 @@ "git_sha": "135d883a43b1a2324b4112cd270cf22e20835533", "installed_by": ["modules"] }, + "modkit/dmr": { + "branch": "master", + "git_sha": "0f9aa70f68367b17b2a92043ff9614aa5792d73f", + "installed_by": ["modules"], + "patch": "modules/nf-core/modkit/dmr/modkit-dmr.diff" + }, "modkit/pileup": { "branch": "master", "git_sha": "3d81317a30d1016b533982d6b84df07713ae520a", diff --git a/modules/local/dmr/haplotype_regions/environment.yml b/modules/local/dmr/haplotype_regions/environment.yml new file mode 100644 index 00000000..8b02e7a2 --- /dev/null +++ b/modules/local/dmr/haplotype_regions/environment.yml @@ -0,0 +1,8 @@ +--- +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/environment-schema.json +channels: + - conda-forge + - bioconda +dependencies: + - bioconda::bedtools=2.31.1 + - bioconda::htslib=1.22.1 diff --git a/modules/local/dmr/haplotype_regions/main.nf b/modules/local/dmr/haplotype_regions/main.nf new file mode 100644 index 00000000..f536deaf --- /dev/null +++ b/modules/local/dmr/haplotype_regions/main.nf @@ -0,0 +1,58 @@ +process DMR_HAPLOTYPE_REGIONS { + tag "$meta.id" + label 'process_single' + + conda "${moduleDir}/environment.yml" + // The bedtools biocontainer lacks tabix; reuse nf-core pints/caller's image (bedtools + htslib) + container "${workflow.containerEngine in ['singularity', 'apptainer'] && !task.ext.singularity_pull_docker_container + ? 'https://community-cr-prod.seqera.io/docker/registry/v2/blobs/sha256/f1/f1a9e30012e1b41baf9acd1ff94e01161138d8aa17f4e97aa32f2dc4effafcd1/data' + : 'community.wave.seqera.io/library/pybedtools_bedtools_htslib_pip_pypints:39699b96998ec5f6'}" + + input: + tuple val(meta), path(hp1_bedmethyl), path(hp1_tbi), path(hp2_bedmethyl), path(hp2_tbi) + path cpg_islands_bed + tuple val(meta2), path(fai) + + output: + tuple val(meta), path("*.regions.bed"), emit: regions_bed + path "versions.yml" , emit: versions + + when: + task.ext.when == null || task.ext.when + + script: + def prefix = task.ext.prefix ?: "${meta.id}" + // Minimum Nvalid_cov (bedMethyl column 10) for a CpG to count as covered + def min_valid_coverage = task.ext.min_valid_coverage ?: 1 + """ + cut -f1,2 ${fai} > genome.txt + + tabix -R ${cpg_islands_bed} ${hp1_bedmethyl} \\ + | awk -v min_cov=${min_valid_coverage} 'BEGIN{OFS="\\t"} \$10>=min_cov {print \$1,\$2,\$3}' \\ + | sort -k1,1 -k2,2n -u > hp1.cpg.bed + tabix -R ${cpg_islands_bed} ${hp2_bedmethyl} \\ + | awk -v min_cov=${min_valid_coverage} 'BEGIN{OFS="\\t"} \$10>=min_cov {print \$1,\$2,\$3}' \\ + | sort -k1,1 -k2,2n -u > hp2.cpg.bed + + bedtools intersect -u -a ${cpg_islands_bed} -b hp1.cpg.bed | sort -k1,1 -k2,2n > islands.hp1.bed + bedtools intersect -u -a ${cpg_islands_bed} -b hp2.cpg.bed | sort -k1,1 -k2,2n > islands.hp2.bed + bedtools intersect -u -a islands.hp1.bed -b islands.hp2.bed | bedtools sort -g genome.txt -i - > ${prefix}.regions.bed + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + bedtools: \$(bedtools --version | sed 's/bedtools v//g') + tabix: \$(tabix --version 2>&1 | head -n1 | sed 's/tabix (htslib) //') + END_VERSIONS + """ + + stub: + def prefix = task.ext.prefix ?: "${meta.id}" + """ + touch ${prefix}.regions.bed + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + bedtools: 2.31.1 + END_VERSIONS + """ +} diff --git a/modules/local/dmr/haplotype_regions/meta.yml b/modules/local/dmr/haplotype_regions/meta.yml new file mode 100644 index 00000000..c3f80db6 --- /dev/null +++ b/modules/local/dmr/haplotype_regions/meta.yml @@ -0,0 +1,76 @@ +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/meta-schema.json +name: "dmr_haplotype_regions" +description: | + Restrict a CpG-islands BED to islands with real read coverage in both of a sample's two + haplotype bedMethyl files, for use as modkit dmr pair's --regions-bed. Without this, an island + covered in one haplotype but not the other would be compared against absence of data rather than + a real methylation difference. +keywords: + - methylation + - dmr + - haplotype + - cpg islands + - bedtools +tools: + - "bedtools": + description: "A powerful toolset for genome arithmetic" + homepage: "https://bedtools.readthedocs.io/" + documentation: "https://bedtools.readthedocs.io/" + licence: ["GPL-2.0-only"] + +input: + - - meta: + type: map + description: | + Groovy Map containing sample information, e.g. `[ id:'sample1' ]` + - hp1_bedmethyl: + type: file + description: Bgzip-compressed bedMethyl for haplotype 1 + pattern: "*.bed.gz" + - hp1_tbi: + type: file + description: Tabix index for hp1_bedmethyl + pattern: "*.bed.gz.tbi" + - hp2_bedmethyl: + type: file + description: Bgzip-compressed bedMethyl for haplotype 2 + pattern: "*.bed.gz" + - hp2_tbi: + type: file + description: Tabix index for hp2_bedmethyl + pattern: "*.bed.gz.tbi" + - - cpg_islands_bed: + type: file + description: Genome-wide CpG islands BED to restrict + pattern: "*.bed" + - - meta2: + type: map + description: | + Groovy Map containing reference information, e.g. `[ id:'GRCh38' ]` + - fai: + type: file + description: FASTA index for the reference genome (used only for its chrom/length columns) + pattern: "*.fai" + +output: + regions_bed: + - meta: + type: map + description: | + Groovy Map containing sample information, e.g. `[ id:'sample1' ]` + - "*.regions.bed": + type: file + description: | + cpg_islands_bed restricted to islands covered in both haplotypes, sorted by + genome/coordinate order, ready to pass to modkit dmr pair's --regions-bed + pattern: "*.regions.bed" + versions: + - "versions.yml": + type: file + description: File containing software versions + pattern: "versions.yml" + +authors: + - "@YannVRB" +maintainers: + - "@YannVRB" diff --git a/modules/local/dmr/nearest_gene/environment.yml b/modules/local/dmr/nearest_gene/environment.yml new file mode 100644 index 00000000..45c307b0 --- /dev/null +++ b/modules/local/dmr/nearest_gene/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: + - bioconda::bedtools=2.31.1 diff --git a/modules/local/dmr/nearest_gene/main.nf b/modules/local/dmr/nearest_gene/main.nf new file mode 100644 index 00000000..578cfdcb --- /dev/null +++ b/modules/local/dmr/nearest_gene/main.nf @@ -0,0 +1,46 @@ +process DMR_NEAREST_GENE { + tag "$meta.id" + label 'process_single' + + conda "${moduleDir}/environment.yml" + container "${workflow.containerEngine in ['singularity', 'apptainer'] && !task.ext.singularity_pull_docker_container + ? 'https://depot.galaxyproject.org/singularity/bedtools:2.31.1--hf5e1c6e_0' + : 'quay.io/biocontainers/bedtools:2.31.1--hf5e1c6e_0'}" + + input: + tuple val(meta), path(dmr_bed) + path gencode_gene_bed + + output: + tuple val(meta), path("*.nearestGene.bed"), emit: nearest_gene_bed + path "versions.yml" , emit: versions + + when: + task.ext.when == null || task.ext.when + + script: + // -D b: signed distance relative to the gene's strand (DMRs are unstranded) + def args = task.ext.args ?: '-D b -t first' + def prefix = task.ext.prefix ?: "${meta.id}" + """ + sort -k1,1 -k2,2n ${dmr_bed} > dmr.sorted.bed + + bedtools closest -a dmr.sorted.bed -b ${gencode_gene_bed} ${args} > ${prefix}.nearestGene.bed + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + bedtools: \$(bedtools --version | sed 's/bedtools v//g') + END_VERSIONS + """ + + stub: + def prefix = task.ext.prefix ?: "${meta.id}" + """ + touch ${prefix}.nearestGene.bed + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + bedtools: 2.31.1 + END_VERSIONS + """ +} diff --git a/modules/local/dmr/nearest_gene/meta.yml b/modules/local/dmr/nearest_gene/meta.yml new file mode 100644 index 00000000..39853526 --- /dev/null +++ b/modules/local/dmr/nearest_gene/meta.yml @@ -0,0 +1,51 @@ +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/meta-schema.json +name: "dmr_nearest_gene" +description: | + Sort a DMR bed (modkit dmr pair's output) by position and annotate each region with its nearest + gene from a GENCODE-derived gene BED, using bedtools closest. +keywords: + - methylation + - dmr + - gene annotation + - bedtools +tools: + - "bedtools": + description: "A powerful toolset for genome arithmetic" + homepage: "https://bedtools.readthedocs.io/" + documentation: "https://bedtools.readthedocs.io/" + licence: ["GPL-2.0-only"] + +input: + - - meta: + type: map + description: | + Groovy Map containing sample information, e.g. `[ id:'sample1' ]` + - dmr_bed: + type: file + description: modkit dmr pair's output BED + pattern: "*.bed" + - - gencode_gene_bed: + type: file + description: Genome-wide, position-sorted gene-level BED (chrom, start, end, gene_name, score, strand) + pattern: "*.bed" + +output: + nearest_gene_bed: + - meta: + type: map + description: | + Groovy Map containing sample information, e.g. `[ id:'sample1' ]` + - "*.nearestGene.bed": + type: file + description: dmr_bed with each region's nearest gene and signed distance appended + pattern: "*.nearestGene.bed" + versions: + - "versions.yml": + type: file + description: File containing software versions + pattern: "versions.yml" + +authors: + - "@YannVRB" +maintainers: + - "@YannVRB" diff --git a/modules/local/sigprofiler/verify/main.nf b/modules/local/sigprofiler/verify/main.nf index 4bf86194..959bcee1 100644 --- a/modules/local/sigprofiler/verify/main.nf +++ b/modules/local/sigprofiler/verify/main.nf @@ -38,7 +38,6 @@ process SIGPROFILER_VERIFY { manager = rgm.ReferenceGenomeManager("genome_volume") if not manager.is_genome_installed(genome): - manager.print_genome_checksum_verification_report(genome) tsb = manager.reference_dir.get_tsb_dir() / genome expected = rgm.CHECKSUMS[genome] missing = [chrom for chrom in expected if not (tsb / f"{chrom}.txt").is_file()] @@ -53,11 +52,15 @@ process SIGPROFILER_VERIFY { "payload (GRCh38 and CHM13-T2T installed for lrsomatic < 1.2.0 are a superseded revision) or the " "copy is corrupted." ) - sys.exit( + message = ( f"ERROR: the {genome} payload in --sigprofiler_genome_dir {cause} Reinstall it with " "--download_sigprofiler_genome (published to /cache/sigprofiler/volume) and pass that directory " "on later runs." ) + # Message before the checksum table, on the same stream: Nextflow's error report keeps only the first lines + print(message, flush=True) + manager.print_genome_checksum_verification_report(genome) + sys.exit(message) PY """ diff --git a/modules/nf-core/htslib/bgziptabix/environment.yml b/modules/nf-core/htslib/bgziptabix/environment.yml new file mode 100644 index 00000000..ec62c057 --- /dev/null +++ b/modules/nf-core/htslib/bgziptabix/environment.yml @@ -0,0 +1,8 @@ +--- +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/environment-schema.json +channels: + - conda-forge + - bioconda +dependencies: + - bioconda::htslib=1.24 + - conda-forge::xz=5.8.3 diff --git a/modules/nf-core/htslib/bgziptabix/main.nf b/modules/nf-core/htslib/bgziptabix/main.nf new file mode 100644 index 00000000..573c05fc --- /dev/null +++ b/modules/nf-core/htslib/bgziptabix/main.nf @@ -0,0 +1,88 @@ +process HTSLIB_BGZIPTABIX { + tag "${meta.id}" + label 'process_low' + + conda "${moduleDir}/environment.yml" + container "${workflow.containerEngine in ['singularity', 'apptainer'] && !task.ext.singularity_pull_docker_container + ? 'https://community-cr-prod.seqera.io/docker/registry/v2/blobs/sha256/86/863ca0dbbba30c8367fa4fbd3fa3a84393532fb7b300a5c5c2e70f0dfc475bbf/data' + : 'community.wave.seqera.io/library/htslib_xz:32f2772a564b3cd2'}" + + input: + tuple val(meta), path(infile), path(infile_tbi), path(regions) + val action + val make_index + val out_ext + + output: + tuple val(meta), path("${outfile}"), emit: output + tuple val(meta), path("${outfile}.{tbi,csi}"), emit: index, optional: true + // all htslib tools have the same version, we use bgzip + tuple val("${task.process}"), val('htslib'), eval("bgzip --version | sed '1! d; s/bgzip (htslib) //'"), topic: versions, emit: versions_htslib + tuple val("${task.process}"), val('xz'), eval("xz --version | sed '1! d; s/xz (XZ Utils) //'"), topic: versions, emit: versions_xz + + when: + task.ext.when == null || task.ext.when + + script: + def allowed_actions = ["compress", "decompress"] + if (action !in allowed_actions) { + error("htslib/bgziptabix: Invalid action: ${action}. Allowed actions are: ${allowed_actions.join(', ')}") + } + + if (action == "decompress" && make_index) { + log.warn("htslib/bgziptabix: Cannot create index when decompressing. Ignoring make_index option.") + } + + def args = task.ext.args ?: '' + def args2 = task.ext.args2 ?: '' + prefix = task.ext.prefix ?: "${meta.id}" + outfile = action == "compress" ? (out_ext ? "${prefix}.${out_ext}.gz" : "${prefix}.gz") : (out_ext ? "${prefix}.${out_ext}" : "${prefix}") + + def compress_cmd = action == "compress" ? "bgzip -c ${args} -@ ${task.cpus}" : "cat" + def bgzip_cmd = action == "compress" ? "[ '\$(basename ${infile})' != '\$(basename ${outfile})' ] && ln -s ${infile} ${outfile}" : "bgzip -c -d ${args} -@ ${task.cpus} ${infile} > ${outfile}" + + def regions_arg = regions ? "-R ${regions}" : "" + def tabix_cmd = (make_index && !infile_tbi) ? "tabix -@ ${task.cpus} ${regions_arg} ${args2} -f ${outfile}" : "" + def link_tabix_cmd = make_index && infile_tbi ? "ln -s ${infile_tbi} ${outfile}.${infile_tbi.extension}" : "" + def uncompressed_cmd = action == "compress" ? "${compress_cmd} ${infile} > ${outfile}" : (infile.getName() == outfile ? "" : "ln -s ${infile} ${outfile}") + """ + ${link_tabix_cmd} + + FILE_TYPE=\$(htsfile ${infile}) + + case "\$FILE_TYPE" in + *BGZF-compressed*) + ${bgzip_cmd} ;; + *gzip-compressed*) + [ "\$(basename ${infile})" == "\$(basename ${outfile})" ] && echo "Input and output names cannot be the same" && exit 1 + bgzip -d -c -@ ${task.cpus} ${infile} | ${compress_cmd} > ${outfile} ;; + *bzip2-compressed*) + bzcat ${infile} | ${compress_cmd} > ${outfile} ;; + *XZ-compressed*) + xzcat ${infile} | ${compress_cmd} > ${outfile} ;; + *) + ${uncompressed_cmd} ;; + esac + + ${tabix_cmd} + """ + + stub: + def args = task.ext.args ?: '' + def args2 = task.ext.args2 ?: '' + prefix = task.ext.prefix ?: "${meta.id}" + outfile = action == "compress" ? (out_ext ? "${prefix}.${out_ext}.gz" : "${prefix}.gz") : (out_ext ? "${prefix}.${out_ext}" : "${prefix}") + + def touch_cmd = action == "compress" ? "echo | bgzip -c" : "echo" + def index_fmt = args2.contains('-C') ? 'csi' : 'tbi' + def tabix_cmd = make_index ? "touch ${outfile}.${index_fmt}" : "" + def link_tabix_cmd = make_index && infile_tbi ? "ln -s ${infile_tbi} ${outfile}.${infile_tbi.extension}" : "" + """ + echo ${args} + + ${touch_cmd} > ${outfile} + + ${tabix_cmd} + ${link_tabix_cmd} + """ +} diff --git a/modules/nf-core/htslib/bgziptabix/meta.yml b/modules/nf-core/htslib/bgziptabix/meta.yml new file mode 100644 index 00000000..4cdefd0e --- /dev/null +++ b/modules/nf-core/htslib/bgziptabix/meta.yml @@ -0,0 +1,125 @@ +name: "htslib_bgziptabix" +description: "Multi-purpose module to compress, decompress and index files using bgzip + and tabix." +keywords: + - compress + - decompress + - index + - bgzip + - tabix + - gzip + - bzip + - xz +tools: + - "htslib": + description: "C library for high-throughput sequencing data formats." + homepage: "http://www.htslib.org/" + documentation: "http://www.htslib.org/doc/" + tool_dev_url: "https://github.com/samtools/htslib" + doi: "10.1093/gigascience/giab007" + licence: + - "MIT" + identifier: biotools:htslib +input: + - - meta: + type: map + description: | + Groovy Map containing sample information + e.g. [ id:'sample1' ] + - infile: + type: file + description: Input file to compress or decompress + pattern: "*" + ontologies: [] + - infile_tbi: + type: file + description: Optional tabix index for the input file. + pattern: "*.{tbi,csi}" + ontologies: + - edam: http://edamontology.org/format_3616 # tabix + - regions: + type: file + description: Optional file of regions to extract (BED or chr:start-end format). + Only used when creating an index for the output file. + pattern: "*.{bed,txt,tsv}" + ontologies: + - edam: http://edamontology.org/format_3475 # TSV + - edam: http://edamontology.org/format_3003 # BED + - action: + type: string + description: Action to perform, either `compress` or `decompress` + - make_index: + type: boolean + description: Whether to create a tabix index for the output file; only used + if `action` is `compress` + - out_ext: + type: string + description: Output file extension without `.gz` suffix (for example `vcf`) +output: + output: + - - meta: + type: map + description: | + Groovy Map containing sample information + e.g. [ id:'sample1' ] + - ${outfile}: + type: file + description: Compressed or decompressed output file + pattern: "*" + ontologies: [] + index: + - - meta: + type: map + description: | + Groovy Map containing sample information + e.g. [ id:'sample1' ] + - ${outfile}.{tbi,csi}: + type: file + description: Tabix index file for the compressed output file + pattern: "*.{tbi,csi}" + ontologies: + - edam: http://edamontology.org/format_3616 # tabix + versions_htslib: + - - ${task.process}: + type: string + description: The name of the process + - htslib: + type: string + description: The name of the tool + - bgzip --version | sed '1! d; s/bgzip (htslib) //': + type: eval + description: The expression to obtain the version of the tool + versions_xz: + - - ${task.process}: + type: string + description: The name of the process + - xz: + type: string + description: The name of the tool + - xz --version | sed '1! d; s/xz (XZ Utils) //': + type: eval + description: The expression to obtain the version of the tool +topics: + versions: + - - ${task.process}: + type: string + description: The name of the process + - htslib: + type: string + description: The name of the tool + - bgzip --version | sed '1! d; s/bgzip (htslib) //': + type: eval + description: The expression to obtain the version of the tool + - - ${task.process}: + type: string + description: The name of the process + - xz: + type: string + description: The name of the tool + - xz --version | sed '1! d; s/xz (XZ Utils) //': + type: eval + description: The expression to obtain the version of the tool +authors: + - "@itrujnara" +maintainers: + - "@itrujnara" diff --git a/modules/nf-core/htslib/bgziptabix/tests/main.nf.test b/modules/nf-core/htslib/bgziptabix/tests/main.nf.test new file mode 100644 index 00000000..a7346506 --- /dev/null +++ b/modules/nf-core/htslib/bgziptabix/tests/main.nf.test @@ -0,0 +1,435 @@ +nextflow_process { + + name "Test Process HTSLIB_BGZIPTABIX" + script "../main.nf" + process "HTSLIB_BGZIPTABIX" + + tag "modules" + tag "modules_nfcore" + tag "htslib" + tag "htslib/bgziptabix" + + test("sarscov2 - vcf - decompress") { + + when { + process { + """ + input[0] = [ + [ id:'test' ], + file(params.modules_testdata_base_path + 'genomics/sarscov2/illumina/vcf/test.vcf', checkIfExists: true), + [], + [] + ] + input[1] = 'decompress' // action + input[2] = false // make_index + input[3] = 'vcf' // out_ext + """ + } + } + + then { + assert process.success + assertAll( + { assert snapshot(sanitizeOutput(process.out)).match() }, + { assert process.out.output.get(0).get(1).endsWith('.vcf') }, + { assert process.out.index.size() == 0 } + ) + } + + } + + test("sarscov2 - vcf - compress - index") { + + when { + process { + """ + input[0] = [ + [ id:'test' ], + file(params.modules_testdata_base_path + 'genomics/sarscov2/illumina/vcf/test.vcf', checkIfExists: true), + [], + [] + ] + input[1] = 'compress' // action + input[2] = true // make_index + input[3] = 'vcf' // out_ext + """ + } + } + + then { + assert process.success + assertAll( + { assert snapshot(sanitizeOutput(process.out)).match() }, + { assert process.out.output.get(0).get(1).endsWith('.vcf.gz') }, + { assert process.out.index.get(0).get(1).endsWith('.vcf.gz.tbi') } + ) + } + + } + + test("sarscov2 - vcf + regions - compress - index") { + when { + process { + """ + input[0] = [ + [ id:'example' ], + file(params.modules_testdata_base_path + 'genomics/sarscov2/illumina/vcf/test.vcf.gz', checkIfExists: true), + file(params.modules_testdata_base_path + 'genomics/sarscov2/illumina/vcf/test.vcf.gz.tbi', checkIfExists: true), + file('https://raw.githubusercontent.com/luisas/test-datasets/refs/heads/add-bedgraph-subset-illumina/data/genomics/sarscov2/illumina/bed/test.bed', checkIfExists: true) + ] + input[1] = 'compress' // action + input[2] = true // make_index + input[3] = 'vcf' // out_ext + """ + } + } + + then { + assertAll ( + { assert process.success }, + { assert snapshot( + sanitizeOutput(process.out), + path(process.out.output[0][1]).vcf.getVariantsMD5(), + ).match() } + ) + } + } + + test("sarscov2 - bgzip - decompress") { + + when { + process { + """ + input[0] = [ + [ id:'test' ], + file(params.modules_testdata_base_path + 'genomics/sarscov2/illumina/vcf/test.vcf.gz', checkIfExists: true), + [], + [] + ] + input[1] = 'decompress' // action + input[2] = false // make_index + input[3] = 'vcf' // out_ext + """ + } + } + + then { + assert process.success + assertAll( + { assert snapshot(sanitizeOutput(process.out)).match() }, + { assert process.out.output.get(0).get(1).endsWith('.vcf') }, + { assert process.out.index.size() == 0 } + ) + } + } + + test("sarscov2 - bgzip - compress - no index") { + + when { + process { + """ + input[0] = [ + [ id:'test' ], + file(params.modules_testdata_base_path + 'genomics/sarscov2/illumina/vcf/test.vcf.gz', checkIfExists: true), + [], + [] + ] + input[1] = 'compress' // action + input[2] = false // make_index + input[3] = 'vcf' // out_ext + """ + } + } + + then { + assert process.success + assertAll( + { assert snapshot(sanitizeOutput(process.out)).match() }, + { assert process.out.output.get(0).get(1).endsWith('.vcf.gz') }, + { assert process.out.index.size() == 0 } + ) + } + + } + + test("sarscov2 - gzip - decompress") { + + when { + process { + """ + input[0] = [ + [ id:'test' ], + file(params.modules_testdata_base_path + 'genomics/sarscov2/illumina/fastq/test_1.fastq.gz', checkIfExists: true), + [], + [] + ] + input[1] = 'decompress' // action + input[2] = false // make_index + input[3] = 'fastq' // out_ext + """ + } + } + + then { + assert process.success + assertAll( + { assert snapshot(sanitizeOutput(process.out)).match() }, + { assert process.out.output.get(0).get(1).endsWith('.fastq') }, + { assert process.out.index.size() == 0 } + ) + } + } + + test("sarscov2 - gzip - (re)compress - no index") { + + when { + process { + """ + input[0] = [ + [ id:'test' ], + file(params.modules_testdata_base_path + 'genomics/sarscov2/illumina/fastq/test_1.fastq.gz', checkIfExists: true), + [], + [] + ] + input[1] = 'compress' // action + input[2] = false // make_index + input[3] = 'fastq' // out_ext + """ + } + } + + then { + assert process.success + assertAll( + { assert snapshot(sanitizeOutput(process.out)).match() }, + { assert process.out.output.get(0).get(1).endsWith('.fastq.gz') }, + { assert process.out.index.size() == 0 } + ) + } + } + + test("sarscov2 - gzip - name clash") { + + when { + process { + """ + input[0] = [ + [ id:'test_1' ], + file(params.modules_testdata_base_path + 'genomics/sarscov2/illumina/fastq/test_1.fastq.gz', checkIfExists: true), + [], + [] + ] + input[1] = 'compress' // action + input[2] = false // make_index + input[3] = 'fastq' // out_ext + """ + } + } + + then { + assert process.failed + assertAll( + { assert process.errorReport.contains("Input and output names cannot be the same") } + ) + } + } + + test("metagenome - bz2 - decompress") { + + when { + process { + """ + input[0] = [ + [ id:'test' ], + file(params.modules_testdata_base_path + 'genomics/prokaryotes/metagenome/rgi/card-data.tar.bz2', checkIfExists: true), + [], + [] + ] + input[1] = 'decompress' // action + input[2] = false // make_index + input[3] = 'tar' // out_ext + """ + } + } + + then { + assert process.success + assertAll( + { assert snapshot(sanitizeOutput(process.out)).match() }, + { assert process.out.output.get(0).get(1).endsWith('.tar') }, + { assert process.out.index.size() == 0 } + ) + } + } + + test("metagenome - bz2 - (re)compress - no index") { + + when { + process { + """ + input[0] = [ + [ id:'test' ], + file(params.modules_testdata_base_path + 'genomics/prokaryotes/metagenome/rgi/card-data.tar.bz2', checkIfExists: true), + [], + [] + ] + input[1] = 'compress' // action + input[2] = false // make_index + input[3] = 'tar' // out_ext + """ + } + } + + then { + assert process.success + assertAll( + { assert snapshot(sanitizeOutput(process.out)).match() }, + { assert process.out.output.get(0).get(1).endsWith('.tar.gz') }, + { assert process.out.index.size() == 0 } + ) + } + } + + test("metagenome - xz - decompress") { + + when { + process { + """ + input[0] = [ + [ id:'test' ], + file(params.modules_testdata_base_path + 'genomics/prokaryotes/metagenome/taxonomy/misc/taxa_sqlite.xz', checkIfExists: true), + [], + [] + ] + input[1] = 'decompress' // action + input[2] = false // make_index + input[3] = '' // out_ext + """ + } + } + + then { + assert process.success + assertAll( + { assert snapshot( + process.out, + process.out.findAll { key, val -> key.startsWith('versions') } + ).match() }, + { assert process.out.output.get(0).get(1).endsWith('test') }, + { assert process.out.index.size() == 0 } + ) + } + } + + test("metagenome - xz - (re)compress - no index") { + + when { + process { + """ + input[0] = [ + [ id:'test' ], + file(params.modules_testdata_base_path + 'genomics/prokaryotes/metagenome/taxonomy/misc/taxa_sqlite.xz', checkIfExists: true), + [], + [] + ] + input[1] = 'compress' // action + input[2] = false // make_index + input[3] = '' // out_ext + """ + } + } + + then { + assert process.success + assertAll( + { assert snapshot(sanitizeOutput(process.out)).match() }, + { assert process.out.output.get(0).get(1).endsWith('.gz') }, + { assert process.out.index.size() == 0 } + ) + } + } + + test("sarscov2 - vcf - compress - index - stub") { + + options "-stub" + + when { + process { + """ + input[0] = [ + [ id:'test' ], + file(params.modules_testdata_base_path + 'genomics/sarscov2/illumina/vcf/test.vcf', checkIfExists: true), + [], + [] + ] + input[1] = 'compress' // action + input[2] = true // make_index + input[3] = 'vcf' // out_ext + """ + } + } + + then { + assert process.success + assertAll( + { assert snapshot(sanitizeOutput(process.out)).match() } + ) + } + + } + + test("sarscov2 - vcf - decompress - stub") { + + options "-stub" + + when { + process { + """ + input[0] = [ + [ id:'test' ], + file(params.modules_testdata_base_path + 'genomics/sarscov2/illumina/vcf/test.vcf.gz', checkIfExists: true), + [], + [] + ] + input[1] = 'decompress' // action + input[2] = false // make_index + input[3] = 'vcf' // out_ext + """ + } + } + + then { + assert process.success + assertAll( + { assert snapshot(sanitizeOutput(process.out)).match() } + ) + } + + } + + test("illegal action") { + + when { + process { + """ + input[0] = [ + [ id:'test' ], + file(params.modules_testdata_base_path + 'genomics/sarscov2/illumina/vcf/test.vcf', checkIfExists: true), + [], + [] + ] + input[1] = 'invalid_action' // action + input[2] = true // make_index + input[3] = 'vcf' // out_ext + """ + } + } + + then { + assert process.failed + assert process.errorReport.contains("Invalid action: invalid_action. Allowed actions are: compress, decompress") + } + + } + +} diff --git a/modules/nf-core/htslib/bgziptabix/tests/main.nf.test.snap b/modules/nf-core/htslib/bgziptabix/tests/main.nf.test.snap new file mode 100644 index 00000000..4d8b3d25 --- /dev/null +++ b/modules/nf-core/htslib/bgziptabix/tests/main.nf.test.snap @@ -0,0 +1,527 @@ +{ + "sarscov2 - gzip - (re)compress - no index": { + "content": [ + { + "index": [ + + ], + "output": [ + [ + { + "id": "test" + }, + "test.fastq.gz:md5,4161df271f9bfcd25d5845a1e220dbec" + ] + ], + "versions_htslib": [ + [ + "HTSLIB_BGZIPTABIX", + "htslib", + "1.24" + ] + ], + "versions_xz": [ + [ + "HTSLIB_BGZIPTABIX", + "xz", + "5.8.3" + ] + ] + } + ], + "timestamp": "2026-07-10T08:07:32.981182", + "meta": { + "nf-test": "0.9.5", + "nextflow": "26.04.6" + } + }, + "metagenome - xz - (re)compress - no index": { + "content": [ + { + "index": [ + + ], + "output": [ + [ + { + "id": "test" + }, + "test.gz:md5,b8d852a2b1ee52ed64d83046dcdb9de2" + ] + ], + "versions_htslib": [ + [ + "HTSLIB_BGZIPTABIX", + "htslib", + "1.24" + ] + ], + "versions_xz": [ + [ + "HTSLIB_BGZIPTABIX", + "xz", + "5.8.3" + ] + ] + } + ], + "timestamp": "2026-07-10T08:08:57.748138", + "meta": { + "nf-test": "0.9.5", + "nextflow": "26.04.6" + } + }, + "metagenome - bz2 - decompress": { + "content": [ + { + "index": [ + + ], + "output": [ + [ + { + "id": "test" + }, + "test.tar:md5,39e9e71fd16cfd09ceca12cd46e6abce" + ] + ], + "versions_htslib": [ + [ + "HTSLIB_BGZIPTABIX", + "htslib", + "1.24" + ] + ], + "versions_xz": [ + [ + "HTSLIB_BGZIPTABIX", + "xz", + "5.8.3" + ] + ] + } + ], + "timestamp": "2026-07-10T08:07:44.215383", + "meta": { + "nf-test": "0.9.5", + "nextflow": "26.04.6" + } + }, + "sarscov2 - vcf - decompress - stub": { + "content": [ + { + "index": [ + + ], + "output": [ + [ + { + "id": "test" + }, + "test.vcf:md5,68b329da9893e34099c7d8ad5cb9c940" + ] + ], + "versions_htslib": [ + [ + "HTSLIB_BGZIPTABIX", + "htslib", + "1.24" + ] + ], + "versions_xz": [ + [ + "HTSLIB_BGZIPTABIX", + "xz", + "5.8.3" + ] + ] + } + ], + "timestamp": "2026-07-10T08:10:11.941036", + "meta": { + "nf-test": "0.9.5", + "nextflow": "26.04.6" + } + }, + "sarscov2 - gzip - decompress": { + "content": [ + { + "index": [ + + ], + "output": [ + [ + { + "id": "test" + }, + "test.fastq:md5,4161df271f9bfcd25d5845a1e220dbec" + ] + ], + "versions_htslib": [ + [ + "HTSLIB_BGZIPTABIX", + "htslib", + "1.24" + ] + ], + "versions_xz": [ + [ + "HTSLIB_BGZIPTABIX", + "xz", + "5.8.3" + ] + ] + } + ], + "timestamp": "2026-07-10T08:07:28.301585", + "meta": { + "nf-test": "0.9.5", + "nextflow": "26.04.6" + } + }, + "sarscov2 - vcf - compress - index - stub": { + "content": [ + { + "index": [ + [ + { + "id": "test" + }, + "test.vcf.gz.tbi:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "output": [ + [ + { + "id": "test" + }, + "test.vcf.gz:md5,68b329da9893e34099c7d8ad5cb9c940" + ] + ], + "versions_htslib": [ + [ + "HTSLIB_BGZIPTABIX", + "htslib", + "1.24" + ] + ], + "versions_xz": [ + [ + "HTSLIB_BGZIPTABIX", + "xz", + "5.8.3" + ] + ] + } + ], + "timestamp": "2026-07-10T08:09:24.961486", + "meta": { + "nf-test": "0.9.5", + "nextflow": "26.04.6" + } + }, + "sarscov2 - vcf - decompress": { + "content": [ + { + "index": [ + + ], + "output": [ + [ + { + "id": "test" + }, + "test.vcf:md5,8e722884ffb75155212a3fc053918766" + ] + ], + "versions_htslib": [ + [ + "HTSLIB_BGZIPTABIX", + "htslib", + "1.24" + ] + ], + "versions_xz": [ + [ + "HTSLIB_BGZIPTABIX", + "xz", + "5.8.3" + ] + ] + } + ], + "timestamp": "2026-07-10T08:07:05.41219", + "meta": { + "nf-test": "0.9.5", + "nextflow": "26.04.6" + } + }, + "metagenome - bz2 - (re)compress - no index": { + "content": [ + { + "index": [ + + ], + "output": [ + [ + { + "id": "test" + }, + "test.tar.gz:md5,39e9e71fd16cfd09ceca12cd46e6abce" + ] + ], + "versions_htslib": [ + [ + "HTSLIB_BGZIPTABIX", + "htslib", + "1.24" + ] + ], + "versions_xz": [ + [ + "HTSLIB_BGZIPTABIX", + "xz", + "5.8.3" + ] + ] + } + ], + "timestamp": "2026-07-10T08:07:51.772557", + "meta": { + "nf-test": "0.9.5", + "nextflow": "26.04.6" + } + }, + "sarscov2 - vcf - compress - index": { + "content": [ + { + "index": [ + [ + { + "id": "test" + }, + "test.vcf.gz.tbi:md5,7f005943c935f2b55ba3f9d4802aa09f" + ] + ], + "output": [ + [ + { + "id": "test" + }, + "test.vcf.gz:md5,8e722884ffb75155212a3fc053918766" + ] + ], + "versions_htslib": [ + [ + "HTSLIB_BGZIPTABIX", + "htslib", + "1.24" + ] + ], + "versions_xz": [ + [ + "HTSLIB_BGZIPTABIX", + "xz", + "5.8.3" + ] + ] + } + ], + "timestamp": "2026-07-10T08:07:09.68169", + "meta": { + "nf-test": "0.9.5", + "nextflow": "26.04.6" + } + }, + "metagenome - xz - decompress": { + "content": [ + { + "0": [ + [ + { + "id": "test" + }, + "test:md5,b8d852a2b1ee52ed64d83046dcdb9de2" + ] + ], + "1": [ + + ], + "2": [ + [ + "HTSLIB_BGZIPTABIX", + "htslib", + "1.24" + ] + ], + "3": [ + [ + "HTSLIB_BGZIPTABIX", + "xz", + "5.8.3" + ] + ], + "index": [ + + ], + "output": [ + [ + { + "id": "test" + }, + "test:md5,b8d852a2b1ee52ed64d83046dcdb9de2" + ] + ], + "versions_htslib": [ + [ + "HTSLIB_BGZIPTABIX", + "htslib", + "1.24" + ] + ], + "versions_xz": [ + [ + "HTSLIB_BGZIPTABIX", + "xz", + "5.8.3" + ] + ] + }, + { + "versions_htslib": [ + [ + "HTSLIB_BGZIPTABIX", + "htslib", + "1.24" + ] + ], + "versions_xz": [ + [ + "HTSLIB_BGZIPTABIX", + "xz", + "5.8.3" + ] + ] + } + ], + "timestamp": "2026-07-10T08:08:19.920765", + "meta": { + "nf-test": "0.9.5", + "nextflow": "26.04.6" + } + }, + "sarscov2 - bgzip - compress - no index": { + "content": [ + { + "index": [ + + ], + "output": [ + [ + { + "id": "test" + }, + "test.vcf.gz:md5,8e722884ffb75155212a3fc053918766" + ] + ], + "versions_htslib": [ + [ + "HTSLIB_BGZIPTABIX", + "htslib", + "1.24" + ] + ], + "versions_xz": [ + [ + "HTSLIB_BGZIPTABIX", + "xz", + "5.8.3" + ] + ] + } + ], + "timestamp": "2026-07-10T08:07:23.405306", + "meta": { + "nf-test": "0.9.5", + "nextflow": "26.04.6" + } + }, + "sarscov2 - vcf + regions - compress - index": { + "content": [ + { + "index": [ + [ + { + "id": "example" + }, + "example.vcf.gz.tbi:md5,d22e5b84e4fcd18792179f72e6da702e" + ] + ], + "output": [ + [ + { + "id": "example" + }, + "example.vcf.gz:md5,8e722884ffb75155212a3fc053918766" + ] + ], + "versions_htslib": [ + [ + "HTSLIB_BGZIPTABIX", + "htslib", + "1.24" + ] + ], + "versions_xz": [ + [ + "HTSLIB_BGZIPTABIX", + "xz", + "5.8.3" + ] + ] + }, + "bc7bf3ee9e8430e064c539eb81e59bf9" + ], + "timestamp": "2026-07-10T08:07:15.187648", + "meta": { + "nf-test": "0.9.5", + "nextflow": "26.04.6" + } + }, + "sarscov2 - bgzip - decompress": { + "content": [ + { + "index": [ + + ], + "output": [ + [ + { + "id": "test" + }, + "test.vcf:md5,8e722884ffb75155212a3fc053918766" + ] + ], + "versions_htslib": [ + [ + "HTSLIB_BGZIPTABIX", + "htslib", + "1.24" + ] + ], + "versions_xz": [ + [ + "HTSLIB_BGZIPTABIX", + "xz", + "5.8.3" + ] + ] + } + ], + "timestamp": "2026-07-10T08:07:19.313945", + "meta": { + "nf-test": "0.9.5", + "nextflow": "26.04.6" + } + } +} \ No newline at end of file diff --git a/modules/nf-core/modkit/dmr/environment.yml b/modules/nf-core/modkit/dmr/environment.yml new file mode 100644 index 00000000..62b97863 --- /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 00000000..41ccec1c --- /dev/null +++ b/modules/nf-core/modkit/dmr/main.nf @@ -0,0 +1,60 @@ +process MODKIT_DMR { + tag "${meta.id}" + label 'process_medium' + + // Conda is not supported: matches MODKIT_PILEUP's own patch (the guard in `script:` stops + // conda/mamba runs) -- kept on the same patched build as pileup for a single consistent + // modkit binary across the whole DMR chain, rather than mixing stock 0.6.1 (dmr pair) with + // the patched 0.6.4 (pileup). Revert to the biocontainer once the patch is upstreamed/released. + container "${(workflow.containerEngine == 'singularity' || workflow.containerEngine == 'apptainer') && !task.ext.singularity_pull_docker_container + ? 'oras://ghcr.io/ljwharbers/modkit-sif:0.6.4-pacbiofix-6e0afa2' + : 'ghcr.io/ljwharbers/modkit:0.6.4-pacbiofix-6e0afa2'}" + + 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: + // Exit if running this module with -profile conda / -profile mamba + if (workflow.profile.tokenize(',').intersect(['conda', 'mamba']).size() >= 1) { + error "MODKIT_DMR does not support Conda: the patched container is Docker/Singularity/Apptainer-only. Use one of those, or --skip_dmr." + } + 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 00000000..63c958da --- /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/modkit-dmr.diff b/modules/nf-core/modkit/dmr/modkit-dmr.diff new file mode 100644 index 00000000..8a99ae4f --- /dev/null +++ b/modules/nf-core/modkit/dmr/modkit-dmr.diff @@ -0,0 +1,148 @@ +Changes in component 'nf-core/modkit/dmr' +'modules/nf-core/modkit/dmr/environment.yml' is unchanged +'modules/nf-core/modkit/dmr/meta.yml' is unchanged +Changes in 'modkit/dmr/main.nf': +--- modules/nf-core/modkit/dmr/main.nf ++++ modules/nf-core/modkit/dmr/main.nf +@@ -2,10 +2,13 @@ + 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'}" ++ // Conda is not supported: matches MODKIT_PILEUP's own patch (the guard in `script:` stops ++ // conda/mamba runs) -- kept on the same patched build as pileup for a single consistent ++ // modkit binary across the whole DMR chain, rather than mixing stock 0.6.1 (dmr pair) with ++ // the patched 0.6.4 (pileup). Revert to the biocontainer once the patch is upstreamed/released. ++ container "${(workflow.containerEngine == 'singularity' || workflow.containerEngine == 'apptainer') && !task.ext.singularity_pull_docker_container ++ ? 'oras://ghcr.io/ljwharbers/modkit-sif:0.6.4-pacbiofix-6e0afa2' ++ : 'ghcr.io/ljwharbers/modkit:0.6.4-pacbiofix-6e0afa2'}" + + input: + tuple val(meta), path(bedmethyl_a), path(bedmethyl_a_tbi) +@@ -22,6 +25,10 @@ + task.ext.when == null || task.ext.when + + script: ++ // Exit if running this module with -profile conda / -profile mamba ++ if (workflow.profile.tokenize(',').intersect(['conda', 'mamba']).size() >= 1) { ++ error "MODKIT_DMR does not support Conda: the patched container is Docker/Singularity/Apptainer-only. Use one of those, or --skip_dmr." ++ } + def args = task.ext.args ?: '' + def prefix = task.ext.prefix ?: "${meta.id}" + def regions = regions_bed ? "-r ${regions_bed}" : '' + +'modules/nf-core/modkit/dmr/tests/nextflow.config' is unchanged +Changes in 'modkit/dmr/tests/main.nf.test': +--- modules/nf-core/modkit/dmr/tests/main.nf.test ++++ modules/nf-core/modkit/dmr/tests/main.nf.test +@@ -21,12 +21,16 @@ + 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) + ] ++ // lrsomatic's patched modkit/pileup takes fasta and fai as separate inputs + 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', checkIfExists: true) ++ ] ++ input[2] = [ ++ [ id: 'test_ref' ], + file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta.fai', checkIfExists: true) + ] +- input[2] = [[],[]] ++ input[3] = [[:],[]] + """ + } + } + +Changes in 'modkit/dmr/tests/main.nf.test.snap': +--- modules/nf-core/modkit/dmr/tests/main.nf.test.snap ++++ modules/nf-core/modkit/dmr/tests/main.nf.test.snap +@@ -6,7 +6,7 @@ + { + "id": "test" + }, +- "test.bed:md5,c57c0b35ce975dbafa19f3b4b5604dec" ++ "test.bed:md5,88e4b55db4eaf1e77382a81bf89f361c" + ] + ], + { +@@ -14,16 +14,16 @@ + [ + "MODKIT_DMR", + "modkit", +- "0.6.1" ++ "0.6.4" + ] + ] + } + ], ++ "timestamp": "2026-10-26T13:39:52.386460283", + "meta": { +- "nf-test": "0.9.0", +- "nextflow": "25.10.2" +- }, +- "timestamp": "2026-09-29T11:00:29.834449326" ++ "nf-test": "0.9.4", ++ "nextflow": "26.04.3" ++ } + }, + "[bedmethyl_a, tbi], [bedmethyl_b, tbi], [], fasta": { + "content": [ +@@ -32,7 +32,7 @@ + { + "id": "test" + }, +- "test.bed:md5,e1ed04436fd4aee1c05fa4a612846417" ++ "test.bed:md5,64d0d1131027300a5d9766bf29c35c67" + ] + ], + { +@@ -40,16 +40,16 @@ + [ + "MODKIT_DMR", + "modkit", +- "0.6.1" ++ "0.6.4" + ] + ] + } + ], ++ "timestamp": "2026-10-26T13:40:10.750011704", + "meta": { +- "nf-test": "0.9.0", +- "nextflow": "25.10.2" +- }, +- "timestamp": "2026-09-29T11:00:38.70209544" ++ "nf-test": "0.9.4", ++ "nextflow": "26.04.3" ++ } + }, + "[bedmethyl_a, tbi], [bedmethyl_b, tbi], regions_bed, fasta - stub": { + "content": [ +@@ -66,15 +66,15 @@ + [ + "MODKIT_DMR", + "modkit", +- "0.6.1" ++ "0.6.4" + ] + ] + } + ], ++ "timestamp": "2026-10-26T13:40:27.769775341", + "meta": { +- "nf-test": "0.9.0", +- "nextflow": "25.10.2" +- }, +- "timestamp": "2026-09-29T11:00:46.52021922" ++ "nf-test": "0.9.4", ++ "nextflow": "26.04.3" ++ } + } +-}+} + +************************************************************ 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 00000000..96195470 --- /dev/null +++ b/modules/nf-core/modkit/dmr/tests/main.nf.test @@ -0,0 +1,198 @@ +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) + ] + // lrsomatic's patched modkit/pileup takes fasta and fai as separate inputs + input[1] = [ + [ id: 'test_ref' ], + file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta', checkIfExists: true) + ] + input[2] = [ + [ id: 'test_ref' ], + file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta.fai', checkIfExists: true) + ] + input[3] = [[:],[]] + """ + } + } + + 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 00000000..b25206c5 --- /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,88e4b55db4eaf1e77382a81bf89f361c" + ] + ], + { + "versions_modkit": [ + [ + "MODKIT_DMR", + "modkit", + "0.6.4" + ] + ] + } + ], + "timestamp": "2026-10-26T13:39:52.386460283", + "meta": { + "nf-test": "0.9.4", + "nextflow": "26.04.3" + } + }, + "[bedmethyl_a, tbi], [bedmethyl_b, tbi], [], fasta": { + "content": [ + [ + [ + { + "id": "test" + }, + "test.bed:md5,64d0d1131027300a5d9766bf29c35c67" + ] + ], + { + "versions_modkit": [ + [ + "MODKIT_DMR", + "modkit", + "0.6.4" + ] + ] + } + ], + "timestamp": "2026-10-26T13:40:10.750011704", + "meta": { + "nf-test": "0.9.4", + "nextflow": "26.04.3" + } + }, + "[bedmethyl_a, tbi], [bedmethyl_b, tbi], regions_bed, fasta - stub": { + "content": [ + [ + [ + { + "id": "test" + }, + "test.bed:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + { + "versions_modkit": [ + [ + "MODKIT_DMR", + "modkit", + "0.6.4" + ] + ] + } + ], + "timestamp": "2026-10-26T13:40:27.769775341", + "meta": { + "nf-test": "0.9.4", + "nextflow": "26.04.3" + } + } +} 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 00000000..a6b69568 --- /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 + } +} diff --git a/nextflow.config b/nextflow.config index 3a2a0f3d..836884c5 100644 --- a/nextflow.config +++ b/nextflow.config @@ -30,6 +30,11 @@ params { modkit_args = '--cpg --modified-bases 5mC' modkit_phased = false + // DMR options (only with modkit_phased); the BEDs default per --genome via igenomes.config + skip_dmr = false + dmr_cpg_islands_bed = null + dmr_gencode_gene_bed = null + // PON Options clairsto_pon_vcfs = null clairsto_pon_flags = null diff --git a/nextflow_schema.json b/nextflow_schema.json index a7234ecf..46c7fe25 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -135,7 +135,7 @@ "modkit_args": { "type": "string", "description": "Command line arguments passed to `modkit pileup`. The value replaces the default rather than extending it; the default restricts the bedMethyl output to 5mC calls at CpG sites.", - "help_text": "Repeat `--cpg --modified-bases 5mC` when adding flags (e.g. `--cpg --modified-bases 5mC --combine-strands`). An empty value gives the unfiltered pileup of every modification code at every position: use `--modkit_args=''` (the `=` form) or an empty `modkit_args` entry in a params file, since `--modkit_args ''` reaches the pipeline as `true`. `--modified-bases` only filters the output: on PacBio data, positions where the 5mC and 5hmC probabilities sum above 1 are still dropped and non-conflicting 5hmC calls are counted in `N_other`.", + "help_text": "Repeat `--cpg --modified-bases 5mC` when adding flags (e.g. `--cpg --modified-bases 5mC --combine-strands`). An empty value gives the unfiltered pileup of every modification code at every position: use `--modkit_args=''` (the `=` form) or an empty `modkit_args` entry in a params file, since `--modkit_args ''` reaches the pipeline as `true`. `--modified-bases` only filters the output: on PacBio data, positions where the 5mC and 5hmC probabilities sum above 1 are still dropped and non-conflicting 5hmC calls are counted in `N_other`. When `--modkit_phased` is set and DMR will run, `5hmC` is appended automatically if `--modified-bases` names `5mC` without it, and the published `modkit_pileup` files then include the 5hmC rows: `modkit dmr pair` rejects bedMethyl rows where an uncounted modification code (e.g. real dual 5mC/5hmC SUP calls) pushes `N_other` above zero. `--skip_dmr` keeps the value as given.", "fa_icon": "fas fa-terminal", "default": "--cpg --modified-bases 5mC" }, @@ -179,6 +179,20 @@ "description": "Verdict CNA resource directory for ClairS-TO, holding loci_files/, allele_files/ and one GC_*.txt. Read only with --skip_ascat.", "help_text": "Read only on `--skip_ascat` runs; otherwise Verdict's germline tagging comes from the pipeline's ASCAT run and this is ignored. The image ships the GRCh38 set and `--genome CHM13 --skip_ascat` builds a CHM13 one, so this is only needed for another assembly or as an override. The directory must be self-contained, with no links reaching outside it. Add an `RT_.txt` to correct LogR for replication timing as well; without it, correction is GC-only.", "fa_icon": "fas fa-folder-open" + }, + "dmr_cpg_islands_bed": { + "type": "string", + "format": "file-path", + "description": "CpG-islands BED to restrict haplotype DMR calling to. Has a built-in default for --genome GRCh38/CHM13; required as an explicit override for any other genome when --modkit_phased is set and --skip_dmr is not.", + "help_text": "Only used by the DMR subworkflow (see --skip_dmr). Without a value -- built-in or supplied -- DMR is skipped with a warning rather than the run failing.", + "fa_icon": "fas fa-map-marker-alt" + }, + "dmr_gencode_gene_bed": { + "type": "string", + "format": "file-path", + "description": "Gene-model BED (used for nearest-gene annotation of DMR calls) to use for haplotype DMR calling. Has a built-in default for --genome GRCh38/CHM13; required as an explicit override for any other genome when --modkit_phased is set and --skip_dmr is not.", + "help_text": "Only used by the DMR subworkflow (see --skip_dmr). Without a value -- built-in or supplied -- DMR is skipped with a warning rather than the run failing.", + "fa_icon": "fas fa-dna" } } }, @@ -673,6 +687,10 @@ "type": "boolean", "description": "Skip SAVANA (SV + copy-number calling)" }, + "skip_dmr": { + "type": "boolean", + "description": "Skip differential methylation region (DMR) calling between haplotypes. Only takes effect when --modkit_phased is also set; has no effect otherwise." + }, "skip_padfoot": { "type": "boolean", "help_text": "Skips Padfoot annotation of somatic SVs and CNAs (Severus + Wakhan and SAVANA pairs). Padfoot runs only inside its container; skip it under `-profile conda`.", diff --git a/subworkflows/local/dmr.nf b/subworkflows/local/dmr.nf new file mode 100644 index 00000000..5d869c26 --- /dev/null +++ b/subworkflows/local/dmr.nf @@ -0,0 +1,137 @@ +// IMPORT MODULES +include { HTSLIB_BGZIPTABIX } from '../../modules/nf-core/htslib/bgziptabix/main' +include { MODKIT_DMR } from '../../modules/nf-core/modkit/dmr/main' +include { DMR_HAPLOTYPE_REGIONS } from '../../modules/local/dmr/haplotype_regions/main' +include { DMR_NEAREST_GENE } from '../../modules/local/dmr/nearest_gene/main' + +workflow DMR { + + take: + modkit_bedgz // [meta, bed.gz(es)] -- MODKIT_PILEUP.out.bedgz; only meaningful when modkit_phased is true + fasta // [[:], fasta] + fai // [[:], fai] + cpg_islands_bed // path + gencode_gene_bed // path + + main: + ch_versions = channel.empty() + + // Keep the hp1/hp2 pair; unphased output has neither and is dropped + modkit_bedgz + .map { meta, files -> + def flist = files instanceof List ? files : [files] + def hp1 = flist.find { it.name.endsWith('_hp1.bed.gz') } + def hp2 = flist.find { it.name.endsWith('_hp2.bed.gz') } + return [meta, hp1, hp2] + } + .filter { _meta, hp1, hp2 -> hp1 && hp2 } + .set { haplotype_bedmethyl } + // haplotype_bedmethyl: [meta, hp1_bedgz, hp2_bedgz] + + // + // MODULE: HTSLIB_BGZIPTABIX (label: process_low) + // Tabix-index each haplotype (already bgzipped by modkit --bgzf) + // + haplotype_bedmethyl + .flatMap { meta, hp1, hp2 -> + return [ + [meta + [haplotype: 'hp1'], hp1, [], []], + [meta + [haplotype: 'hp2'], hp2, [], []] + ] + } + .set { bgziptabix_input } + // bgziptabix_input: [meta+haplotype, bedmethyl, [], []] + + HTSLIB_BGZIPTABIX ( + bgziptabix_input, + 'compress', + true, + 'bed' + ) + + HTSLIB_BGZIPTABIX.out.output + .join(HTSLIB_BGZIPTABIX.out.index) + .map { meta, bedgz, tbi -> + def haplotype = meta.haplotype + def sample_meta = meta.findAll { it.key != 'haplotype' } + return [sample_meta, haplotype, bedgz, tbi] + } + .branch { _meta, haplotype, _bedgz, _tbi -> + hp1: haplotype == 'hp1' + hp2: haplotype == 'hp2' + } + .set { indexed_branched } + + indexed_branched.hp1 + .map { meta, _haplotype, bedgz, tbi -> [meta, bedgz, tbi] } + .set { hp1_indexed } + indexed_branched.hp2 + .map { meta, _haplotype, bedgz, tbi -> [meta, bedgz, tbi] } + .set { hp2_indexed } + + hp1_indexed + .join(hp2_indexed) + .set { dmr_haplotype_input } + // dmr_haplotype_input: [meta, hp1_bedgz, hp1_tbi, hp2_bedgz, hp2_tbi] + + // + // MODULE: DMR_HAPLOTYPE_REGIONS (label: process_single) + // Restrict cpg_islands_bed to islands covered in both haplotypes. + // + DMR_HAPLOTYPE_REGIONS ( + dmr_haplotype_input, + cpg_islands_bed, + fai.first() + ) + ch_versions = ch_versions.mix(DMR_HAPLOTYPE_REGIONS.out.versions) + + // modkit dmr pair fails on an empty regions file: skip those samples + DMR_HAPLOTYPE_REGIONS.out.regions_bed + .filter { meta, regions_bed -> + if (regions_bed.size() == 0) { + log.warn("DMR: no CpG island is covered in both haplotypes for ${meta.id} (${meta.type}); skipping modkit dmr pair") + return false + } + return true + } + .set { nonempty_regions_bed } + + // + // MODULE: MODKIT_DMR (label: process_medium) + // Compare methylation between the two haplotypes over the restricted regions. + // + dmr_haplotype_input + .join(nonempty_regions_bed) + .multiMap { meta, hp1_bedgz, hp1_tbi, hp2_bedgz, hp2_tbi, regions_bed -> + hp1: [meta, hp1_bedgz, hp1_tbi] + hp2: [meta, hp2_bedgz, hp2_tbi] + regions: [meta, regions_bed] + } + .set { modkit_dmr_input } + + fasta.first() + .combine(fai.first()) + .map { meta, fasta_file, _meta2, fai_file -> [meta, fasta_file, fai_file] } + .set { fasta_with_fai } + + MODKIT_DMR ( + modkit_dmr_input.hp1, + modkit_dmr_input.hp2, + modkit_dmr_input.regions, + fasta_with_fai + ) + + // + // MODULE: DMR_NEAREST_GENE (label: process_single) + // Annotate each DMR with its nearest gene. + // + DMR_NEAREST_GENE ( + MODKIT_DMR.out.bed, + gencode_gene_bed + ) + ch_versions = ch_versions.mix(DMR_NEAREST_GENE.out.versions) + + emit: + nearest_gene_bed = DMR_NEAREST_GENE.out.nearest_gene_bed // [meta, bed] -- DMRs annotated with nearest gene + versions = ch_versions +} diff --git a/subworkflows/local/tests/dmr.nf.test b/subworkflows/local/tests/dmr.nf.test new file mode 100644 index 00000000..d8a88a05 --- /dev/null +++ b/subworkflows/local/tests/dmr.nf.test @@ -0,0 +1,128 @@ +nextflow_workflow { + + name "Test Workflow DMR" + script "../dmr.nf" + workflow "DMR" + config "./nextflow.config" + + tag "subworkflows" + tag "subworkflows_local" + tag "dmr" + // "small" is what .github/workflows/nf-test.yml selects on for pull_request + tag "small" + + // Full DMR chain on the public phased nanopore BAM, with a synthetic chr22 island and gene + setup { + run("MODKIT_PILEUP") { + script "../../../modules/nf-core/modkit/pileup/main.nf" + process { + """ + // lrsomatic's patched pileup takes fasta and fai separately, plus a bed tuple + input[0] = [ + [ id: 'test', type: 'tumor' ], + 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) + ] + input[2] = [ + [ id: 'test_ref' ], + file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta.fai', checkIfExists: true) + ] + input[3] = [[:],[]] + """ + } + } + } + + test("real phased public data -> haplotype DMR, annotated with nearest gene") { + + when { + workflow { + """ + def cpg_islands_bed = Channel.of('chr22\\t0\\t40001') + .collectFile(name: 'cpg_islands.bed', newLine: true) + + def gencode_gene_bed = Channel.of('chr22\\t0\\t40001\\tTEST_GENE\\t.\\t+') + .collectFile(name: 'gencode_gene.bed', newLine: true) + + // Value channels: DMR calls .first() on fasta/fai + input[0] = MODKIT_PILEUP.out.bedgz + input[1] = channel.value([ + [ id: 'test_ref' ], + file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta', checkIfExists: true) + ]) + input[2] = channel.value([ + [ id: 'test_ref' ], + file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta.fai', checkIfExists: true) + ]) + input[3] = cpg_islands_bed + input[4] = gencode_gene_bed + """ + } + } + + then { + def processes = workflow.trace.tasks().collect { task -> task.name } + + assertAll( + { assert workflow.success }, + { assert processes.any { name -> name.contains('HTSLIB_BGZIPTABIX') } }, + { assert processes.any { name -> name.contains('DMR_HAPLOTYPE_REGIONS') } }, + { assert processes.any { name -> name.contains('MODKIT_DMR') } }, + { assert processes.any { name -> name.contains('DMR_NEAREST_GENE') } }, + { assert workflow.out.nearest_gene_bed.size() == 1 }, + // Real content, not just a file existing: at least one DMR call, annotated + // with the synthetic gene above. + { + def bed = workflow.out.nearest_gene_bed[0][1] + def lines = file(bed).readLines() + assert lines.size() > 0 + assert lines.every { line -> line.contains('TEST_GENE') } + } + ) + } + } + + test("no CpG island covered in both haplotypes -> sample skipped, run succeeds") { + + when { + workflow { + """ + // An island on a contig the BAM has no reads on: DMR_HAPLOTYPE_REGIONS writes an + // empty regions BED, which modkit dmr pair would reject. + def cpg_islands_bed = Channel.of('chr1\\t0\\t1000') + .collectFile(name: 'cpg_islands.bed', newLine: true) + + def gencode_gene_bed = Channel.of('chr22\\t0\\t40001\\tTEST_GENE\\t.\\t+') + .collectFile(name: 'gencode_gene.bed', newLine: true) + + input[0] = MODKIT_PILEUP.out.bedgz + input[1] = channel.value([ + [ id: 'test_ref' ], + file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta', checkIfExists: true) + ]) + input[2] = channel.value([ + [ id: 'test_ref' ], + file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta.fai', checkIfExists: true) + ]) + input[3] = cpg_islands_bed + input[4] = gencode_gene_bed + """ + } + } + + then { + def processes = workflow.trace.tasks().collect { task -> task.name } + + assertAll( + { assert workflow.success }, + { assert processes.any { name -> name.contains('DMR_HAPLOTYPE_REGIONS') } }, + { assert !processes.any { name -> name.contains('MODKIT_DMR') } }, + { assert workflow.out.nearest_gene_bed.size() == 0 } + ) + } + } +} diff --git a/subworkflows/local/tests/nextflow.config b/subworkflows/local/tests/nextflow.config new file mode 100644 index 00000000..e9eed735 --- /dev/null +++ b/subworkflows/local/tests/nextflow.config @@ -0,0 +1,14 @@ +// conf/modules.config's DMR selectors don't apply in this nf-test, so set ext.* here +process { + withName: 'MODKIT_PILEUP' { + ext.args = '--phased --modified-bases 5mC 5hmC' + } + withName: '.*HTSLIB_BGZIPTABIX' { + // Avoid an hp1/hp2 filename collision + ext.prefix = { "${meta.id}_${meta.haplotype}" } + ext.args2 = '-p bed' + } + withName: '.*MODKIT_DMR' { + ext.args = '--base C' + } +} diff --git a/tests/default.nf.test.snap b/tests/default.nf.test.snap index f6f60bbc..810824dd 100644 --- a/tests/default.nf.test.snap +++ b/tests/default.nf.test.snap @@ -611,7 +611,7 @@ "sample1_whatshap_stats.tsv:md5,cf8917be389dfef1344eeb6b7e99c7c5", "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50", - "breakpoints_double.csv:md5,003e4e80ec4233315762233e5f1ec0d8", + "breakpoints_double.csv:md5,e307fa52370d89e72ea09f84d4eff2f7", "read_qual.txt:md5,92321a03aae1168fdce8ea921b635036", "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50", @@ -630,7 +630,7 @@ "sample2_whatshap_stats.tsv:md5,31f4d9c3228d60ac463695ca19d05f3a", "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50", - "breakpoints_double.csv:md5,ab7ccd32508029426ca1365053fdc23b", + "breakpoints_double.csv:md5,fa0c1072b74d3dcf914901f73154bbad", "read_qual.txt:md5,57aaac5f18f81b1cbde13e49c012651a", "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50", @@ -644,7 +644,7 @@ "sample3_whatshap_stats.tsv:md5,f7cc79156f23e884ead18e50b8434dbf", "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50", - "breakpoints_double.csv:md5,c18b66eebc1714279e685ab1c56d6cf3", + "breakpoints_double.csv:md5,0c3021e35f1e8b0ae78fd2783e678c9f", "read_qual.txt:md5,f6d803b54f3aa4a36f21af1a9a97923f", "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50" diff --git a/workflows/lrsomatic.nf b/workflows/lrsomatic.nf index 6fd130a8..e6f56c9e 100644 --- a/workflows/lrsomatic.nf +++ b/workflows/lrsomatic.nf @@ -61,6 +61,7 @@ include { PAIRED_SMALLVAR_GERMLINE } from '../subworkflows/local/paired/p include { PHASING_HAPLOTYPING } from '../subworkflows/local/phasing_haplotyping' include { TUMORONLY_SAVANA } from '../subworkflows/local/tumor_only/tumoronly_savana' include { PAIRED_SAVANA } from '../subworkflows/local/paired/paired_savana' +include { DMR } from '../subworkflows/local/dmr' include { PADFOOT_ANNOTATION } from '../subworkflows/local/padfoot_annotation' include { RECONPLOT_FIGURES } from '../subworkflows/local/reconplot_figures' @@ -114,6 +115,9 @@ workflow LRSOMATIC { params.bed_file = getGenomeAttribute('bed_file') params.savana_contigs = getGenomeAttribute('savana_contigs') params.savana_g1000_vcf = getGenomeAttribute('savana_g1000_vcf') + // Locals, not params: reassigning a param declared in nextflow.config is a no-op + def dmr_cpg_islands_bed = params.dmr_cpg_islands_bed ?: getGenomeAttribute('dmr_cpg_islands_bed') + def dmr_gencode_gene_bed = params.dmr_gencode_gene_bed ?: getGenomeAttribute('dmr_gencode_gene_bed') params.vep_genome = getGenomeAttribute('vep_genome') params.vep_species = getGenomeAttribute('vep_species') params.sigprofiler_genome = getGenomeAttribute('sigprofiler_genome') @@ -698,6 +702,31 @@ workflow LRSOMATIC { : ch_index_minimap // ch_modkit_input: [meta, bam, bai] -- BAM to pile up; meta.type selects the publish directory MODKIT_PILEUP(ch_modkit_input, ch_fasta, ch_fai, [[:],[]]) + + // + // SUBWORKFLOW: DMR + // hp1 vs hp2 methylation; BEDs default only for GRCh38/CHM13, otherwise warn and skip + if (params.modkit_phased && !params.skip_dmr) { + if (dmr_cpg_islands_bed && dmr_gencode_gene_bed) { + DMR ( + MODKIT_PILEUP.out.bedgz, + ch_fasta, + ch_fai, + file(dmr_cpg_islands_bed, checkIfExists: true), + file(dmr_gencode_gene_bed, checkIfExists: true) + ) + ch_versions = ch_versions.mix(DMR.out.versions) + } + else { + def genome_desc = params.genome ? "--genome '${params.genome}'" : "a custom --fasta with no --genome" + log.warn( + "--modkit_phased is set but DMR is being skipped: --dmr_cpg_islands_bed/" + + "--dmr_gencode_gene_bed have no default for ${genome_desc} " + + "(only GRCh38/CHM13 ship one). Pass both explicitly to enable DMR calling, " + + "or set --skip_dmr to silence this warning." + ) + } + } } // Prepare phased VCFs for VEP: add empty 'extra' list required by ENSEMBLVEP_VEP