Add bowtie2 as a mapper

README.md

@@ -89,6 +89,7 @@ and one of these for mapping:
 * [strobealign](https://github.com/ksahlin/StrobeAlign) v0.14.0
 * [minimap2](https://github.com/lh3/minimap2) v2.21
 * [bwa-mem2](https://github.com/bwa-mem2/bwa-mem2) v2.0
+* [bowtie2](https://github.com/BenLangmead/bowtie2) v2.4.0
 
 and one of these for genome dereplication:
 * [skani](https://github.com/bluenote-1577/skani) v0.1.1
@@ -119,10 +120,58 @@ CoverM operates in several modes. Detailed usage information including examples
 
 There are several utility modes as well:
 * [make](https://wwood.github.io/CoverM/coverm-make.html) - Generate BAM files through alignment
+* [makedb](https://wwood.github.io/CoverM/coverm-makedb.html) - Generate mapping database(s) from reference FASTA files
 * [filter](https://wwood.github.io/CoverM/coverm-filter.html) - Remove (or only keep) alignments with insufficient identity
 * [cluster](https://wwood.github.io/CoverM/coverm-cluster.html) - Dereplicate and cluster genomes
 * shell-completion - Generate shell completion scripts
 
+The `makedb` mode pre-generates mapping indexes (databases) from reference
+genome or contig FASTA files, so that they can be reused across multiple
+`coverm contig`/`coverm genome` runs without re-indexing each time. The kind of
+database is selected with `--mapper`, and multiple `--mapper` values may be
+given to generate several databases at once. minimap2 (all presets),
+`bwa-mem`/`bwa-mem2`, `bowtie2` and `strobealign` are supported:
+
+```bash
+# Generate a short-read minimap2 database
+coverm makedb -r combined_genomes.fna -p minimap2-sr -o db_dir
+
+# Use the generated minimap2 database when calculating coverage
+coverm contig \
+  -r db_dir/combined_genomes.fna.minimap2-sr.mmi \
+  --minimap2-reference-is-index \
+  -1 read1.fq -2 read2.fq
+
+# Generate several databases at once, one per mapper
+coverm makedb -r combined_genomes.fna -p minimap2-sr minimap2-ont bwa-mem -o db_dir
+
+# Generate a bowtie2 database
+coverm makedb -r combined_genomes.fna -p bowtie2 -o db_dir
+
+# Use the generated bowtie2 database when calculating coverage (pass the prefix)
+coverm contig \
+  -r db_dir/combined_genomes.fna.bowtie2 \
+  -p bowtie2 \
+  -1 read1.fq -2 read2.fq
+```
+
+strobealign indexes are read-length specific and require the reference FASTA at
+mapping time, so for `strobealign` the reference is copied into the output
+directory alongside the index. Set the read length with `--strobealign-params
+'-r <length>'`, or estimate it from an example read dataset by passing a reads
+file:
+
+```bash
+# Generate a strobealign database for 150bp reads
+coverm makedb -r combined_genomes.fna -p strobealign --strobealign-params '-r 150' -o db_dir
+
+# Use the generated strobealign database when calculating coverage
+coverm contig \
+  -r db_dir/combined_genomes.fna \
+  --strobealign-use-index \
+  -1 read1.fq -2 read2.fq
+```
+
 ## Demo
 
 A common use case for CoverM is to calculate the coverage or relative abundance of a set of genomes in a metagenomic sample.

pixi.toml

@@ -19,6 +19,7 @@ samtools = ">=1.9"
 coreutils = "*"
 minimap2 = ">=2.24"
 bwa = ">=0.7.17"
+bowtie2 = ">=2.4.0"
 skani = ">=0.2.2"
 fastani = ">=1.3"
 extern = "*"

release.sh

@@ -32,7 +32,7 @@ tar czf coverm-x86_64-unknown-linux-musl-$VERSION.tar.gz coverm-x86_64-unknown-l
 cd ..
 
 echo "Building HTML versions of man pages .."
-for SUBCOMMAND in genome cluster contig filter make
+for SUBCOMMAND in genome cluster contig filter make makedb
 do
     echo "Documenting $SUBCOMMAND .."
     cargo run -- $SUBCOMMAND --full-help-roff |pandoc - -t markdown -f man |sed 's/\\\[/[/g; s/\\\]/]/g' |cat <(sed s/SUBCOMMAND/$SUBCOMMAND/ prelude) - >docs/coverm-$SUBCOMMAND.Rmd

src/bam_generator.rs

@@ -54,6 +54,7 @@ pub enum MappingProgram {
     MINIMAP2_LR_HQ,
     MINIMAP2_NO_PRESET,
     STROBEALIGN,
+    BOWTIE2,
 }
 
 pub struct BamFileNamedReader {
@@ -416,7 +417,10 @@ pub fn generate_named_bam_readers_from_reads(
 
     // Required because of https://github.com/wwood/CoverM/issues/58
     let minimap2_log_file_index = match mapping_program {
-        MappingProgram::BWA_MEM | MappingProgram::BWA_MEM2 | MappingProgram::STROBEALIGN => None,
+        MappingProgram::BWA_MEM
+        | MappingProgram::BWA_MEM2
+        | MappingProgram::STROBEALIGN
+        | MappingProgram::BOWTIE2 => None,
         // Required because of https://github.com/lh3/minimap2/issues/527
         MappingProgram::MINIMAP2_SR
         | MappingProgram::MINIMAP2_ONT
@@ -857,6 +861,20 @@ impl NamedBamMaker {
     }
 }
 
+/// Peek at the first record of a reads file to determine whether it is FASTA
+/// (as opposed to FASTQ). Compression (gzip/bzip2/xz) is handled transparently
+/// by needletail. Defaults to false (treat as FASTQ) if the file cannot be
+/// parsed, leaving any resulting error to be reported by the mapper itself.
+fn reads_are_fasta(read_path: &str) -> bool {
+    match needletail::parse_fastx_file(read_path) {
+        Ok(mut reader) => match reader.next() {
+            Some(Ok(record)) => record.format() == needletail::parser::Format::Fasta,
+            _ => false,
+        },
+        Err(_) => false,
+    }
+}
+
 pub fn build_mapping_command(
     mapping_program: MappingProgram,
     read_format: ReadFormat,
@@ -866,6 +884,35 @@ pub fn build_mapping_command(
     read2_path: Option<&str>,
     mapping_options: Option<&str>,
 ) -> String {
+    // bowtie2 has a sufficiently different command-line structure (it takes the
+    // index via -x, the reads via -1/-2/-U/--interleaved, and threads via -p)
+    // that it is simplest to build its command separately.
+    if let MappingProgram::BOWTIE2 = mapping_program {
+        let read_args = match read_format {
+            ReadFormat::Interleaved => format!("--interleaved '{read1_path}'"),
+            ReadFormat::Coupled => {
+                format!("-1 '{}' -2 '{}'", read1_path, read2_path.unwrap())
+            }
+            ReadFormat::Single => format!("-U '{read1_path}'"),
+        };
+        // Unlike minimap2/bwa/strobealign, bowtie2 does not auto-detect whether
+        // the reads are FASTA or FASTQ (it defaults to FASTQ), so pass '-f' when
+        // the reads are FASTA.
+        let format_flag = if reads_are_fasta(read1_path) {
+            "-f "
+        } else {
+            ""
+        };
+        return format!(
+            "bowtie2 {}{} -p {} -x '{}' {}",
+            format_flag,
+            mapping_options.unwrap_or(""),
+            threads,
+            reference.index_path(),
+            read_args
+        );
+    }
+
     let read_params1 = match mapping_program {
         // minimap2 auto-detects interleaved based on read names
         MappingProgram::MINIMAP2_SR
@@ -882,6 +929,7 @@ pub fn build_mapping_command(
             ReadFormat::Interleaved => "--interleaved",
             ReadFormat::Coupled | ReadFormat::Single => "",
         },
+        MappingProgram::BOWTIE2 => unreachable!(),
     };
 
     let read_params2 = match read_format {
@@ -915,7 +963,8 @@ pub fn build_mapping_command(
                     match mapping_program {
                         MappingProgram::BWA_MEM
                         | MappingProgram::BWA_MEM2
-                        | MappingProgram::STROBEALIGN => unreachable!(),
+                        | MappingProgram::STROBEALIGN
+                        | MappingProgram::BOWTIE2 => unreachable!(),
                         MappingProgram::MINIMAP2_SR => "-x sr",
                         MappingProgram::MINIMAP2_ONT => "-x map-ont",
                         MappingProgram::MINIMAP2_HIFI => "-x map-hifi",