Replace bwa with bwa-mem2 and align samtools versions

.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

.gitignore

@@ -12,3 +12,4 @@ conda/
 singularity/
 .vscode/
 .nf-test/
+.nf-test-conda/

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

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/)
 

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/<sample_name>.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/<sample_name>.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
 

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]

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" }

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' }
 }

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 }

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' }; }

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
 }

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

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

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

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

envs/bwa.yaml → 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

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