diff --git a/README.md b/README.md index 68d2618..4f2c946 100644 --- a/README.md +++ b/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 '`, 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. diff --git a/pixi.lock b/pixi.lock index c2df8fa..6ceaebe 100644 --- a/pixi.lock +++ b/pixi.lock @@ -8,6 +8,7 @@ environments: linux-64: - conda: https://conda.anaconda.org/conda-forge/linux-64/_libgcc_mutex-0.1-conda_forge.tar.bz2 - conda: https://conda.anaconda.org/conda-forge/linux-64/_openmp_mutex-4.5-2_gnu.tar.bz2 + - conda: https://conda.anaconda.org/bioconda/linux-64/bowtie2-2.5.4-h7071971_4.tar.bz2 - conda: https://conda.anaconda.org/bioconda/linux-64/bwa-0.7.18-he4a0461_1.tar.bz2 - conda: https://conda.anaconda.org/bioconda/linux-64/bwa-mem2-2.2.1-hd03093a_5.tar.bz2 - conda: https://conda.anaconda.org/conda-forge/linux-64/bzip2-1.0.8-h4bc722e_7.conda @@ -69,6 +70,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/linux-64/zlib-1.2.13-h4ab18f5_6.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/zstd-1.5.6-ha6fb4c9_0.conda osx-64: + - conda: https://conda.anaconda.org/bioconda/osx-64/bowtie2-2.5.5-h9566767_0.conda - conda: https://conda.anaconda.org/bioconda/osx-64/bwa-0.7.18-h10309d6_1.tar.bz2 - conda: https://conda.anaconda.org/conda-forge/osx-64/bzip2-1.0.8-hfdf4475_7.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/c-ares-1.34.5-hf13058a_0.conda @@ -78,7 +80,7 @@ environments: - conda: https://conda.anaconda.org/bioconda/osx-64/fastani-1.34-ha9a6269_4.tar.bz2 - conda: https://conda.anaconda.org/conda-forge/osx-64/gsl-2.7-h93259b0_0.tar.bz2 - conda: https://conda.anaconda.org/bioconda/osx-64/htslib-1.21-hec81eee_0.tar.bz2 - - conda: https://conda.anaconda.org/bioconda/osx-64/k8-1.2-hdf58011_0.tar.bz2 + - conda: https://conda.anaconda.org/bioconda/osx-64/k8-1.2-h2ec61ea_6.tar.bz2 - conda: https://conda.anaconda.org/conda-forge/osx-64/krb5-1.21.3-h37d8d59_0.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/libblas-3.9.0-32_h7f60823_openblas.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/libcblas-3.9.0-32_hff6cab4_openblas.conda @@ -96,14 +98,13 @@ environments: - conda: https://conda.anaconda.org/conda-forge/osx-64/libopenblas-0.3.30-openmp_hbf64a52_0.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/libsqlite-3.46.0-h1b8f9f3_0.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/libssh2-1.11.0-hd019ec5_0.conda - - conda: https://conda.anaconda.org/conda-forge/osx-64/libzlib-1.2.13-h87427d6_6.conda + - conda: https://conda.anaconda.org/conda-forge/osx-64/libzlib-1.3.2-hbb4bfdb_2.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/llvm-openmp-20.1.8-hf4e0ed4_0.conda - conda: https://conda.anaconda.org/bioconda/osx-64/minimap2-2.28-h10309d6_3.tar.bz2 - conda: https://conda.anaconda.org/conda-forge/osx-64/ncurses-6.5-h0622a9a_3.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/openssl-3.5.1-hc426f3f_0.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/perl-5.32.1-7_h10d778d_perl5.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/python-3.8.19-h5ba8234_0_cpython.conda - - conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.8-8_cp38.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/readline-8.2-h7cca4af_2.conda - conda: https://conda.anaconda.org/bioconda/osx-64/samtools-1.21-h94387ee_0.tar.bz2 - conda: https://conda.anaconda.org/bioconda/osx-64/skani-0.2.2-hdc87acc_2.tar.bz2 @@ -111,9 +112,10 @@ environments: - conda: https://conda.anaconda.org/conda-forge/osx-64/xz-5.8.1-h357f2ed_2.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/xz-gpl-tools-5.8.1-h357f2ed_2.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/xz-tools-5.8.1-hd471939_2.conda - - conda: https://conda.anaconda.org/conda-forge/osx-64/zlib-1.2.13-h87427d6_6.conda - - conda: https://conda.anaconda.org/conda-forge/osx-64/zstd-1.5.6-h915ae27_0.conda + - conda: https://conda.anaconda.org/conda-forge/osx-64/zlib-1.3.2-hbb4bfdb_2.conda + - conda: https://conda.anaconda.org/conda-forge/osx-64/zstd-1.5.7-h3eecb57_6.conda osx-arm64: + - conda: https://conda.anaconda.org/bioconda/osx-arm64/bowtie2-2.5.5-h9e91881_0.conda - conda: https://conda.anaconda.org/bioconda/osx-arm64/bwa-0.7.19-hba9b596_1.tar.bz2 - conda: https://conda.anaconda.org/conda-forge/osx-arm64/bzip2-1.0.8-h99b78c6_7.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/c-ares-1.34.5-h5505292_0.conda @@ -176,6 +178,50 @@ packages: license_family: BSD size: 23621 timestamp: 1650670423406 +- conda: https://conda.anaconda.org/bioconda/linux-64/bowtie2-2.5.4-h7071971_4.tar.bz2 + sha256: d02662a31dbca71f6230065e154b3e00666f387b5083486f16b2f7b563a6f77e + md5: 69822858766e6c8b12ae90d78d54d8ea + depends: + - _openmp_mutex >=4.5 + - libgcc-ng >=12 + - libgomp + - libstdcxx-ng >=12 + - libzlib >=1.2.13,<2.0a0 + - perl + - python + - zstd >=1.5.6,<1.6.0a0 + license: GPL-3.0-only + license_family: GPL3 + size: 14966002 + timestamp: 1724162485205 +- conda: https://conda.anaconda.org/bioconda/osx-64/bowtie2-2.5.5-h9566767_0.conda + sha256: 88c1113ff9a0636f22dd82ea7065bfc89d0baaf0d44614fa65dfa15102070436 + md5: d7afdb37f545f01cfc53013b380e9c1b + depends: + - libcxx >=19 + - libzlib >=1.3.1,<2.0a0 + - llvm-openmp >=19.1.7 + - perl + - python + - zstd >=1.5.7,<1.6.0a0 + license: GPL-3.0-or-later + license_family: GPL3 + size: 942370 + timestamp: 1771304854228 +- conda: https://conda.anaconda.org/bioconda/osx-arm64/bowtie2-2.5.5-h9e91881_0.conda + sha256: d8d48ada756702685b36e63c27e79463e5ee11aea31f84239e11e1a52dba6e9b + md5: 0381ecdbc2697f14bbff354c82edc92d + depends: + - libcxx >=19 + - libzlib >=1.3.1,<2.0a0 + - llvm-openmp >=19.1.7 + - perl + - python + - zstd >=1.5.7,<1.6.0a0 + license: GPL-3.0-or-later + license_family: GPL3 + size: 842128 + timestamp: 1771302725933 - conda: https://conda.anaconda.org/bioconda/linux-64/bwa-0.7.18-he4a0461_1.tar.bz2 sha256: f8a9482e07214cd1595f021ff66c327e4d77a903515d11295a7f5db05d3fdb86 md5: 4ecde7d4a03fabe440a90aea1389d62d @@ -459,17 +505,16 @@ packages: license_family: MIT size: 7561524 timestamp: 1727909086631 -- conda: https://conda.anaconda.org/bioconda/osx-64/k8-1.2-hdf58011_0.tar.bz2 - sha256: 5e689f331993af7be0822456d957693778818ac30bfdf9696228ee91dd09a149 - md5: 955639d6044f3d9946c9e2de78ca75e5 +- conda: https://conda.anaconda.org/bioconda/osx-64/k8-1.2-h2ec61ea_6.tar.bz2 + sha256: 596c13ecad0e679a907c229b9fcd07f649235ad25abbf783e8910a449ec60f2d + md5: 917525300d56b7a7ddac9e1bde2dd4ba depends: - - libcxx >=14 - - libzlib >=1.2.13,<1.3.0a0 - - python_abi 3.8.* *_cp38 + - libcxx >=18 + - libzlib >=1.3.1,<2.0a0 license: MIT license_family: MIT - size: 8445584 - timestamp: 1716846952727 + size: 8484801 + timestamp: 1748948933576 - conda: https://conda.anaconda.org/bioconda/osx-arm64/k8-1.2-hda5e58c_6.tar.bz2 sha256: 71a5f0bcdb892760556b741fcba8b52d9824a414097cc5e877f36e89b093e3c0 md5: 7564d9ace7ab6eed5f972bf2a799b166 @@ -1194,17 +1239,17 @@ packages: license_family: Other size: 61571 timestamp: 1716874066944 -- conda: https://conda.anaconda.org/conda-forge/osx-64/libzlib-1.2.13-h87427d6_6.conda - sha256: 1c70fca0720685242b5c68956f310665c7ed43f04807aa4227322eee7925881c - md5: c0ef3c38a80c02ae1d86588c055184fc +- conda: https://conda.anaconda.org/conda-forge/osx-64/libzlib-1.3.2-hbb4bfdb_2.conda + sha256: 4c6da089952b2d70150c74234679d6f7ac04f4a98f9432dec724968f912691e7 + md5: 30439ff30578e504ee5e0b390afc8c65 depends: - - __osx >=10.13 + - __osx >=11.0 constrains: - - zlib 1.2.13 *_6 + - zlib 1.3.2 *_2 license: Zlib license_family: Other - size: 57373 - timestamp: 1716874185419 + size: 59000 + timestamp: 1774073052242 - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libzlib-1.3.1-h8359307_2.conda sha256: ce34669eadaba351cd54910743e6a2261b67009624dbc7daeeafdef93616711b md5: 369964e85dc26bfe78f41399b366c435 @@ -1435,16 +1480,6 @@ packages: license_family: BSD size: 7002 timestamp: 1752805902938 -- conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.8-8_cp38.conda - build_number: 8 - sha256: 83c22066a672ce0b16e693c84aa6d5efb68e02eff037a55e047d7095d0fdb5ca - md5: 4f7b6e3de4f15cc44e0f93b39f07205d - constrains: - - python 3.8.* *_cpython - license: BSD-3-Clause - license_family: BSD - size: 6960 - timestamp: 1752805923703 - conda: https://conda.anaconda.org/conda-forge/linux-64/readline-8.2-h8c095d6_2.conda sha256: 2d6d0c026902561ed77cd646b5021aef2d4db22e57a5b0178dfc669231e06d2c md5: 283b96675859b20a825f8fa30f311446 @@ -1674,16 +1709,16 @@ packages: license_family: Other size: 92883 timestamp: 1716874088980 -- conda: https://conda.anaconda.org/conda-forge/osx-64/zlib-1.2.13-h87427d6_6.conda - sha256: 3091d48a579c08ba20885bc8856def925e9dee9d1a7d8713e3ce002eb29fcd19 - md5: 700b922d6d22e7deb5fb2964d0c8cf6a +- conda: https://conda.anaconda.org/conda-forge/osx-64/zlib-1.3.2-hbb4bfdb_2.conda + sha256: 5dd728cebca2e96fa48d41661f1a35ed0ee3cb722669eee4e2d854c6745655eb + md5: 6276aa61ffc361cbf130d78cfb88a237 depends: - - __osx >=10.13 - - libzlib 1.2.13 h87427d6_6 + - __osx >=11.0 + - libzlib 1.3.2 hbb4bfdb_2 license: Zlib license_family: Other - size: 88732 - timestamp: 1716874218187 + size: 92411 + timestamp: 1774073075482 - conda: https://conda.anaconda.org/conda-forge/linux-64/zstd-1.5.6-ha6fb4c9_0.conda sha256: c558b9cc01d9c1444031bd1ce4b9cff86f9085765f17627a6cd85fc623c8a02b md5: 4d056880988120e29d75bfff282e0f45 @@ -1695,16 +1730,16 @@ packages: license_family: BSD size: 554846 timestamp: 1714722996770 -- conda: https://conda.anaconda.org/conda-forge/osx-64/zstd-1.5.6-h915ae27_0.conda - sha256: efa04a98cb149643fa54c4dad5a0179e36a5fbc88427ea0eec88ceed87fd0f96 - md5: 4cb2cd56f039b129bb0e491c1164167e +- conda: https://conda.anaconda.org/conda-forge/osx-64/zstd-1.5.7-h3eecb57_6.conda + sha256: 47101a4055a70a4876ffc87b750ab2287b67eca793f21c8224be5e1ee6394d3f + md5: 727109b184d680772e3122f40136d5ca depends: - - __osx >=10.9 - - libzlib >=1.2.13,<2.0.0a0 + - __osx >=10.13 + - libzlib >=1.3.1,<2.0a0 license: BSD-3-Clause license_family: BSD - size: 498900 - timestamp: 1714723303098 + size: 528148 + timestamp: 1764777156963 - conda: https://conda.anaconda.org/conda-forge/osx-arm64/zstd-1.5.7-h6491c7d_2.conda sha256: 0d02046f57f7a1a3feae3e9d1aa2113788311f3cf37a3244c71e61a93177ba67 md5: e6f69c7bcccdefa417f056fa593b40f0 diff --git a/pixi.toml b/pixi.toml index d4a4526..7ea1a5c 100644 --- a/pixi.toml +++ b/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 = "*" diff --git a/release.sh b/release.sh index 89daeaf..19cce60 100755 --- a/release.sh +++ b/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 diff --git a/src/bam_generator.rs b/src/bam_generator.rs index 5c1b953..9d7cdac 100644 --- a/src/bam_generator.rs +++ b/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", diff --git a/src/bin/coverm.rs b/src/bin/coverm.rs index 1fa08e8..49b7ce4 100644 --- a/src/bin/coverm.rs +++ b/src/bin/coverm.rs @@ -715,6 +715,147 @@ fn main() { } } } + Some("makedb") => { + let m = matches.subcommand_matches("makedb").unwrap(); + bird_tool_utils::clap_utils::print_full_help_if_needed(m, makedb_full_help()); + set_log_level(m, true); + manually_check_args_at_runtime(m); + + let output_directory = m.get_one::("output-directory").unwrap(); + setup_bam_cache_directory(output_directory); + + let references: Vec<&String> = m + .get_many::("reference") + .expect("No reference provided") + .collect(); + + // Validate that each reference is an existing FASTA file. We + // deliberately do not use check_reference_existence here: for BWA + // it inspects pre-existing index files beside the reference, which + // is inappropriate when makedb is *creating* a new index (and would + // spuriously fail when, say, a bwa-mem2 index already sits next to a + // FASTA from which a bwa-mem database is being built). + for reference in &references { + let ref_path = std::path::Path::new(reference.as_str()); + if !ref_path.exists() { + error!("The reference specified '{reference}' does not appear to exist"); + process::exit(1); + } else if !ref_path.is_file() { + error!( + "The reference specified '{reference}' should be a file, \ + not e.g. a directory" + ); + process::exit(1); + } + } + + // Guard against multiple references that share a file name (e.g. + // refs in different directories both named ref.fna), which would + // otherwise generate databases with colliding output paths. + if references.len() > 1 { + let mut seen_stems = HashSet::new(); + for reference in &references { + let stem = std::path::Path::new(reference.as_str()) + .file_name() + .expect("Failed to glean file name from reference path") + .to_string_lossy() + .to_string(); + if !seen_stems.insert(stem.clone()) { + error!( + "Multiple references share the file name '{stem}', which would \ + generate databases with colliding output paths. Please rename or \ + generate them in separate output directories." + ); + process::exit(1); + } + } + } + + // Parse mappers, de-duplicating while preserving order, and check + // that the underlying mapping software is installed. + let mut mapping_programs = vec![]; + let mut seen_mappers = HashSet::new(); + for mapper_name in m.get_many::("mapper").expect("No mapper provided") { + if !seen_mappers.insert(mapper_name.clone()) { + warn!("Ignoring duplicate --mapper value {mapper_name}"); + continue; + } + let mapping_program = mapping_program_from_name(Some(mapper_name)); + check_mapping_program_dependencies(mapping_program); + mapping_programs.push(mapping_program); + } + + let num_threads = *m.get_one::("threads").unwrap(); + + let mut generated_dbs = vec![]; + for &mapping_program in &mapping_programs { + let index_creation_params = match mapping_program { + MappingProgram::BWA_MEM | MappingProgram::BWA_MEM2 => { + m.get_one::("bwa-params") + } + MappingProgram::MINIMAP2_SR + | MappingProgram::MINIMAP2_ONT + | MappingProgram::MINIMAP2_PB + | MappingProgram::MINIMAP2_HIFI + | MappingProgram::MINIMAP2_LR_HQ + | MappingProgram::MINIMAP2_NO_PRESET => m.get_one::("minimap2-params"), + MappingProgram::STROBEALIGN => m.get_one::("strobealign-params"), + MappingProgram::BOWTIE2 => m.get_one::("bowtie2-params"), + }; + for reference in &references { + let db_path = coverm::mapping_index_maintenance::generate_persistent_index( + mapping_program, + reference, + output_directory, + Some(num_threads), + index_creation_params.map(|x| x.as_str()), + ); + info!("Generated {mapping_program:?} database at {db_path}"); + generated_dbs.push((mapping_program, db_path)); + } + } + + info!("Finished generating {} database(s).", generated_dbs.len()); + for (mapping_program, db_path) in &generated_dbs { + match mapping_program { + MappingProgram::MINIMAP2_SR + | MappingProgram::MINIMAP2_ONT + | MappingProgram::MINIMAP2_PB + | MappingProgram::MINIMAP2_HIFI + | MappingProgram::MINIMAP2_LR_HQ + | MappingProgram::MINIMAP2_NO_PRESET => { + info!( + "To use the minimap2 database, run e.g.: coverm contig \ + --reference {db_path} --minimap2-reference-is-index -1 read1.fq -2 read2.fq" + ); + } + MappingProgram::BWA_MEM => { + info!( + "To use the BWA database, run e.g.: coverm contig \ + --reference {db_path} -p bwa-mem -1 read1.fq -2 read2.fq" + ); + } + MappingProgram::BWA_MEM2 => { + info!( + "To use the BWA-MEM2 database, run e.g.: coverm contig \ + --reference {db_path} -p bwa-mem2 -1 read1.fq -2 read2.fq" + ); + } + MappingProgram::STROBEALIGN => { + info!( + "To use the strobealign database, run e.g.: coverm contig \ + --reference {db_path} --strobealign-use-index -1 read1.fq -2 read2.fq" + ); + } + MappingProgram::BOWTIE2 => { + info!( + "To use the bowtie2 database, run e.g.: coverm contig \ + --reference {db_path} -p bowtie2 -1 read1.fq -2 read2.fq" + ); + } + } + } + } Some("shell-completion") => { let m = matches.subcommand_matches("shell-completion").unwrap(); set_log_level(m, true); @@ -800,6 +941,11 @@ fn setup_mapping_index( ) } } + MappingProgram::BOWTIE2 => coverm::mapping_index_maintenance::generate_bowtie2_index( + reference_wise_params.reference, + Some(*m.get_one::("threads").unwrap()), + m.get_one::("bowtie2-params").map(|x| x.as_str()), + ), MappingProgram::STROBEALIGN => { // Indexing once for a batch of readsets is not yet supported for strobealign info!("Not pre-generating strobealign index"); @@ -871,7 +1017,13 @@ fn dereplicate(m: &clap::ArgMatches, genome_fasta_files: &Vec) -> Vec MappingProgram { - let mapping_program = match m.get_one::("mapper").map(|x| &**x) { + let mapping_program = mapping_program_from_name(m.get_one::("mapper").map(|x| &**x)); + check_mapping_program_dependencies(mapping_program); + mapping_program +} + +fn mapping_program_from_name(name: Option<&str>) -> MappingProgram { + match name { Some("bwa-mem") => MappingProgram::BWA_MEM, Some("bwa-mem2") => MappingProgram::BWA_MEM2, Some("minimap2-sr") => MappingProgram::MINIMAP2_SR, @@ -881,12 +1033,13 @@ fn parse_mapping_program(m: &clap::ArgMatches) -> MappingProgram { Some("minimap2-lr-hq") => MappingProgram::MINIMAP2_LR_HQ, Some("minimap2-no-preset") => MappingProgram::MINIMAP2_NO_PRESET, Some("strobealign") => MappingProgram::STROBEALIGN, + Some("bowtie2") => MappingProgram::BOWTIE2, None => DEFAULT_MAPPING_SOFTWARE_ENUM, - _ => panic!( - "Unexpected definition for --mapper: {:?}", - m.get_one::("mapper") - ), - }; + _ => panic!("Unexpected definition for --mapper: {:?}", name), + } +} + +fn check_mapping_program_dependencies(mapping_program: MappingProgram) { match mapping_program { MappingProgram::BWA_MEM => { external_command_checker::check_for_bwa(); @@ -905,8 +1058,10 @@ fn parse_mapping_program(m: &clap::ArgMatches) -> MappingProgram { MappingProgram::STROBEALIGN => { external_command_checker::check_for_strobealign(); } + MappingProgram::BOWTIE2 => { + external_command_checker::check_for_bowtie2(); + } } - mapping_program } struct EstimatorsAndTaker { @@ -1494,7 +1649,7 @@ fn setup_bam_cache_directory(cache_directory: &str) { ); process::exit(1); } else { - info!("Writing BAM files to already existing directory {cache_directory}") + info!("Writing output files to already existing directory {cache_directory}") } } else { match path.parent() { diff --git a/src/cli.rs b/src/cli.rs index 1793e36..ade1670 100644 --- a/src/cli.rs +++ b/src/cli.rs @@ -21,6 +21,7 @@ const MAPPING_SOFTWARE_LIST: &[&str] = &[ "minimap2-lr-hq", "minimap2-no-preset", "strobealign", + "bowtie2", ]; const DEFAULT_MAPPING_SOFTWARE: &str = "strobealign"; @@ -94,6 +95,10 @@ fn add_mapping_options(manual: Manual) -> Manual { &monospace_roff("minimap2-no-preset"), &format!("minimap2 with no '{}' option", &monospace_roff("-x")) ], + &[ + &monospace_roff("bowtie2"), + "bowtie2 using default parameters" + ], ]) ))) .option(Opt::new("PARAMS").long("--minimap2-params").help(&format!( @@ -119,6 +124,13 @@ fn add_mapping_options(manual: Manual) -> Manual { that usage of this parameter has security \ implications if untrusted input is specified. \ [default: none]", + )) + .option(Opt::new("PARAMS").long("--bowtie2-params").help( + "Extra parameters to provide to bowtie2, both the \ + bowtie2-build indexing command (if used) and \ + the mapping command. Note that usage of this \ + parameter has security implications if untrusted \ + input is specified. [default: none]", )) .flag(Flag::new().long("--strobealign-use-index").help( "Use a pregenerated index (one that has been created with 'strobealign --create-index'). The --reference option should be specified as the original FASTA file i.e. 'ref.fna' not 'ref.fna.r100.sti' [default: not set]", @@ -478,6 +490,173 @@ pub fn make_full_help() -> Manual { manual } +pub fn makedb_full_help() -> Manual { + let mut manual = Manual::new("coverm makedb") + .about(format!( + "Generate mapping database(s) from reference FASTA files (version: {})", + crate_version!() + )) + .custom_synopsis_expansion("-r -p .. -o ") + .author(Author::new(crate::AUTHOR).email("benjwoodcroft near gmail.com")) + .description( + "coverm makedb pre-generates one or more mapping databases (indexes) from \ + reference genome or contig FASTA files. The generated database can then be \ + supplied to 'coverm contig' or 'coverm genome' (or 'coverm make') to avoid \ + regenerating the index each time reads are mapped.\n\n\ + For minimap2 databases, pass the generated '.mmi' file as the reference \ + together with '--minimap2-reference-is-index'. For BWA databases, pass the \ + generated prefix as the reference together with the matching '-p bwa-mem' or \ + '-p bwa-mem2'. For bowtie2 databases, pass the generated prefix as the reference \ + together with '-p bowtie2'. For strobealign databases, the reference FASTA is copied into the \ + output directory next to the index (strobealign reads the sequences from it at \ + mapping time); pass that copied FASTA as the reference together with \ + '--strobealign-use-index'.\n\n\ + Multiple '-p/--mapper' values may be specified to create several databases in \ + one invocation, one per mapper. Likewise, multiple references may be given, in \ + which case a database is created for each combination of reference and mapper.\n\n", + ); + + manual = manual.custom(Section::new("Input").option( + Opt::new("PATH ..").short("-r").long("--reference").help( + "FASTA file(s) of contigs e.g. concatenated genomes or metagenome assembly. \ + May be gzip-compressed. [required]", + ), + )); + + manual = + manual.custom( + Section::new("Database type") + .option( + Opt::new("NAME ..") + .short("-p") + .long("--mapper") + .help(&format!( + "Kind(s) of database to generate, one per mapping software. Specify more \ + than once (or as a space-separated list) to generate several databases. \ + {}. One of: {}", + default_roff("minimap2-sr"), + bird_tool_utils::clap_utils::table_roff(&[ + &["name", "description"], + &[ + &monospace_roff("minimap2-sr"), + &format!("minimap2 index built with '{}'", &monospace_roff("-x sr")) + ], + &[ + &monospace_roff("minimap2-lr-hq"), + &format!("minimap2 index built with '{}'", &monospace_roff("-x lr:hq")) + ], + &[ + &monospace_roff("minimap2-ont"), + &format!("minimap2 index built with '{}'", &monospace_roff("-x map-ont")) + ], + &[ + &monospace_roff("minimap2-pb"), + &format!("minimap2 index built with '{}'", &monospace_roff("-x map-pb")) + ], + &[ + &monospace_roff("minimap2-hifi"), + &format!("minimap2 index built with '{}'", &monospace_roff("-x map-hifi")) + ], + &[ + &monospace_roff("minimap2-no-preset"), + &format!("minimap2 index built with no '{}' option", &monospace_roff("-x")) + ], + &[&monospace_roff("bwa-mem"), "BWA index (bwa index)"], + &[&monospace_roff("bwa-mem2"), "BWA-MEM2 index (bwa-mem2 index)"], + &[ + &monospace_roff("bowtie2"), + "bowtie2 index (bowtie2-build). The generated prefix is used \ + as the reference together with '-p bowtie2'" + ], + &[ + &monospace_roff("strobealign"), + &format!( + "strobealign index (strobealign --create-index). The reference \ + FASTA is copied into the output directory alongside the index, \ + since strobealign requires it at mapping time. Use it via '{}'", + &monospace_roff("--strobealign-use-index") + ) + ], + ]) + )), + ) + .option(Opt::new("PARAMS").long("--minimap2-params").help( + "Extra parameters to provide to the minimap2 indexing command. \ + Note that usage of this parameter has security implications if \ + untrusted input is specified. [default: none]", + )) + .option(Opt::new("PARAMS").long("--bwa-params").help( + "Extra parameters to provide to the BWA or BWA-MEM2 indexing \ + command. Note that usage of this parameter has security \ + implications if untrusted input is specified. [default: none]", + )) + .option(Opt::new("PARAMS").long("--strobealign-params").help( + "Extra parameters to provide to the 'strobealign --create-index' \ + command. Strobealign indexes are read-length specific: set the \ + canonical read length with e.g. '-r 150', or estimate it from an \ + example read dataset by passing a reads file (e.g. 'reads.fq'). \ + Note that usage of this parameter has security implications if \ + untrusted input is specified. [default: none]", + )) + .option(Opt::new("PARAMS").long("--bowtie2-params").help( + "Extra parameters to provide to the 'bowtie2-build' indexing \ + command. Note that usage of this parameter has security \ + implications if untrusted input is specified. [default: none]", + )), + ); + + manual = manual.custom(Section::new("Output").option( + Opt::new("DIR").short("-o").long("--output-directory").help( + "Where the generated database(s) will be written. The directory will \ + be created if it does not exist. [required]", + ), + )); + + manual = manual.example( + Example::new() + .text("Generate a short-read minimap2 database from a set of genomes") + .command("coverm makedb -r combined_genomes.fna -p minimap2-sr -o db_dir"), + ); + manual = manual.example( + Example::new() + .text("Use the generated minimap2 database when calculating contig coverage") + .command( + "coverm contig -r db_dir/combined_genomes.fna.minimap2-sr.mmi \ + --minimap2-reference-is-index -1 read1.fq -2 read2.fq", + ), + ); + manual = manual.example( + Example::new() + .text("Generate several databases at once, one per mapper") + .command("coverm makedb -r combined_genomes.fna -p minimap2-sr minimap2-ont bwa-mem -o db_dir"), + ); + manual = manual.example( + Example::new() + .text( + "Generate a strobealign database for 150bp reads (the reference FASTA is \ + copied into db_dir alongside the index)", + ) + .command( + "coverm makedb -r combined_genomes.fna -p strobealign \ + --strobealign-params '-r 150' -o db_dir", + ), + ); + + let mut general_section = Section::new("General options").option( + Opt::new("INT").short("-t").long("--threads").help(&format!( + "Number of threads used to generate the database(s). {}", + default_roff("1") + )), + ); + general_section = add_help_options_to_section(general_section); + general_section = add_verbosity_flags_to_section(general_section); + manual = manual.custom(general_section); + + manual = manual.custom(faq_section()); + + manual +} + pub fn contig_full_help() -> Manual { let mut manual = Manual::new("coverm contig") .about(format!("Calculate read coverage per-contig (version {})",crate_version!())) @@ -1023,6 +1202,25 @@ See coverm make --full-help for further options and further detail. storing sorted BAM files in output_dir/" ), ); + static ref MAKEDB_HELP: String = format!( + " + {} + {} + +{} + + coverm makedb -r combined_genomes.fna -p minimap2-sr -o db_dir + +See coverm makedb --full-help for further options and further detail. +", + ansi_term::Colour::Green.paint("coverm makedb"), + ansi_term::Colour::Green + .paint("Generate mapping database(s) from reference FASTA files"), + ansi_term::Colour::Purple.paint( + "Example: Generate a short-read minimap2 database from a set of genomes,\n\ + storing it in db_dir/" + ), + ); } let mut app = Command::new("coverm") @@ -1045,6 +1243,7 @@ Main subcommands: Less used utility subcommands: \tmake\tGenerate BAM files through alignment +\tmakedb\tGenerate mapping database(s) from reference FASTA files \tfilter\tRemove (or only keep) alignments with insufficient identity \tcluster\tDereplicate and cluster genomes \tshell-completion @@ -1238,6 +1437,13 @@ Ben J. Woodcroft .allow_hyphen_values(true) .requires("reference"), ) + .arg( + Arg::new("bowtie2-params") + .long("bowtie2-params") + .alias("bowtie2-parameters") + .allow_hyphen_values(true) + .requires("reference"), + ) .arg( Arg::new("strobealign-use-index") .long("strobealign-use-index") @@ -1773,6 +1979,13 @@ Ben J. Woodcroft .allow_hyphen_values(true) .requires("reference"), ) + .arg( + Arg::new("bowtie2-params") + .long("bowtie2-params") + .alias("bowtie2-parameters") + .allow_hyphen_values(true) + .requires("reference"), + ) .arg( Arg::new("strobealign-use-index") .long("strobealign-use-index") @@ -2151,6 +2364,13 @@ Ben J. Woodcroft .allow_hyphen_values(true) .requires("reference"), ) + .arg( + Arg::new("bowtie2-params") + .long("bowtie2-params") + .alias("bowtie2-parameters") + .allow_hyphen_values(true) + .requires("reference"), + ) .arg( Arg::new("strobealign-use-index") .long("strobealign-use-index") @@ -2158,6 +2378,75 @@ Ben J. Woodcroft .action(clap::ArgAction::SetTrue), ), ) + .subcommand( + add_clap_verbosity_flags(Command::new("makedb")) + .about("Generate mapping database(s) from reference FASTA files") + .override_help(MAKEDB_HELP.as_str()) + .arg( + Arg::new("full-help") + .long("full-help") + .action(clap::ArgAction::SetTrue), + ) + .arg( + Arg::new("full-help-roff") + .long("full-help-roff") + .action(clap::ArgAction::SetTrue), + ) + .arg( + Arg::new("reference") + .short('r') + .long("reference") + .action(clap::ArgAction::Append) + .num_args(1..) + .required_unless_present_any(["full-help", "full-help-roff"]), + ) + .arg( + Arg::new("output-directory") + .short('o') + .long("output-directory") + .required_unless_present_any(["full-help", "full-help-roff"]), + ) + .arg( + Arg::new("mapper") + .short('p') + .long("mapper") + .action(clap::ArgAction::Append) + .num_args(1..) + .value_parser(MAPPING_SOFTWARE_LIST.iter().collect::>()) + .default_value("minimap2-sr"), + ) + .arg( + Arg::new("threads") + .short('t') + .long("threads") + .default_value("1") + .value_parser(clap::value_parser!(u16)), + ) + .arg( + Arg::new("minimap2-params") + .long("minimap2-params") + .alias("minimap2-parameters") + .allow_hyphen_values(true), + ) + .arg( + Arg::new("bwa-params") + .long("bwa-params") + .alias("bwa-parameters") + .allow_hyphen_values(true), + ) + .arg( + Arg::new("strobealign-params") + .long("strobealign-params") + .alias("strobealign-parameters") + .allow_hyphen_values(true), + ) + .arg( + Arg::new("bowtie2-params") + .long("bowtie2-params") + .alias("bowtie2-parameters") + .allow_hyphen_values(true), + ), + ) .subcommand( add_clap_verbosity_flags(Command::new("shell-completion")) .about("Generate a shell completion script for coverm") diff --git a/src/external_command_checker.rs b/src/external_command_checker.rs index 54976ce..34ee9b8 100644 --- a/src/external_command_checker.rs +++ b/src/external_command_checker.rs @@ -24,6 +24,14 @@ pub fn check_for_samtools() { .expect("Failed to find sufficient version of samtools"); } +pub fn check_for_bowtie2() { + check_for_external_command_presence_with_which("bowtie2") + .expect("Failed to find installed bowtie2"); + // bowtie2-build is a separate executable used to generate the index. + check_for_external_command_presence_with_which("bowtie2-build") + .expect("Failed to find installed bowtie2-build"); +} + pub fn check_for_strobealign() { check_for_external_command_presence_with_which("strobealign") .expect("Failed to find installed strobealign"); diff --git a/src/mapping_index_maintenance.rs b/src/mapping_index_maintenance.rs index 64f6c3f..87c2a1e 100644 --- a/src/mapping_index_maintenance.rs +++ b/src/mapping_index_maintenance.rs @@ -57,104 +57,224 @@ impl TemporaryIndexStruct { .expect("Failed to glean file stem from reference DB. Strange."), ); - info!("Generating {mapping_program:?} index for {reference_path} .."); - let mut cmd = match mapping_program { - MappingProgram::BWA_MEM => std::process::Command::new("bwa"), - MappingProgram::BWA_MEM2 => std::process::Command::new("bwa-mem2"), - MappingProgram::MINIMAP2_SR - | MappingProgram::MINIMAP2_ONT - | MappingProgram::MINIMAP2_PB - | MappingProgram::MINIMAP2_HIFI - | MappingProgram::MINIMAP2_LR_HQ - | MappingProgram::MINIMAP2_NO_PRESET => std::process::Command::new("minimap2"), - MappingProgram::STROBEALIGN => std::process::Command::new("strobealign"), - }; - match &mapping_program { - MappingProgram::BWA_MEM | MappingProgram::BWA_MEM2 => { - cmd.arg("index") - .arg("-p") - .arg(&index_path) - .arg(reference_path); + run_index_command( + mapping_program, + reference_path, + &index_path, + num_threads, + index_creation_options, + ); + + TemporaryIndexStruct { + index_path_internal: index_path.to_string_lossy().to_string(), + tempdir: td, + } + } +} + +/// Build the index-generation command for the given mapping program, writing +/// the resulting index to `index_path`. Returns None for STROBEALIGN, which +/// does not support standalone pre-indexing in this manner. +fn build_index_command( + mapping_program: MappingProgram, + reference_path: &str, + index_path: &Path, + num_threads: Option, + index_creation_options: Option<&str>, +) -> Option { + let mut cmd = match mapping_program { + MappingProgram::BWA_MEM => std::process::Command::new("bwa"), + MappingProgram::BWA_MEM2 => std::process::Command::new("bwa-mem2"), + MappingProgram::MINIMAP2_SR + | MappingProgram::MINIMAP2_ONT + | MappingProgram::MINIMAP2_PB + | MappingProgram::MINIMAP2_HIFI + | MappingProgram::MINIMAP2_LR_HQ + | MappingProgram::MINIMAP2_NO_PRESET => std::process::Command::new("minimap2"), + MappingProgram::BOWTIE2 => std::process::Command::new("bowtie2-build"), + MappingProgram::STROBEALIGN => { + warn!("STROBEALIGN pre-indexing is not supported currently, so skipping index generation."); + return None; + } + }; + match &mapping_program { + MappingProgram::BWA_MEM | MappingProgram::BWA_MEM2 => { + cmd.arg("index") + .arg("-p") + .arg(index_path) + .arg(reference_path); + } + MappingProgram::BOWTIE2 => { + if let Some(t) = num_threads { + cmd.arg("--threads").arg(format!("{t}")); } - MappingProgram::MINIMAP2_SR - | MappingProgram::MINIMAP2_ONT - | MappingProgram::MINIMAP2_HIFI - | MappingProgram::MINIMAP2_PB - | MappingProgram::MINIMAP2_LR_HQ - | MappingProgram::MINIMAP2_NO_PRESET => { - match &mapping_program { - MappingProgram::MINIMAP2_SR => { - cmd.arg("-x").arg("sr"); - } - MappingProgram::MINIMAP2_ONT => { - cmd.arg("-x").arg("map-ont"); - } - MappingProgram::MINIMAP2_HIFI => { - cmd.arg("-x").arg("map-hifi"); - } - MappingProgram::MINIMAP2_PB => { - cmd.arg("-x").arg("map-pb"); - } - MappingProgram::MINIMAP2_LR_HQ => { - cmd.arg("-x").arg("lr:hq"); - } - MappingProgram::MINIMAP2_NO_PRESET - | MappingProgram::BWA_MEM - | MappingProgram::BWA_MEM2 - | MappingProgram::STROBEALIGN => {} - }; - if let Some(t) = num_threads { - cmd.arg("-t").arg(format!("{t}")); + // bowtie2-build takes the reference FASTA followed by the index + // prefix it should write the .bt2 files to. + cmd.arg(reference_path).arg(index_path); + } + MappingProgram::MINIMAP2_SR + | MappingProgram::MINIMAP2_ONT + | MappingProgram::MINIMAP2_HIFI + | MappingProgram::MINIMAP2_PB + | MappingProgram::MINIMAP2_LR_HQ + | MappingProgram::MINIMAP2_NO_PRESET => { + match &mapping_program { + MappingProgram::MINIMAP2_SR => { + cmd.arg("-x").arg("sr"); } - cmd.arg("-d").arg(&index_path).arg(reference_path); - } - MappingProgram::STROBEALIGN => { - warn!("STROBEALIGN pre-indexing is not supported currently, so skipping index generation."); - } - }; - if let Some(params) = index_creation_options { - for s in params.split_whitespace() { - cmd.arg(s); + MappingProgram::MINIMAP2_ONT => { + cmd.arg("-x").arg("map-ont"); + } + MappingProgram::MINIMAP2_HIFI => { + cmd.arg("-x").arg("map-hifi"); + } + MappingProgram::MINIMAP2_PB => { + cmd.arg("-x").arg("map-pb"); + } + MappingProgram::MINIMAP2_LR_HQ => { + cmd.arg("-x").arg("lr:hq"); + } + MappingProgram::MINIMAP2_NO_PRESET + | MappingProgram::BWA_MEM + | MappingProgram::BWA_MEM2 + | MappingProgram::STROBEALIGN + | MappingProgram::BOWTIE2 => {} + }; + if let Some(t) = num_threads { + cmd.arg("-t").arg(format!("{t}")); } - }; - // Some BWA versions output log info to stdout. Ignore this. - cmd.stdout(std::process::Stdio::piped()); - cmd.stderr(std::process::Stdio::piped()); - debug!("Running DB indexing command: {cmd:?}"); - - let mut process = cmd - .spawn() - .unwrap_or_else(|_| panic!("Failed to start {:?} index process", mapping_program)); - let es = process.wait().unwrap_or_else(|_| { - panic!( - "Failed to glean exitstatus from failing {:?} index process", - mapping_program - ) - }); - if !es.success() { - error!("Error when running {mapping_program:?} index process."); - let mut err = String::new(); - process - .stderr - .unwrap_or_else(|| { - panic!( - "Failed to grab stderr from failed {:?} index process", - mapping_program - ) - }) - .read_to_string(&mut err) - .expect("Failed to read stderr into string"); - error!("The STDERR was: {err:?}"); - error!("Cannot continue after {mapping_program:?} index failed."); - process::exit(1); + cmd.arg("-d").arg(index_path).arg(reference_path); } - info!("Finished generating {mapping_program:?} index."); - TemporaryIndexStruct { - index_path_internal: index_path.to_string_lossy().to_string(), - tempdir: td, + MappingProgram::STROBEALIGN => unreachable!(), + }; + if let Some(params) = index_creation_options { + for s in params.split_whitespace() { + cmd.arg(s); } + }; + Some(cmd) +} + +/// Run the index-generation command for the given mapping program, writing the +/// resulting index to `index_path`. Exits the process on failure. +fn run_index_command( + mapping_program: MappingProgram, + reference_path: &str, + index_path: &Path, + num_threads: Option, + index_creation_options: Option<&str>, +) { + info!("Generating {mapping_program:?} index for {reference_path} .."); + let cmd = match build_index_command( + mapping_program, + reference_path, + index_path, + num_threads, + index_creation_options, + ) { + Some(cmd) => cmd, + None => return, + }; + execute_index_command(cmd, mapping_program); +} + +/// Spawn and wait on a pre-built index-generation command, exiting the process +/// on failure. +fn execute_index_command(mut cmd: std::process::Command, mapping_program: MappingProgram) { + // Some BWA versions output log info to stdout. Ignore this. + cmd.stdout(std::process::Stdio::piped()); + cmd.stderr(std::process::Stdio::piped()); + debug!("Running DB indexing command: {cmd:?}"); + + let mut process = cmd + .spawn() + .unwrap_or_else(|_| panic!("Failed to start {:?} index process", mapping_program)); + let es = process.wait().unwrap_or_else(|_| { + panic!( + "Failed to glean exitstatus from failing {:?} index process", + mapping_program + ) + }); + if !es.success() { + error!("Error when running {mapping_program:?} index process."); + let mut err = String::new(); + process + .stderr + .unwrap_or_else(|| { + panic!( + "Failed to grab stderr from failed {:?} index process", + mapping_program + ) + }) + .read_to_string(&mut err) + .expect("Failed to read stderr into string"); + error!("The STDERR was: {err:?}"); + error!("Cannot continue after {mapping_program:?} index failed."); + process::exit(1); } + info!("Finished generating {mapping_program:?} index."); } + +/// Generate a strobealign index for `reference_path` inside `output_directory`. +/// +/// Unlike minimap2/BWA, `strobealign --create-index` writes its `.sti` index +/// alongside the reference FASTA (and the FASTA is still required at mapping +/// time, since strobealign reads the sequences from it). To keep the generated +/// database self-contained within `output_directory`, the reference is first +/// copied there and the index is created next to the copy. Returns the path to +/// the copied reference, which is what should be passed to +/// `--reference .. --strobealign-use-index`. +/// +/// Strobealign indexes are read-length specific. The canonical read length can +/// be set by passing `-r ` via `index_creation_options`, or estimated +/// by passing a reads file there (e.g. `reads.fq`), matching the positional +/// argument that `strobealign --create-index` accepts. +fn generate_strobealign_persistent_index( + reference_path: &str, + output_directory: &str, + num_threads: Option, + index_creation_options: Option<&str>, +) -> String { + let reference_file_name = std::path::Path::new(reference_path) + .file_name() + .expect("Failed to glean file name from reference path"); + let copied_reference = std::path::Path::new(output_directory).join(reference_file_name); + + if copied_reference == std::path::Path::new(reference_path) { + error!( + "The strobealign database for {reference_path} would overwrite the reference \ + itself. Please choose an output directory other than the reference's directory." + ); + process::exit(1); + } + + info!( + "Copying reference {} into {} so the strobealign database is self-contained ..", + reference_path, output_directory + ); + std::fs::copy(reference_path, &copied_reference).unwrap_or_else(|e| { + panic!( + "Failed to copy reference {} into output directory: {}", + reference_path, e + ) + }); + + info!("Generating STROBEALIGN index for {reference_path} .."); + let mut cmd = std::process::Command::new("strobealign"); + if let Some(t) = num_threads { + cmd.arg("-t").arg(format!("{t}")); + } + cmd.arg("--create-index").arg(&copied_reference); + if let Some(params) = index_creation_options { + for s in params.split_whitespace() { + cmd.arg(s); + } + }; + execute_index_command(cmd, MappingProgram::STROBEALIGN); + + copied_reference.to_string_lossy().to_string() +} + impl MappingIndex for TemporaryIndexStruct { fn index_path(&self) -> &String { &self.index_path_internal @@ -201,6 +321,14 @@ pub fn check_reference_existence(reference_path: &str, mapping_program: &Mapping return; } } + MappingProgram::BOWTIE2 => { + // The reference may be given either as a FASTA file or as the prefix + // of a pre-generated bowtie2 index (in which case the prefix path + // itself is not a file), so accept either. + if check_for_bowtie2_index_existence(reference_path) { + return; + } + } MappingProgram::MINIMAP2_SR | MappingProgram::MINIMAP2_ONT | MappingProgram::MINIMAP2_HIFI @@ -241,6 +369,37 @@ pub fn generate_bwa_index( } } +/// Return true if a complete set of bowtie2 index files (either the small +/// `.bt2` or large `.bt2l` variants) already exists beside `reference_path`, +/// treating `reference_path` as the bowtie2 index prefix. +fn check_for_bowtie2_index_existence(reference_path: &str) -> bool { + let suffixes = ["1", "2", "3", "4", "rev.1", "rev.2"]; + let all_present = |extension: &str| { + suffixes.iter().all(|suffix| { + std::path::Path::new(&format!("{reference_path}.{suffix}.{extension}")).exists() + }) + }; + all_present("bt2") || all_present("bt2l") +} + +pub fn generate_bowtie2_index( + reference_path: &str, + num_threads: Option, + index_creation_parameters: Option<&str>, +) -> Box { + if check_for_bowtie2_index_existence(reference_path) { + info!("bowtie2 index appears to be complete, so going ahead and using it."); + Box::new(VanillaIndexStruct::new(reference_path)) + } else { + Box::new(TemporaryIndexStruct::new( + MappingProgram::BOWTIE2, + reference_path, + num_threads, + index_creation_parameters, + )) + } +} + pub fn generate_minimap2_index( reference_path: &str, num_threads: Option, @@ -255,6 +414,79 @@ pub fn generate_minimap2_index( )) } +/// Short name used to label the kind of database generated for a given mapping +/// program, e.g. as a filename suffix in `coverm makedb`. +pub fn mapping_program_db_name(mapping_program: MappingProgram) -> &'static str { + match mapping_program { + MappingProgram::BWA_MEM => "bwa-mem", + MappingProgram::BWA_MEM2 => "bwa-mem2", + MappingProgram::MINIMAP2_SR => "minimap2-sr", + MappingProgram::MINIMAP2_ONT => "minimap2-ont", + MappingProgram::MINIMAP2_PB => "minimap2-pb", + MappingProgram::MINIMAP2_HIFI => "minimap2-hifi", + MappingProgram::MINIMAP2_LR_HQ => "minimap2-lr-hq", + MappingProgram::MINIMAP2_NO_PRESET => "minimap2-no-preset", + MappingProgram::STROBEALIGN => "strobealign", + MappingProgram::BOWTIE2 => "bowtie2", + } +} + +/// Generate a mapping index for `reference_path` in `output_directory` that is +/// persisted on disk (unlike [`TemporaryIndexStruct`]), so it can later be fed +/// back into `coverm contig`/`coverm genome`. Returns the path to the generated +/// database. Used by `coverm makedb`. +pub fn generate_persistent_index( + mapping_program: MappingProgram, + reference_path: &str, + output_directory: &str, + num_threads: Option, + index_creation_options: Option<&str>, +) -> String { + // Strobealign writes its index next to the reference (and needs the + // reference at mapping time), so it is handled separately. + if let MappingProgram::STROBEALIGN = mapping_program { + return generate_strobealign_persistent_index( + reference_path, + output_directory, + num_threads, + index_creation_options, + ); + } + + let reference_stem = std::path::Path::new(reference_path) + .file_name() + .expect("Failed to glean file name from reference path") + .to_string_lossy(); + let db_program_name = mapping_program_db_name(mapping_program); + + // For minimap2 the index is a single .mmi file, while for BWA the supplied + // path is used as a prefix for the several files BWA generates. + let index_path = match mapping_program { + MappingProgram::MINIMAP2_SR + | MappingProgram::MINIMAP2_ONT + | MappingProgram::MINIMAP2_PB + | MappingProgram::MINIMAP2_HIFI + | MappingProgram::MINIMAP2_LR_HQ + | MappingProgram::MINIMAP2_NO_PRESET => std::path::Path::new(output_directory) + .join(format!("{reference_stem}.{db_program_name}.mmi")), + MappingProgram::BWA_MEM | MappingProgram::BWA_MEM2 | MappingProgram::BOWTIE2 => { + std::path::Path::new(output_directory) + .join(format!("{reference_stem}.{db_program_name}")) + } + MappingProgram::STROBEALIGN => unreachable!(), + }; + + run_index_command( + mapping_program, + reference_path, + &index_path, + num_threads, + index_creation_options, + ); + + index_path.to_string_lossy().to_string() +} + pub fn generate_concatenated_fasta_file(fasta_file_paths: &Vec) -> NamedTempFile { let tmpfile: NamedTempFile = Builder::new() .prefix("coverm-concatenated-fasta") diff --git a/src/mapping_parameters.rs b/src/mapping_parameters.rs index c762152..043fde4 100644 --- a/src/mapping_parameters.rs +++ b/src/mapping_parameters.rs @@ -125,6 +125,7 @@ impl<'a> MappingParameters<'a> { | MappingProgram::MINIMAP2_LR_HQ | MappingProgram::MINIMAP2_NO_PRESET => "minimap2-params", MappingProgram::STROBEALIGN => "strobealign-params", + MappingProgram::BOWTIE2 => "bowtie2-params", }; let mapping_options = match m.contains_id(mapping_parameters_arg) { true => { diff --git a/src/shard_bam_reader.rs b/src/shard_bam_reader.rs index f123852..c2d6076 100644 --- a/src/shard_bam_reader.rs +++ b/src/shard_bam_reader.rs @@ -236,13 +236,16 @@ where if !second_read_alignments[i].is_unmapped() { score += aux_as(&second_read_alignments[i]); } - if max_score.is_none() || score > max_score.unwrap() { - max_score = Some(score); - winning_indices = vec![i] - } else if score == max_score.unwrap() { - winning_indices.push(i) + match max_score { + // A loser when there was a previous winner + Some(ms) if score < ms => {} + Some(ms) if score == ms => winning_indices.push(i), + // No previous winner, or a new highest score + _ => { + max_score = Some(score); + winning_indices = vec![i] + } } - // Else a loser when there was a previous winner } } } diff --git a/tests/test_cmdline.rs b/tests/test_cmdline.rs index 50752ab..192f289 100644 --- a/tests/test_cmdline.rs +++ b/tests/test_cmdline.rs @@ -502,6 +502,218 @@ mod tests { .is_file()); } + #[test] + fn test_makedb_minimap2() { + let td = tempfile::TempDir::new().unwrap(); + Assert::main_binary() + .with_args(&[ + "makedb", + "--reference", + "tests/data/7seqs.fna", + "--mapper", + "minimap2-sr", + "--output-directory", + td.path().to_str().unwrap(), + ]) + .succeeds() + .unwrap(); + assert!(td.path().join("7seqs.fna.minimap2-sr.mmi").is_file()); + } + + #[test] + fn test_makedb_multiple_mappers() { + let td = tempfile::TempDir::new().unwrap(); + Assert::main_binary() + .with_args(&[ + "makedb", + "--reference", + "tests/data/7seqs.fna", + "--mapper", + "minimap2-sr", + "minimap2-ont", + "--output-directory", + format!("{}/db_dir", td.path().to_str().unwrap()).as_str(), + ]) + .succeeds() + .unwrap(); + assert!(td + .path() + .join("db_dir") + .join("7seqs.fna.minimap2-sr.mmi") + .is_file()); + assert!(td + .path() + .join("db_dir") + .join("7seqs.fna.minimap2-ont.mmi") + .is_file()); + } + + #[test] + fn test_makedb_then_use_as_minimap2_index() { + let td = tempfile::TempDir::new().unwrap(); + Assert::main_binary() + .with_args(&[ + "makedb", + "--reference", + "tests/data/7seqs.fna", + "--mapper", + "minimap2-sr", + "--output-directory", + td.path().to_str().unwrap(), + ]) + .succeeds() + .unwrap(); + let db_path = td.path().join("7seqs.fna.minimap2-sr.mmi"); + assert!(db_path.is_file()); + + // The generated database can be fed back into coverm contig as a + // minimap2 index. + Assert::main_binary() + .with_args(&[ + "contig", + "--coupled", + "tests/data/reads_for_seq1_and_seq2.1.fq.gz", + "tests/data/reads_for_seq1_and_seq2.2.fq.gz", + "--reference", + db_path.to_str().unwrap(), + "--minimap2-reference-is-index", + "-p", + "minimap2-sr", + "-m", + "mean", + ]) + .succeeds() + .stdout() + .contains("genome2~seq1") + .unwrap(); + } + + #[test] + fn test_makedb_strobealign_then_use_as_index() { + let td = tempfile::TempDir::new().unwrap(); + // The test reads are ~100bp, so build the index for read length 100 to + // match the canonical read length strobealign estimates at mapping time. + Assert::main_binary() + .with_args(&[ + "makedb", + "--reference", + "tests/data/7seqs.fna", + "--mapper", + "strobealign", + "--strobealign-params", + "-r 100", + "--output-directory", + td.path().to_str().unwrap(), + ]) + .succeeds() + .unwrap(); + // The reference is copied into the output directory alongside the index. + let copied_reference = td.path().join("7seqs.fna"); + assert!(copied_reference.is_file()); + assert!(td.path().join("7seqs.fna.r100.sti").is_file()); + + // The generated database can be fed back into coverm contig as a + // pregenerated strobealign index. + Assert::main_binary() + .with_args(&[ + "contig", + "--coupled", + "tests/data/reads_for_seq1_and_seq2.1.fq.gz", + "tests/data/reads_for_seq1_and_seq2.2.fq.gz", + "--reference", + copied_reference.to_str().unwrap(), + "--strobealign-use-index", + "-p", + "strobealign", + "-m", + "mean", + ]) + .succeeds() + .stdout() + .contains("genome2~seq1") + .unwrap(); + } + + #[test] + fn test_bowtie2_mapping() { + Assert::main_binary() + .with_args(&[ + "contig", + "--coupled", + "tests/data/reads_for_seq1_and_seq2.1.fq.gz", + "tests/data/reads_for_seq1_and_seq2.2.fq.gz", + "--reference", + "tests/data/7seqs.fna", + "-p", + "bowtie2", + "-m", + "mean", + ]) + .succeeds() + .stdout() + .contains("genome2~seq1") + .unwrap(); + } + + #[test] + fn test_makedb_bowtie2() { + let td = tempfile::TempDir::new().unwrap(); + Assert::main_binary() + .with_args(&[ + "makedb", + "--reference", + "tests/data/7seqs.fna", + "--mapper", + "bowtie2", + "--output-directory", + td.path().to_str().unwrap(), + ]) + .succeeds() + .unwrap(); + // bowtie2-build writes its index files using the generated prefix. + assert!(td.path().join("7seqs.fna.bowtie2.1.bt2").is_file()); + assert!(td.path().join("7seqs.fna.bowtie2.rev.2.bt2").is_file()); + } + + #[test] + fn test_makedb_then_use_as_bowtie2_index() { + let td = tempfile::TempDir::new().unwrap(); + Assert::main_binary() + .with_args(&[ + "makedb", + "--reference", + "tests/data/7seqs.fna", + "--mapper", + "bowtie2", + "--output-directory", + td.path().to_str().unwrap(), + ]) + .succeeds() + .unwrap(); + let db_prefix = td.path().join("7seqs.fna.bowtie2"); + assert!(td.path().join("7seqs.fna.bowtie2.1.bt2").is_file()); + + // The generated database can be fed back into coverm contig as a + // pregenerated bowtie2 index by passing the prefix as the reference. + Assert::main_binary() + .with_args(&[ + "contig", + "--coupled", + "tests/data/reads_for_seq1_and_seq2.1.fq.gz", + "tests/data/reads_for_seq1_and_seq2.2.fq.gz", + "--reference", + db_prefix.to_str().unwrap(), + "-p", + "bowtie2", + "-m", + "mean", + ]) + .succeeds() + .stdout() + .contains("genome2~seq1") + .unwrap(); + } + #[test] fn test_relative_abundance_all_mapped() { Assert::main_binary()