diff --git a/.github/workflows/dryrun.yml b/.github/workflows/dryrun.yml index 071d820..22de860 100644 --- a/.github/workflows/dryrun.yml +++ b/.github/workflows/dryrun.yml @@ -22,7 +22,9 @@ jobs: matrix: NXF_VER: - "23.04.1" - - "latest" + # pinned instead of "latest": the strict syntax of Nextflow >=26.04 + # rejects the top-level statements in clean.nf, see issue #125 + - "25.10.0" inputtype: - fasta - illumina diff --git a/.gitignore b/.gitignore index 0e71847..f43be68 100644 --- a/.gitignore +++ b/.gitignore @@ -12,3 +12,4 @@ conda/ singularity/ .vscode/ .nf-test/ +.nf-test-conda/ diff --git a/CHANGELOG.md b/CHANGELOG.md index eb93954..44f86a4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,21 @@ ## unreleased +### Changed + +- replaced `bwa mem` with [`bwa-mem2`](https://github.com/bwa-mem2/bwa-mem2) as short-read mapper alternative; the `--bwa` parameter is unchanged + - faster index building, which mainly helps for large indices combining several eukaryotic genomes + - note that a `bwa-mem2` index is about 4x larger than a `bwa` index and is held in memory during mapping, so the memory requirements went up + +- all containers and conda environments now ship the same `samtools`/`htslib` version (1.24); this also updates `minimap2` 2.26 -> 2.31, `bedtools` 2.30.0 -> 2.31.1, `seqkit` 2.6.1 -> 2.13.0 and `pigz` 2.3.4 -> 2.8 + +### Fixed + +- the short-read mapper alternative now also gets an environment with the `conda`/`mamba` profiles and CPUs/memory with the `local`/`standard` profiles +- `samclipy` gets its own container (the `samtools` one, it has `python` and `git`) instead of implicitly using the one of the `smallTask` label, and `git` was added to its conda environment +- syntax that the strict parser of Nextflow >=25.10 rejects: the variable declaration in `nextflow.config`, typed `for` loops, `if` blocks around `publishDir`/`storeDir` directives, `env(VAR)` outputs and `addParams()` on `include` statements + - the CI matrix is pinned to Nextflow 25.10.0 instead of `latest`, because the top-level statements in `clean.nf` still need to move into the entry workflow before the pipeline runs on Nextflow >=26.04 + ## [v1.1.0] - 2024-11-08 ### Added diff --git a/CITATIONS.md b/CITATIONS.md index 641d26b..858b56c 100644 --- a/CITATIONS.md +++ b/CITATIONS.md @@ -12,9 +12,9 @@ > Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010 Mar 15;26(6):841-2. doi: 10.1093/bioinformatics/btq033. Epub 2010 Jan 28. PubMed PMID: 20110278; PubMed Central PMCID: PMC2832824. -- [BWA](https://arxiv.org/abs/1303.3997) +- [BWA-MEM2](https://ieeexplore.ieee.org/document/8820962) - > Li H. (2013) Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv:1303.3997v2 [q-bio.GN] + > Vasimuddin Md, Misra S, Li H, Aluru S. (2019) Efficient Architecture-Aware Acceleration of BWA-MEM for Multicore Systems. IEEE International Parallel and Distributed Processing Symposium (IPDPS). doi: 10.1109/IPDPS.2019.00041 - [FastQC](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) diff --git a/README.md b/README.md index e2cf8d9..61df55a 100644 --- a/README.md +++ b/README.md @@ -15,7 +15,7 @@ Technologies ([DNA CS (DCS)](https://assets.ctfassets.net/hkzaxo8a05x5/2IX56YmF5 ## What this workflow does for you -With this workflow you can screen and clean your Illumina, Nanopore, PacBio CLR or any FASTA-formated sequence data. The results are the clean sequences and the sequences identified as contaminated. Per default [minimap2](https://github.com/lh3/minimap2) is used for aligning your sequences to reference sequences (with the `map-ont` settings for Nanopore data, `map-bp` for PacBio CLR data, and `sr` settings for short-read data activated automatically). However, for short-read data, you may want to switch to [BWA](https://github.com/lh3/bwa) (`--bwa`). As another alternative, we provide `bbduk`, part of [BBTools](https://github.com/BioInfoTools/BBMap), as a kmer-based approach (`--bbduk`). However, no mapping file will be produced with `bbduk` and thus some subsequent statistics are not calculated. +With this workflow you can screen and clean your Illumina, Nanopore, PacBio CLR or any FASTA-formated sequence data. The results are the clean sequences and the sequences identified as contaminated. Per default [minimap2](https://github.com/lh3/minimap2) is used for aligning your sequences to reference sequences (with the `map-ont` settings for Nanopore data, `map-bp` for PacBio CLR data, and `sr` settings for short-read data activated automatically). However, for short-read data, you may want to switch to [BWA-MEM2](https://github.com/bwa-mem2/bwa-mem2) (`--bwa`). As another alternative, we provide `bbduk`, part of [BBTools](https://github.com/BioInfoTools/BBMap), as a kmer-based approach (`--bbduk`). However, no mapping file will be produced with `bbduk` and thus some subsequent statistics are not calculated. You can simply specify provided hosts and controls for the cleanup or use your own FASTA files. The reads are then mapped (or kmer-based compared in case of `bbduk`) against the specified host, control, and user defined FASTA files. All reads that match are considered as contamination. In case of Illumina paired-end reads, both mates need to be aligned (singleton files will be produced otherwise). @@ -169,7 +169,7 @@ results/ The most important files you are likely interested in are `results/clean/.fastq.gz`, which are the "cleaned" reads. These are the input reads that *do not* map to the host, control, own fasta or rRNA files (or the subset of these that you provided), plus those reads that map to the "keep" sequence if you used the `--keep` option. Any files that were removed from your input fasta file are placed in `results/removed/.fastq.gz`. -For debugging purposes we also provide various intermediate results in the `intermediate/` folder. For mapping-based approaches (`minimap2`, `bwa`), you will also find a brief summary of mapped/unmapped reads and their proportions. +For debugging purposes we also provide various intermediate results in the `intermediate/` folder. For mapping-based approaches (`minimap2`, `bwa-mem2`), you will also find a brief summary of mapped/unmapped reads and their proportions. ## Acknowledgements diff --git a/clean.nf b/clean.nf index a062fb5..86ed56c 100755 --- a/clean.nf +++ b/clean.nf @@ -10,7 +10,7 @@ Author: hoelzer.martin@gmail.com // Parameters sanity checking -Set valid_params = ['max_cores', 'cores', 'max_memory', 'memory', 'profile', 'help', 'input', 'input_type', 'list', 'host', 'own', 'control', 'keep', 'rm_rrna', 'bwa', 'bbduk', 'bbduk_kmer', 'bbduk_qin', 'reads_rna', 'min_clip', 'dcs_strict', 'output', 'multiqc_dir', 'nf_runinfo_dir', 'databases', 'cleanup_work_dir','condaCacheDir', 'singularityCacheDir', 'singularityCacheDir', 'cloudProcess', 'conda-cache-dir', 'singularity-cache-dir', 'cloud-process', 'publish_dir_mode', 'no_intermediate', 'skip_qc'] // don't ask me why there is also 'conda-cache-dir', 'singularity-cache-dir', 'cloud-process' +Set valid_params = ['max_cores', 'cores', 'max_memory', 'memory', 'profile', 'help', 'input', 'input_type', 'list', 'host', 'own', 'control', 'keep', 'rm_rrna', 'bwa', 'bbduk', 'bbduk_kmer', 'bbduk_qin', 'reads_rna', 'min_clip', 'dcs_strict', 'output', 'multiqc_dir', 'nf_runinfo_dir', 'databases', 'cleanup_work_dir','condaCacheDir', 'singularityCacheDir', 'singularityCacheDir', 'cloudProcess', 'conda-cache-dir', 'singularity-cache-dir', 'cloud-process', 'publish_dir_mode', 'no_intermediate', 'skip_qc', 'trace_timestamp'] // don't ask me why there is also 'conda-cache-dir', 'singularity-cache-dir', 'cloud-process' def parameter_diff = params.keySet() - valid_params if (parameter_diff.size() != 0){ exit 1, "ERROR: Parameter(s) $parameter_diff is/are not valid in the pipeline!\n" @@ -87,9 +87,9 @@ if ( params.input == '' || !params.input_type == '' ) { exit 1, "Missing require if ( params.input_type ) { if ( ! (params.input_type in input_types ) ) { exit 1, "Choose one of the the input types with --input_type: " + input_types } } -if ( params.control ) { for( String ctr : params.control.split(',') ) if ( ! (ctr in controls ) ) { exit 1, "Wrong control defined (" + ctr + "), use one of these: " + controls } } +if ( params.control ) { params.control.split(',').each { ctr -> if ( ! (ctr in controls ) ) { exit 1, "Wrong control defined (" + ctr + "), use one of these: " + controls } } } if ( params.input_type == 'nano' && params.control && 'dcs' in params.control.split(',') && 'eno' in params.control.split(',') ) { exit 1, "Please choose either eno (for ONT dRNA-Seq) or dcs (for ONT DNA-Seq)." } -if ( params.host ) { for( String hst : params.host.split(',') ) if ( ! (hst in hosts ) ) { exit 1, "Wrong host defined (" + hst + "), use one of these: " + hosts } } +if ( params.host ) { params.host.split(',').each { hst -> if ( ! (hst in hosts ) ) { exit 1, "Wrong host defined (" + hst + "), use one of these: " + hosts } } } if ( !params.host && !params.own && !params.control && !params.rm_rrna ) { exit 1, "Please provide a control (--control), a host tag (--host), a FASTA file (--own) or set --rm_rrna for rRNA removal for the clean up."} /************************** @@ -116,6 +116,16 @@ if ( params.input_type == 'illumina' ) { } } +// Without this an empty input channel leaves every process without a task and +// the run reports success without having cleaned anything. checkIfExists does +// not catch it: the glob can match files and still not pair them up. +input_ch = input_ch.ifEmpty { + def hint = ( params.input_type == 'illumina' && !params.list ) + ? " Paired-end reads are collected with fromFilePairs, so the glob needs the read pair group, e.g. '*_R{1,2}.fastq.gz'. For single reads use --input_type illumina_single_end." + : '' + error "No input reads found for --input '${params.input}'.${hint}" +} + // load control fasta sequence if ( params.control ) { if ( 'phix' in params.control.split(',') ) { @@ -172,17 +182,18 @@ if ( params.keep && params.list ) { multiqc_config = Channel.fromPath( workflow.projectDir + '/assets/multiqc_config.yml', checkIfExists: true ) -tool = params.bbduk ? 'bbduk' : 'minimap2' -lib_pairedness = params.input_type == 'illumina' ? 'paired' : 'single' +// `addParams()` on include statements was removed with the strict syntax of +// Nextflow >=25.10, so the modules read this from the params scope directly +params.lib_pairedness = params.input_type == 'illumina' ? 'paired' : 'single' /************************** * MODULES **************************/ -include { prepare_contamination } from './workflows/prepare_contamination_wf' addParams( tool: tool ) +include { prepare_contamination } from './workflows/prepare_contamination_wf' include { check_own as prepare_keep } from './modules/prepare_contamination' -include { clean } from './workflows/clean_wf' addParams( tool: tool, lib_pairedness: lib_pairedness ) -include { keep } from './workflows/keep_wf' addParams( tool: tool, lib_pairedness: lib_pairedness ) +include { clean } from './workflows/clean_wf' +include { keep } from './workflows/keep_wf' include { summarize } from './workflows/summarize_wf' include { qc } from './workflows/qc_wf' @@ -235,7 +246,7 @@ def helpMSG() { Clean your Illumina, Nanopore, PacBio or any FASTA-formated sequence date. The output are the clean and as contaminated identified sequences. Per default minimap2 is used for aligning your sequences - to a host but we recommend using BWA for mapping short reads ${c_dim}--bwa${c_reset} or the ${c_dim}--bbduk${c_reset} flag + to a host but we recommend using BWA-MEM2 for mapping short reads ${c_dim}--bwa${c_reset} or the ${c_dim}--bbduk${c_reset} flag to switch to bbduk to clean short-read data. Use the ${c_dim}--host${c_reset} and ${c_dim}--control${c_reset} flag to download a host database or specify your ${c_dim}--own${c_reset} FASTA. @@ -277,7 +288,7 @@ def helpMSG() { Reads are assigned to a combined index for decontamination and keeping. The use of this parameter can prevent false positive hits and the accidental removal of reads due to (poor quality) mappings. ${c_green}--rm_rrna ${c_reset} Clean your data from rRNA [default: $params.rm_rrna] - ${c_green}--bwa${c_reset} Add this flag to use BAW MEM instead of minimap2 for decontamination of short reads [default: $params.bwa] + ${c_green}--bwa${c_reset} Add this flag to use BWA-MEM2 instead of minimap2 for decontamination of short reads [default: $params.bwa] ${c_green}--bbduk${c_reset} Add this flag to use bbduk instead of minimap2 for decontamination of short reads [default: $params.bbduk] ${c_green}--bbduk_kmer${c_reset} Set kmer for bbduk [default: $params.bbduk_kmer] ${c_green}--bbduk_qin${c_reset} Set quality ASCII encoding for bbduk [default: $params.bbduk_qin; options are: 64, 33, auto] diff --git a/configs/conda.config b/configs/conda.config index d13d30d..19406b3 100644 --- a/configs/conda.config +++ b/configs/conda.config @@ -1,6 +1,9 @@ process { withLabel: minimap2 { conda = "$baseDir/envs/minimap2.yaml" } + withLabel: samtools { conda = "$baseDir/envs/samtools.yaml" } withLabel: bbmap { conda = "$baseDir/envs/bbmap.yaml" } + withLabel: bwamem2 { conda = "$baseDir/envs/bwamem2.yaml" } + withLabel: bwamem2_index { conda = "$baseDir/envs/bwamem2.yaml" } withLabel: pysam { conda = "$baseDir/envs/pysam.yaml" } withLabel: samclipy { conda = "$baseDir/envs/samclipy.yaml" } withLabel: multiqc { conda = "$baseDir/envs/multiqc.yaml" } diff --git a/configs/container.config b/configs/container.config index 0f746c8..52f0b5e 100644 --- a/configs/container.config +++ b/configs/container.config @@ -1,12 +1,15 @@ process { - withLabel: smallTask { container = 'nanozoo/samtools:1.14--d8fb865' } - withLabel: minimap2 { container = 'nanozoo/minimap2:2.26--d9ef6b6' } - withLabel: bwa { container = 'nanozoo/bwa:0.7.18--0eff897' } + withLabel: smallTask { container = 'nanozoo/samtools:1.24--62fe65d' } + withLabel: samclipy { container = 'nanozoo/samtools:1.24--62fe65d' } + withLabel: minimap2 { container = 'nanozoo/minimap2:2.31--c2b4c91' } + withLabel: samtools { container = 'nanozoo/samtools:1.24--62fe65d' } + withLabel: bwamem2 { container = 'nanozoo/bwa-mem2:2.3--a42e075' } + withLabel: bwamem2_index { container = 'nanozoo/bwa-mem2:2.3--a42e075' } withLabel: bbmap { container = 'nanozoo/bbmap:38.79--8e915d7' } withLabel: multiqc { container = 'nanozoo/multiqc:1.9--aba729b' } withLabel: fastqc { container = 'nanozoo/fastqc:0.11.9--f61b8b4' } withLabel: nanoplot { container = 'nanozoo/nanoplot:1.32.0--1ae6f5d' } withLabel: quast { container = 'nanozoo/quast:5.0.2--e7f0cfe' } - withLabel: bed_samtools { container = 'nanozoo/bed_samtools:2.30.0--cc7d1b9' } - withLabel: seqkit { container = 'nanozoo/seqkit:2.6.1--022e008' } + withLabel: bed_samtools { container = 'nanozoo/bed_samtools:2.31.1--852feef' } + withLabel: seqkit { container = 'nanozoo/seqkit:2.13.0--22ab5cc' } } diff --git a/configs/local.config b/configs/local.config index 6b92feb..2d99d9c 100644 --- a/configs/local.config +++ b/configs/local.config @@ -1,6 +1,9 @@ process { withLabel: minimap2 { cpus = params.cores } + withLabel: samtools { cpus = 1 } withLabel: bbmap { cpus = params.cores ; memory = params.memory } + withLabel: bwamem2 { cpus = params.cores ; memory = params.memory } + withLabel: bwamem2_index { cpus = 1 ; memory = params.memory } withLabel: samclipy { cpus = 1 } withLabel: pysam { cpus = 1 } withLabel: fastqc { cpus = 2 } diff --git a/configs/node.config b/configs/node.config index c9705b5..3bafc28 100644 --- a/configs/node.config +++ b/configs/node.config @@ -1,7 +1,22 @@ +// Without a time directive every task inherits the default time limit of the +// queue, which can be far below what mapping a large sample needs; the job is +// then killed mid-run (exit 140 on SLURM) and the memory retry ladder cannot +// help because more RAM does not buy more wall time. Lower these if a queue +// rejects the request. process { - withLabel: minimap2 { cpus = 24; memory = {24.GB * task.attempt}; maxRetries = 4 ; errorStrategy = { task.exitStatus in 1 || 130..140 ? 'retry' : 'terminate' }; } - withLabel: bwa { cpus = 24; memory = {24.GB * task.attempt}; maxRetries = 3 ; errorStrategy = { task.exitStatus in 130..140 ? 'retry' : 'terminate' } } - withLabel: bbmap { cpus = 24; memory = {24.GB * task.attempt}; maxRetries = 6 ; errorStrategy = { task.exitStatus in 1 || 130..140 ? 'retry' : 'terminate' }; } + withLabel: minimap2 { cpus = 24; memory = {24.GB * task.attempt}; time = {12.h * task.attempt}; maxRetries = 4 ; errorStrategy = { task.exitStatus in 1 || 130..140 ? 'retry' : 'terminate' }; } + // samtools index/idxstats/flagstats are minutes even on a large BAM, so they + // ask for a fraction of the mapping request and do not queue behind it + withLabel: samtools { cpus = 4; memory = {8.GB * task.attempt}; time = {2.h * task.attempt}; maxRetries = 2 ; errorStrategy = { task.exitStatus in 130..140 ? 'retry' : 'terminate' } } + // bwa-mem2 holds the full index in memory, which is ~4x larger than a bwa index; + // the index build needs even more memory than the mapping and uses up to ~10 + // cores despite the largely serial suffix array construction. Multi-genome + // hosts need much more than these defaults: measured on a 8 Gbp host (human + + // grey seal + harbor seal) the build peaked at 188 GB and wrote a 42 GB + // index, so raise the memory via -c instead of paying for the retry ladder. + withLabel: bwamem2 { cpus = 24; memory = {32.GB * task.attempt}; time = {12.h * task.attempt}; maxRetries = 4 ; errorStrategy = { task.exitStatus in 130..140 ? 'retry' : 'terminate' } } + withLabel: bwamem2_index { cpus = 8; memory = {48.GB * task.attempt}; time = {12.h * task.attempt}; maxRetries = 4 ; errorStrategy = { task.exitStatus in 130..140 ? 'retry' : 'terminate' } } + withLabel: bbmap { cpus = 24; memory = {24.GB * task.attempt}; time = {12.h * task.attempt}; maxRetries = 6 ; errorStrategy = { task.exitStatus in 1 || 130..140 ? 'retry' : 'terminate' }; } withLabel: smallTask { cpus = 1; memory = 2.GB } withLabel: pysam { cpus = 2; memory = 4.GB } withLabel: fastqc { cpus = 2; memory = {4.GB * task.attempt } ; maxRetries = 3 ; errorStrategy = { task.exitStatus in 130..140 ? 'retry' : 'terminate' }; } diff --git a/configs/test.config b/configs/test.config index beced6e..d1406b9 100644 --- a/configs/test.config +++ b/configs/test.config @@ -9,4 +9,24 @@ params { cores = 1 max_memory = 6.GB memory = 6.GB + + // One shared environment cache for the whole test suite. The default is + // relative to the launch dir, which is a fresh directory per nf-test test, + // so every environment was built again for every single test. + condaCacheDir = "$baseDir/.nf-test-conda" +} + +// The `conda` profile sets this from params.condaCacheDir while nextflow.config +// is parsed, which is before this file is merged, so it has to be set again +conda { + cacheDir = "$baseDir/.nf-test-conda" +} + +// Build one Conda environment at a time. Parallel `conda env create` calls +// fight over the lock of the shared package cache, which fails the tests with +// `LockError: Failed to acquire lock` or a half-extracted package archive. +// The tasks themselves are stubs, so there is nothing to gain from running +// them in parallel here anyway. +executor { + queueSize = 1 } \ No newline at end of file diff --git a/envs/bed_samtools.yaml b/envs/bed_samtools.yaml index f012a93..47a482c 100644 --- a/envs/bed_samtools.yaml +++ b/envs/bed_samtools.yaml @@ -4,5 +4,6 @@ channels: - bioconda - nodefaults dependencies: - - bedtools=2.30.0 - - samtools=1.14 + - bedtools=2.31.1 + - samtools=1.24 + - htslib=1.24 diff --git a/envs/bwamem2.yaml b/envs/bwamem2.yaml new file mode 100644 index 0000000..1c3bc1f --- /dev/null +++ b/envs/bwamem2.yaml @@ -0,0 +1,9 @@ +name: bwamem2 +channels: + - conda-forge + - bioconda + - nodefaults +dependencies: + - bwa-mem2=2.3 + - samtools=1.24 + - htslib=1.24 diff --git a/envs/minimap2.yaml b/envs/minimap2.yaml index 69ed488..0d7c163 100644 --- a/envs/minimap2.yaml +++ b/envs/minimap2.yaml @@ -4,6 +4,7 @@ channels: - bioconda - nodefaults dependencies: - - minimap2=2.26 - - samtools=1.17 - - pigz=2.3.4 + - minimap2=2.31 + - samtools=1.24 + - htslib=1.24 + - pigz=2.8 diff --git a/envs/samclipy.yaml b/envs/samclipy.yaml index 6e428b9..2b2c1d9 100644 --- a/envs/samclipy.yaml +++ b/envs/samclipy.yaml @@ -4,5 +4,7 @@ channels: - bioconda - nodefaults dependencies: - - samtools=1.11 - - python=3.9.5 + - samtools=1.24 + - htslib=1.24 + - python=3.14 + - git diff --git a/envs/bwa.yaml b/envs/samtools.yaml similarity index 53% rename from envs/bwa.yaml rename to envs/samtools.yaml index f85da32..2a6ba42 100644 --- a/envs/bwa.yaml +++ b/envs/samtools.yaml @@ -1,9 +1,8 @@ -name: bwa +name: samtools channels: - conda-forge - bioconda - nodefaults dependencies: - - bwa=0.7.18 - - samtools=1.20 - - htslib=1.20 + - samtools=1.24 + - htslib=1.24 diff --git a/envs/seqkit.yaml b/envs/seqkit.yaml index 1b042a3..8f46965 100644 --- a/envs/seqkit.yaml +++ b/envs/seqkit.yaml @@ -4,6 +4,6 @@ channels: - bioconda - nodefaults dependencies: - - seqkit==2.6.1 - - tabix==1.11 - - samtools==1.18 + - seqkit=2.13.0 + - samtools=1.24 + - htslib=1.24 diff --git a/modules/alignment_processing.nf b/modules/alignment_processing.nf index 06bdc82..f0c414c 100644 --- a/modules/alignment_processing.nf +++ b/modules/alignment_processing.nf @@ -140,12 +140,12 @@ process fastq_from_bam { // When using --keep, cleaned fastq files are generated by the // filter_fastq_by_name process and not here - if ( !params.keep ) { - publishDir ( - path: params.output, - mode: params.publish_dir_mode, - pattern: "*.gz", - saveAs: { fn -> + publishDir ( + path: params.output, + mode: params.publish_dir_mode, + pattern: "*.gz", + enabled: !params.keep, + saveAs: { fn -> fn.matches('.*.unmapped.fast[aq].gz$') ? "clean/${fn}".replaceAll(~'.unmapped(.fast[aq].gz)$', '$1') : fn.matches('.*.mapped.fast[aq].gz$') ? "removed/${fn}".replaceAll(~'.mapped(.fast[aq].gz)$', '$1') : fn.matches('.*.unmapped_merged.fast[aq].gz$') ? "clean/${fn}".replaceAll(~'.unmapped_merged(.fast[aq].gz)$', '$1') : @@ -157,9 +157,8 @@ process fastq_from_bam { fn.matches('.*.unmapped_merged_(1|2|singleton).fast[aq].gz$') ? "clean/${fn}".replaceAll(~'.unmapped_(1|2|singleton)(.fast[aq].gz)$', '_$1$2') : fn.matches('.*.mapped_merged_(1|2|singleton).fast[aq].gz$') ? "removed/${fn}".replaceAll(~'.mapped_(1|2|singleton)(.fast[aq].gz)$', '_$1$2') : fn - } - ) - } + } + ) input: tuple val(name), val(type), path(bam) @@ -188,7 +187,7 @@ process fastq_from_bam { } process idxstats_from_bam { - label 'minimap2' + label 'samtools' publishDir ( path: "${params.output}/intermediate", @@ -217,7 +216,7 @@ process idxstats_from_bam { } process flagstats_from_bam { - label 'minimap2' + label 'samtools' publishDir ( path: "${params.output}/intermediate", @@ -266,7 +265,7 @@ process sort_bam { } process index_bam { - label 'minimap2' + label 'samtools' publishDir ( path: "${params.output}/intermediate", diff --git a/modules/bbmap.nf b/modules/bbmap.nf index 46c56f9..99f17cc 100644 --- a/modules/bbmap.nf +++ b/modules/bbmap.nf @@ -10,18 +10,17 @@ process bbduk { // When using `--keep`, we need to do further processing before we have // the final clean and removed data sets. - if ( !params.keep ) { - publishDir ( - path: params.output, - mode: params.publish_dir_mode, - pattern: "*.gz", - saveAs: { fn -> - fn.endsWith('.clean.fastq.gz') ? "clean/${fn}".replaceAll(~'.fastq.clean.fastq.gz$', '.fastq.gz') : - fn.endsWith('.contamination.fastq.gz') ? "removed/${fn}".replaceAll(~'.fastq.contamination.fastq.gz$', '.fastq.gz') : - fn - } - ) - } + publishDir ( + path: params.output, + mode: params.publish_dir_mode, + pattern: "*.gz", + enabled: !params.keep, + saveAs: { fn -> + fn.endsWith('.clean.fastq.gz') ? "clean/${fn}".replaceAll(~'.fastq.clean.fastq.gz$', '.fastq.gz') : + fn.endsWith('.contamination.fastq.gz') ? "removed/${fn}".replaceAll(~'.fastq.contamination.fastq.gz$', '.fastq.gz') : + fn + } + ) input: tuple val(name), path(reads) diff --git a/modules/bwa.nf b/modules/bwamem2.nf similarity index 64% rename from modules/bwa.nf rename to modules/bwamem2.nf index d6b5179..43fce58 100644 --- a/modules/bwa.nf +++ b/modules/bwamem2.nf @@ -1,33 +1,33 @@ -process bwa_index { - label 'bwa' - +process bwamem2_index { + label 'bwamem2_index' + input: path(fasta) - + output: - path(bwa) , emit: index - + path('bwamem2') , emit: index + script: """ - mkdir bwa - bwa \\ + mkdir bwamem2 + bwa-mem2 \\ index \\ - -p bwa/db \\ + -p bwamem2/db \\ $fasta """ - + stub: """ - mkdir bwa - - touch bwa/db.{amb,ann,bwt,pac,sa} + mkdir bwamem2 + + touch bwamem2/db.{0123,amb,ann,bwt.2bit.64,pac} """ } -process bwa { - label 'bwa' +process bwamem2 { + label 'bwamem2' - input: + input: tuple val(name), path(input) path(db_index) path(db) @@ -39,7 +39,7 @@ process bwa { script: """ INDEX=`find -L ./ -name "*.amb" | sed 's/\\.amb\$//'` - bwa mem \\ + bwa-mem2 mem \\ -t $task.cpus \\ \$INDEX \\ $input \\ diff --git a/modules/minimap2.nf b/modules/minimap2.nf index 3d3d575..5d2cf9d 100644 --- a/modules/minimap2.nf +++ b/modules/minimap2.nf @@ -1,7 +1,44 @@ +// -k and -w are baked into the index, so the index has to be built with the +// same preset that is later used for the mapping +def minimap2_preset() { + if ( params.input_type == 'nano' ) { + return params.reads_rna ? '-x splice -k14' : '-x map-ont' + } else if ( params.input_type == 'pacbio' ) { + return params.reads_rna ? '-x splice -k14' : '-x map-pb' + } else if ( params.input_type.contains('illumina') ) { + return '-x sr' + } else if ( params.input_type == 'fasta' ) { + return '-x asm5' + } else { + error "Unknown input_type: ${params.input_type}" + } +} + +process minimap2_index { + label 'minimap2' + + input: + path (fasta) + + output: + path 'db.mmi' + + script: + // without this every sample rebuilds the index, which dominates the runtime + // for large hosts + """ + minimap2 ${minimap2_preset()} -t ${task.cpus} -d db.mmi ${fasta} + """ + stub: + """ + touch db.mmi + """ +} + process minimap2 { label 'minimap2' - input: + input: tuple val(name), path(input) path (db) @@ -10,27 +47,10 @@ process minimap2 { script: // -N is an internal algorithm option. It controls how many candidates alignment to extend. --secondary is an output option. - if ( params.input_type == 'nano' ) { - params = params.reads_rna ? "-ax splice -k14" : "-ax map-ont" - """ - minimap2 ${params} -N 5 --split-prefix tmp --secondary=no -t ${task.cpus} ${db} ${input} | samtools view -bhS -@ ${task.cpus} > ${name}.bam - """ - } else if ( params.input_type == 'pacbio' ) { - params = params.reads_rna ? "-ax splice -k14" : "-ax map-pb" - """ - minimap2 ${params} -N 5 --split-prefix tmp --secondary=no -t ${task.cpus} ${db} ${input} | samtools view -bhS -@ ${task.cpus} > ${name}.bam - """ - } else if ( params.input_type.contains('illumina') ) { - """ - minimap2 -ax sr -N 5 --split-prefix tmp --secondary=no -t ${task.cpus} ${db} ${input} | samtools view -bhS -@ ${task.cpus} > ${name}.bam - """ - } else if ( params.input_type == 'fasta' ){ - """ - minimap2 -ax asm5 -N 5 --split-prefix tmp --secondary=no -t ${task.cpus} ${db} ${input} | samtools view -bhS -@ ${task.cpus} > ${name}.bam - """ - } else { - error "Unknown input_type: ${params.input_type}" - } + // --split-prefix keeps the output correct when the index is split into parts + """ + minimap2 -a ${minimap2_preset()} -N 5 --split-prefix tmp --secondary=no -t ${task.cpus} ${db} ${input} | samtools view -bhS -@ ${task.cpus} > ${name}.bam + """ stub: """ touch ${name}.bam diff --git a/modules/prepare_contamination.nf b/modules/prepare_contamination.nf index dbe3b70..9f03962 100644 --- a/modules/prepare_contamination.nf +++ b/modules/prepare_contamination.nf @@ -1,12 +1,14 @@ process download_host { label 'minimap2' - if (params.cloudProcess) { - publishDir "${params.databases}/hosts", mode: params.publish_dir_mode, pattern: "*.fa.gz" - } - else { - storeDir "${params.databases}/hosts" - } + // in the cloud we cannot store the downloaded genomes, so we publish them + publishDir ( + path: "${params.databases}/hosts", + mode: params.publish_dir_mode, + pattern: "*.fa.gz", + enabled: params.cloudProcess + ) + storeDir ( params.cloudProcess ? null : "${params.databases}/hosts" ) input: val host @@ -62,15 +64,17 @@ process check_own { path fasta output: - path 'checked.fa.gz' + path "${fasta.baseName}.checked.fa.gz" script: + // the output is named after the input: with more than one --own/--keep FASTA + // a fixed name would collide when the checked files are staged together """ - seqkit seq ${fasta} -o checked.fa.gz + seqkit seq ${fasta} -o ${fasta.baseName}.checked.fa.gz """ stub: """ - touch checked.fa.gz + touch ${fasta.baseName}.checked.fa.gz """ } diff --git a/modules/utils.nf b/modules/utils.nf index 2a81130..4942b6a 100644 --- a/modules/utils.nf +++ b/modules/utils.nf @@ -5,7 +5,7 @@ process get_number_of_records { tuple val(name), path(reads) output: - tuple val(name), env(TOTALRECORDS), emit: TOTALRECORDS + tuple val(name), env('TOTALRECORDS'), emit: TOTALRECORDS script: if ( params.lib_pairedness == 'paired' ) { @@ -92,12 +92,16 @@ process filter_fastq_by_name { // When using --keep, this is where the final cleaned fastq file is // generated. If not, it's generated by the fastq_from_bam process or // bbduk mapping - if ( params.keep ) { - publishDir ( - path: params.output, - mode: params.publish_dir_mode, - pattern: "*.gz", - saveAs: { fn -> + publishDir ( + path: params.output, + mode: params.publish_dir_mode, + pattern: "*.gz", + // `as boolean` is needed: params.keep is a file path, and publishDir runs + // the value through Boolean.parseBoolean(), which turns any other string + // into false. Elsewhere we use `!params.keep`, where `!` already gives a + // boolean, so the coercion is only needed here. + enabled: params.keep as boolean, + saveAs: { fn -> fn.matches('.*.unmapped.fast[aq].gz$') ? "clean/${fn}".replaceAll(~'.unmapped(.fast[aq].gz)$', '$1') : fn.matches('.*.mapped.fast[aq].gz$') ? "removed/${fn}".replaceAll(~'.mapped(.fast[aq].gz)$', '$1') : fn.matches('.*.unmapped_merged.fast[aq].gz$') ? "clean/${fn}".replaceAll(~'.unmapped_merged(.fast[aq].gz)$', '$1') : @@ -111,9 +115,8 @@ process filter_fastq_by_name { fn.matches('.*.unmapped_merged_(1|2|singleton).fast[aq].gz$') ? "clean/${fn}".replaceAll(~'.unmapped_(1|2|singleton)(.fast[aq].gz)$', '_$1$2') : fn.matches('.*.mapped_merged_(1|2|singleton).fast[aq].gz$') ? "removed/${fn}".replaceAll(~'.mapped_(1|2|singleton)(.fast[aq].gz)$', '_$1$2') : fn - } - ) - } + } + ) input: tuple val(name), path(keep_read_name_list), val(mapped), path(reads_mapped), val(unmapped), path(reads_unmapped) diff --git a/nextflow.config b/nextflow.config index 03d0e09..06fac10 100644 --- a/nextflow.config +++ b/nextflow.config @@ -24,7 +24,7 @@ params { control = false keep = false rm_rrna = false - bwa = false + bwa = false // use BWA-MEM2 instead of minimap2 for short reads bbduk = false bbduk_kmer = 27 bbduk_qin = 'auto' @@ -51,20 +51,23 @@ params { // cleanup cleanup_work_dir = false no_intermediate = false + + // timestamp for the timeline and report file names; a param because the + // config syntax of Nextflow >=25.10 does not allow variable declarations + trace_timestamp = new java.util.Date().format('yyyy-MM-dd_HH-mm-ss') } // see https://www.nextflow.io/docs/latest/config.html?highlight=cleanup#miscellaneous cleanup = params.cleanup_work_dir -def trace_timestamp = new java.util.Date().format('yyyy-MM-dd_HH-mm-ss') timeline { enabled = true - file = "${params.output}/${params.nf_runinfo_dir}/execution_timeline_${trace_timestamp}.html" + file = "${params.output}/${params.nf_runinfo_dir}/execution_timeline_${params.trace_timestamp}.html" } report { enabled = true - file = "${params.output}/${params.nf_runinfo_dir}/execution_report_${trace_timestamp}.html" + file = "${params.output}/${params.nf_runinfo_dir}/execution_report_${params.trace_timestamp}.html" } profiles { diff --git a/tests/illumina/main.nf.test b/tests/illumina/main.nf.test index c62b7c0..3d2fde3 100644 --- a/tests/illumina/main.nf.test +++ b/tests/illumina/main.nf.test @@ -41,8 +41,8 @@ nextflow_pipeline { } } - test("Stub paired-end bwa") { - tag "bwa" + test("Stub paired-end bwa-mem2") { + tag "bwamem2" options "-stub-run" when { @@ -61,8 +61,8 @@ nextflow_pipeline { } } - test("Stub single-end bwa") { - tag "bwa" + test("Stub single-end bwa-mem2") { + tag "bwamem2" options "-stub-run" when { diff --git a/workflows/clean_wf.nf b/workflows/clean_wf.nf index 904ce06..71f13f0 100644 --- a/workflows/clean_wf.nf +++ b/workflows/clean_wf.nf @@ -1,5 +1,5 @@ -include { minimap2 } from '../modules/minimap2' -include { bwa_index; bwa } from '../modules/bwa' +include { minimap2_index; minimap2 } from '../modules/minimap2' +include { bwamem2_index; bwamem2 } from '../modules/bwamem2' include { bbduk } from '../modules/bbmap' include { bbduk_stats } from '../modules/utils' include { split_bam; fastq_from_bam ; idxstats_from_bam ; flagstats_from_bam ; index_bam as index_bam; index_bam as index_bam2; sort_bam ; filter_true_dcs_alignments ; merge_bam as merge_bam1 ; merge_bam as merge_bam2 ; merge_bam as merge_bam3 ; merge_bam as merge_bam4 ; filter_soft_clipped_alignments } from '../modules/alignment_processing' @@ -26,11 +26,12 @@ workflow clean { } else { if ( params.bwa ) { - bwa_index(contamination) - bwa(input, bwa_index.out, contamination) | sort_bam | index_bam | ( idxstats_from_bam & flagstats_from_bam ) - split_bam(bwa.out.bam) + bwamem2_index(contamination) + bwamem2(input, bwamem2_index.out, contamination) | sort_bam | index_bam | ( idxstats_from_bam & flagstats_from_bam ) + split_bam(bwamem2.out.bam) } else { - minimap2(input, contamination) | sort_bam | index_bam | ( idxstats_from_bam & flagstats_from_bam ) + minimap2_index(contamination) + minimap2(input, minimap2_index.out) | sort_bam | index_bam | ( idxstats_from_bam & flagstats_from_bam ) split_bam(minimap2.out.bam) } contamination_bam = split_bam.out.mapped