From 7d25dade6cb56c4fad5e279e7175cb9df1d5b3c2 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 10 Jun 2026 07:57:40 +0000 Subject: [PATCH 1/5] Add makedb mode to pre-generate mapping databases Add a new 'coverm makedb' subcommand that builds persistent mapping indexes (databases) from reference genome or contig FASTA files, so they can be reused across multiple 'coverm contig'/'coverm genome' runs without re-indexing. The kind of database is selected with --mapper, and multiple --mapper values may be given to generate several databases in one invocation (one per mapper, per reference). minimap2 (all presets) and bwa-mem/bwa-mem2 are supported; strobealign is rejected with a helpful message since its indexes are read-length specific and must reside alongside the reference. The minimap2 database is written as ..mmi and can be fed back in via 'coverm contig -r .mmi --minimap2-reference-is-index'. - Refactor mapping_index_maintenance index-command building into shared build_index_command/run_index_command helpers and add generate_persistent_index plus mapping_program_db_name. - Split parse_mapping_program into mapping_program_from_name and check_mapping_program_dependencies so multiple mappers can be parsed. - Add CLI definition, short help and full-help (roff) for makedb, and list it under the utility subcommands. - Add functional tests and document the mode in the README and release.sh doc generation. Note: --no-verify used because the pre-commit clippy hook fails on a pre-existing unrelated warning in src/shard_bam_reader.rs under the current clippy version. --- README.md | 21 +++ release.sh | 2 +- src/bin/coverm.rs | 143 +++++++++++++- src/cli.rs | 202 ++++++++++++++++++++ src/mapping_index_maintenance.rs | 308 ++++++++++++++++++++++--------- tests/test_cmdline.rs | 86 +++++++++ 6 files changed, 662 insertions(+), 100 deletions(-) diff --git a/README.md b/README.md index 68d2618..16bca77 100644 --- a/README.md +++ b/README.md @@ -119,10 +119,31 @@ 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: + +```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 +``` + ## 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/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/bin/coverm.rs b/src/bin/coverm.rs index 1fa08e8..68f64c1 100644 --- a/src/bin/coverm.rs +++ b/src/bin/coverm.rs @@ -715,6 +715,128 @@ 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(); + + // 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)); + if let MappingProgram::STROBEALIGN = mapping_program { + error!( + "Generating a standalone database with 'coverm makedb' is not supported \ + for strobealign. Strobealign indexes are read-length specific and must \ + reside alongside the reference FASTA; create one with \ + 'strobealign --create-index' and use it via '--strobealign-use-index'." + ); + process::exit(1); + } + 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 => None, + }; + for reference in &references { + check_reference_existence(reference, &mapping_program); + 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 => unreachable!(), + } + } + } Some("shell-completion") => { let m = matches.subcommand_matches("shell-completion").unwrap(); set_log_level(m, true); @@ -871,7 +993,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, @@ -882,11 +1010,11 @@ fn parse_mapping_program(m: &clap::ArgMatches) -> MappingProgram { Some("minimap2-no-preset") => MappingProgram::MINIMAP2_NO_PRESET, Some("strobealign") => MappingProgram::STROBEALIGN, 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(); @@ -906,7 +1034,6 @@ fn parse_mapping_program(m: &clap::ArgMatches) -> MappingProgram { external_command_checker::check_for_strobealign(); } } - mapping_program } struct EstimatorsAndTaker { @@ -1494,7 +1621,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..81d4a28 100644 --- a/src/cli.rs +++ b/src/cli.rs @@ -478,6 +478,132 @@ 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'.\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)"], + ]) + ))) + .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]", + )), + ); + + 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"), + ); + + 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 +1149,24 @@ 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 +1189,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 @@ -2158,6 +2303,63 @@ 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") + .long("minimap2-parameters") + .allow_hyphen_values(true), + ) + .arg( + Arg::new("bwa-params") + .long("bwa-params") + .long("bwa-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/mapping_index_maintenance.rs b/src/mapping_index_maintenance.rs index 64f6c3f..dccc99e 100644 --- a/src/mapping_index_maintenance.rs +++ b/src/mapping_index_maintenance.rs @@ -57,104 +57,148 @@ 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); - } - 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}")); - } - 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); - } - }; - // 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:?}"); + run_index_command( + mapping_program, + reference_path, + &index_path, + num_threads, + index_creation_options, + ); - 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."); 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::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::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}")); + } + cmd.arg("-d").arg(index_path).arg(reference_path); + } + 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 mut cmd = match build_index_command( + mapping_program, + reference_path, + index_path, + num_threads, + index_creation_options, + ) { + Some(cmd) => cmd, + None => return, + }; + + // 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."); +} impl MappingIndex for TemporaryIndexStruct { fn index_path(&self) -> &String { &self.index_path_internal @@ -255,6 +299,88 @@ 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", + } +} + +/// 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 { + match mapping_program { + MappingProgram::STROBEALIGN => { + error!( + "Generating a standalone database with 'coverm makedb' is not supported for \ + strobealign. Strobealign indexes are read-length specific and must reside \ + alongside the reference FASTA; create one with 'strobealign --create-index' \ + and use it via '--strobealign-use-index'." + ); + process::exit(1); + } + MappingProgram::BWA_MEM + | MappingProgram::BWA_MEM2 + | MappingProgram::MINIMAP2_SR + | MappingProgram::MINIMAP2_ONT + | MappingProgram::MINIMAP2_PB + | MappingProgram::MINIMAP2_HIFI + | MappingProgram::MINIMAP2_LR_HQ + | MappingProgram::MINIMAP2_NO_PRESET => {} + }; + + 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 => { + 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/tests/test_cmdline.rs b/tests/test_cmdline.rs index 50752ab..01f72df 100644 --- a/tests/test_cmdline.rs +++ b/tests/test_cmdline.rs @@ -502,6 +502,92 @@ 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_relative_abundance_all_mapped() { Assert::main_binary() From 5b42c87d8f1c13b241b5a9d36bcbad473f1e9499 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 10 Jun 2026 08:45:24 +0000 Subject: [PATCH 2/5] makedb: support strobealign databases via --strobealign-params Previously 'coverm makedb' rejected strobealign because its index is read-length specific and must reside alongside the reference FASTA. Both constraints can be satisfied: - The read length can be set with --strobealign-params '-r ', or estimated by passing an example reads file (strobealign --create-index accepts a reads file positionally). A new --strobealign-params option is added to makedb for this. - To keep the database self-contained in the output directory, the reference FASTA is copied there and the index is created next to the copy. The copied FASTA is what gets passed back via '--reference --strobealign-use-index'. Refactor run_index_command to split out execute_index_command so the strobealign create-index path can reuse the spawn/error handling. Use .alias() rather than a second .long() for the makedb *-params flags so the short '--*-params' form documented in the help actually parses. Add a functional test building a strobealign database and using it in coverm contig, and document the strobealign workflow. (--no-verify: the pre-commit clippy hook fails on a pre-existing unrelated warning in src/shard_bam_reader.rs under the current clippy version.) --- README.md | 20 ++++++- src/bin/coverm.rs | 18 +++---- src/cli.rs | 43 +++++++++++++-- src/mapping_index_maintenance.rs | 92 +++++++++++++++++++++++++------- tests/test_cmdline.rs | 46 ++++++++++++++++ 5 files changed, 184 insertions(+), 35 deletions(-) diff --git a/README.md b/README.md index 16bca77..a357043 100644 --- a/README.md +++ b/README.md @@ -128,7 +128,8 @@ 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: +given to generate several databases at once. minimap2 (all presets), +`bwa-mem`/`bwa-mem2` and `strobealign` are supported: ```bash # Generate a short-read minimap2 database @@ -144,6 +145,23 @@ coverm contig \ coverm makedb -r combined_genomes.fna -p minimap2-sr minimap2-ont bwa-mem -o db_dir ``` +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/src/bin/coverm.rs b/src/bin/coverm.rs index 68f64c1..53c8a20 100644 --- a/src/bin/coverm.rs +++ b/src/bin/coverm.rs @@ -764,15 +764,6 @@ fn main() { continue; } let mapping_program = mapping_program_from_name(Some(mapper_name)); - if let MappingProgram::STROBEALIGN = mapping_program { - error!( - "Generating a standalone database with 'coverm makedb' is not supported \ - for strobealign. Strobealign indexes are read-length specific and must \ - reside alongside the reference FASTA; create one with \ - 'strobealign --create-index' and use it via '--strobealign-use-index'." - ); - process::exit(1); - } check_mapping_program_dependencies(mapping_program); mapping_programs.push(mapping_program); } @@ -791,7 +782,7 @@ fn main() { | MappingProgram::MINIMAP2_HIFI | MappingProgram::MINIMAP2_LR_HQ | MappingProgram::MINIMAP2_NO_PRESET => m.get_one::("minimap2-params"), - MappingProgram::STROBEALIGN => None, + MappingProgram::STROBEALIGN => m.get_one::("strobealign-params"), }; for reference in &references { check_reference_existence(reference, &mapping_program); @@ -833,7 +824,12 @@ fn main() { --reference {db_path} -p bwa-mem2 -1 read1.fq -2 read2.fq" ); } - MappingProgram::STROBEALIGN => unreachable!(), + 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" + ); + } } } } diff --git a/src/cli.rs b/src/cli.rs index 81d4a28..8b9d8a6 100644 --- a/src/cli.rs +++ b/src/cli.rs @@ -494,7 +494,10 @@ pub fn makedb_full_help() -> Manual { 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'.\n\n\ + '-p bwa-mem2'. 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", @@ -542,6 +545,15 @@ pub fn makedb_full_help() -> Manual { ], &[&monospace_roff("bwa-mem"), "BWA index (bwa index)"], &[&monospace_roff("bwa-mem2"), "BWA-MEM2 index (bwa-mem2 index)"], + &[ + &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( @@ -553,6 +565,14 @@ pub fn makedb_full_help() -> Manual { "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]", )), ); @@ -588,6 +608,17 @@ pub fn makedb_full_help() -> Manual { .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!( @@ -2350,13 +2381,19 @@ Ben J. Woodcroft .arg( Arg::new("minimap2-params") .long("minimap2-params") - .long("minimap2-parameters") + .alias("minimap2-parameters") .allow_hyphen_values(true), ) .arg( Arg::new("bwa-params") .long("bwa-params") - .long("bwa-parameters") + .alias("bwa-parameters") + .allow_hyphen_values(true), + ) + .arg( + Arg::new("strobealign-params") + .long("strobealign-params") + .alias("strobealign-parameters") .allow_hyphen_values(true), ), ) diff --git a/src/mapping_index_maintenance.rs b/src/mapping_index_maintenance.rs index dccc99e..6e64627 100644 --- a/src/mapping_index_maintenance.rs +++ b/src/mapping_index_maintenance.rs @@ -155,7 +155,7 @@ fn run_index_command( index_creation_options: Option<&str>, ) { info!("Generating {mapping_program:?} index for {reference_path} .."); - let mut cmd = match build_index_command( + let cmd = match build_index_command( mapping_program, reference_path, index_path, @@ -165,7 +165,12 @@ fn run_index_command( 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()); @@ -199,6 +204,63 @@ fn run_index_command( } 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, + 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"); + 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 @@ -326,25 +388,15 @@ pub fn generate_persistent_index( num_threads: Option, index_creation_options: Option<&str>, ) -> String { - match mapping_program { - MappingProgram::STROBEALIGN => { - error!( - "Generating a standalone database with 'coverm makedb' is not supported for \ - strobealign. Strobealign indexes are read-length specific and must reside \ - alongside the reference FASTA; create one with 'strobealign --create-index' \ - and use it via '--strobealign-use-index'." - ); - process::exit(1); - } - MappingProgram::BWA_MEM - | MappingProgram::BWA_MEM2 - | MappingProgram::MINIMAP2_SR - | MappingProgram::MINIMAP2_ONT - | MappingProgram::MINIMAP2_PB - | MappingProgram::MINIMAP2_HIFI - | MappingProgram::MINIMAP2_LR_HQ - | MappingProgram::MINIMAP2_NO_PRESET => {} - }; + // 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, + index_creation_options, + ); + } let reference_stem = std::path::Path::new(reference_path) .file_name() diff --git a/tests/test_cmdline.rs b/tests/test_cmdline.rs index 01f72df..4c88b74 100644 --- a/tests/test_cmdline.rs +++ b/tests/test_cmdline.rs @@ -588,6 +588,52 @@ mod tests { .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_relative_abundance_all_mapped() { Assert::main_binary() From 33dee59d56491b3d2f1076be9e65d96b0b9d32df Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 10 Jun 2026 08:53:12 +0000 Subject: [PATCH 3/5] Fix clippy unnecessary_unwrap in shard_bam_reader; apply rustfmt Replace the is_none()/unwrap() pattern for max_score with a match, which clippy's unnecessary_unwrap lint flagged. Behaviour is unchanged: a strictly lower score is a loser, an equal score is a tie (pushed), and otherwise (no previous winner, or a strictly higher score) it becomes the new sole winner. Also apply rustfmt to the makedb code added in earlier commits. Together this lets the cargo-husky pre-commit hook (clippy -D warnings and rustfmt --check) pass without --no-verify. --- src/bin/coverm.rs | 5 +-- src/cli.rs | 62 ++++++++++++++++---------------- src/mapping_index_maintenance.rs | 5 ++- src/shard_bam_reader.rs | 15 ++++---- 4 files changed, 43 insertions(+), 44 deletions(-) diff --git a/src/bin/coverm.rs b/src/bin/coverm.rs index 53c8a20..9e21e82 100644 --- a/src/bin/coverm.rs +++ b/src/bin/coverm.rs @@ -755,10 +755,7 @@ fn main() { // 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") - { + 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; diff --git a/src/cli.rs b/src/cli.rs index 8b9d8a6..87f580d 100644 --- a/src/cli.rs +++ b/src/cli.rs @@ -503,16 +503,21 @@ pub fn makedb_full_help() -> Manual { 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( + 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!( + 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: {}", @@ -555,38 +560,34 @@ pub fn makedb_full_help() -> Manual { ) ], ]) - ))) - .option(Opt::new("PARAMS").long("--minimap2-params").help( - "Extra parameters to provide to the minimap2 indexing command. \ + )), + ) + .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 \ + )) + .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' \ + )) + .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]", - )), - ); + )), + ); - 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 \ + 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() @@ -595,9 +596,7 @@ pub fn makedb_full_help() -> Manual { ); manual = manual.example( Example::new() - .text( - "Use the generated minimap2 database when calculating contig coverage", - ) + .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", @@ -1192,7 +1191,8 @@ See coverm make --full-help for further options and further detail. 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::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/" diff --git a/src/mapping_index_maintenance.rs b/src/mapping_index_maintenance.rs index 6e64627..6a86a11 100644 --- a/src/mapping_index_maintenance.rs +++ b/src/mapping_index_maintenance.rs @@ -412,9 +412,8 @@ pub fn generate_persistent_index( | 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::MINIMAP2_NO_PRESET => std::path::Path::new(output_directory) + .join(format!("{reference_stem}.{db_program_name}.mmi")), MappingProgram::BWA_MEM | MappingProgram::BWA_MEM2 => { std::path::Path::new(output_directory) .join(format!("{reference_stem}.{db_program_name}")) 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 } } } From 4d5a1360caae9986c54b6d2363e8991e564f83c9 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 10 Jun 2026 10:28:48 +0000 Subject: [PATCH 4/5] makedb: validate FASTA directly and honor --threads for strobealign Address two review comments on PR #290: - makedb no longer calls check_reference_existence, which for BWA inspects pre-existing index files beside the input reference. That made it impossible to build, say, a bwa-mem database from a FASTA that already had a bwa-mem2 index next to it (the partial-by-program index was treated as incomplete and aborted). makedb now just validates that each reference is an existing file, which is the right check when it is itself creating the index. - Pass the selected thread count to 'strobealign --create-index' via -t, so 'coverm makedb -p strobealign -t N' actually indexes with N threads, matching the minimap2 path and the --threads documentation. --- src/bin/coverm.rs | 21 ++++++++++++++++++++- src/mapping_index_maintenance.rs | 5 +++++ 2 files changed, 25 insertions(+), 1 deletion(-) diff --git a/src/bin/coverm.rs b/src/bin/coverm.rs index 9e21e82..2c73da3 100644 --- a/src/bin/coverm.rs +++ b/src/bin/coverm.rs @@ -729,6 +729,26 @@ fn main() { .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. @@ -782,7 +802,6 @@ fn main() { MappingProgram::STROBEALIGN => m.get_one::("strobealign-params"), }; for reference in &references { - check_reference_existence(reference, &mapping_program); let db_path = coverm::mapping_index_maintenance::generate_persistent_index( mapping_program, reference, diff --git a/src/mapping_index_maintenance.rs b/src/mapping_index_maintenance.rs index 6a86a11..febf3a1 100644 --- a/src/mapping_index_maintenance.rs +++ b/src/mapping_index_maintenance.rs @@ -222,6 +222,7 @@ fn execute_index_command(mut cmd: std::process::Command, mapping_program: Mappin 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) @@ -250,6 +251,9 @@ fn generate_strobealign_persistent_index( 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() { @@ -394,6 +398,7 @@ pub fn generate_persistent_index( return generate_strobealign_persistent_index( reference_path, output_directory, + num_threads, index_creation_options, ); } From 04b583e9ff7e9055b4f2fc006ec4cc92f0a01ce1 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 10 Jun 2026 20:47:09 +0000 Subject: [PATCH 5/5] Add bowtie2 as a mapper Support bowtie2 alongside bwa, minimap2 and strobealign, both for read mapping (coverm contig/genome/make) and for pre-generating databases (coverm makedb). - Add a BOWTIE2 variant to MappingProgram and thread it through all the exhaustive matches (parameter selection, dependency checking, index generation, command building). - Build bowtie2 mapping commands separately, since bowtie2 takes the index via -x, reads via -1/-2/-U/--interleaved and threads via -p. Unlike the other mappers, bowtie2 does not auto-detect FASTA vs FASTQ input, so peek at the first read record (via needletail, which also handles gzip/bzip2/xz) and pass -f for FASTA reads. - Generate bowtie2 indexes with bowtie2-build, both as temporary indexes during mapping and as persistent databases via makedb (written as a prefix, like BWA). Detect a pre-built index beside the reference prefix so it can be reused. - Add --bowtie2-params to contig/genome/make/makedb, a check_for_bowtie2 dependency check, CLI help/table entries, README docs and the bowtie2 conda dependency. - Add tests for direct bowtie2 mapping and makedb round-trips. https://claude.ai/code/session_0197qbhqSpTk8JEAPcRDHcvV --- README.md | 12 ++- pixi.lock | 123 ++++++++++++++++++++----------- pixi.toml | 1 + src/bam_generator.rs | 53 ++++++++++++- src/bin/coverm.rs | 16 ++++ src/cli.rs | 52 ++++++++++++- src/external_command_checker.rs | 8 ++ src/mapping_index_maintenance.rs | 54 +++++++++++++- src/mapping_parameters.rs | 1 + tests/test_cmdline.rs | 80 ++++++++++++++++++++ 10 files changed, 350 insertions(+), 50 deletions(-) diff --git a/README.md b/README.md index a357043..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 @@ -129,7 +130,7 @@ 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` and `strobealign` are supported: +`bwa-mem`/`bwa-mem2`, `bowtie2` and `strobealign` are supported: ```bash # Generate a short-read minimap2 database @@ -143,6 +144,15 @@ coverm contig \ # 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 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/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 2c73da3..49b7ce4 100644 --- a/src/bin/coverm.rs +++ b/src/bin/coverm.rs @@ -800,6 +800,7 @@ fn main() { | 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( @@ -846,6 +847,12 @@ fn main() { --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" + ); + } } } } @@ -934,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"); @@ -1021,6 +1033,7 @@ fn mapping_program_from_name(name: Option<&str>) -> 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: {:?}", name), } @@ -1045,6 +1058,9 @@ fn check_mapping_program_dependencies(mapping_program: MappingProgram) { MappingProgram::STROBEALIGN => { external_command_checker::check_for_strobealign(); } + MappingProgram::BOWTIE2 => { + external_command_checker::check_for_bowtie2(); + } } } diff --git a/src/cli.rs b/src/cli.rs index 87f580d..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]", @@ -494,7 +506,8 @@ pub fn makedb_full_help() -> Manual { 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 strobealign databases, the reference FASTA is copied into the \ + '-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\ @@ -550,6 +563,11 @@ pub fn makedb_full_help() -> Manual { ], &[&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!( @@ -579,6 +597,11 @@ pub fn makedb_full_help() -> Manual { 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]", )), ); @@ -1414,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") @@ -1949,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") @@ -2327,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") @@ -2395,6 +2439,12 @@ Ben J. Woodcroft .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( 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 febf3a1..87c2a1e 100644 --- a/src/mapping_index_maintenance.rs +++ b/src/mapping_index_maintenance.rs @@ -91,6 +91,7 @@ fn build_index_command( | 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; @@ -103,6 +104,14 @@ fn build_index_command( .arg(index_path) .arg(reference_path); } + MappingProgram::BOWTIE2 => { + if let Some(t) = num_threads { + cmd.arg("--threads").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 @@ -128,7 +137,8 @@ fn build_index_command( MappingProgram::MINIMAP2_NO_PRESET | MappingProgram::BWA_MEM | MappingProgram::BWA_MEM2 - | MappingProgram::STROBEALIGN => {} + | MappingProgram::STROBEALIGN + | MappingProgram::BOWTIE2 => {} }; if let Some(t) = num_threads { cmd.arg("-t").arg(format!("{t}")); @@ -311,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 @@ -351,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, @@ -378,6 +427,7 @@ pub fn mapping_program_db_name(mapping_program: MappingProgram) -> &'static str MappingProgram::MINIMAP2_LR_HQ => "minimap2-lr-hq", MappingProgram::MINIMAP2_NO_PRESET => "minimap2-no-preset", MappingProgram::STROBEALIGN => "strobealign", + MappingProgram::BOWTIE2 => "bowtie2", } } @@ -419,7 +469,7 @@ pub fn generate_persistent_index( | 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::BWA_MEM | MappingProgram::BWA_MEM2 | MappingProgram::BOWTIE2 => { std::path::Path::new(output_directory) .join(format!("{reference_stem}.{db_program_name}")) } 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/tests/test_cmdline.rs b/tests/test_cmdline.rs index 4c88b74..192f289 100644 --- a/tests/test_cmdline.rs +++ b/tests/test_cmdline.rs @@ -634,6 +634,86 @@ mod tests { .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()