extern crate coverm;
use coverm::bam_generator::*;
use coverm::cli::*;
use coverm::coverage_printer::*;
use coverm::coverage_takers::*;
use coverm::external_command_checker;
use coverm::filter;
use coverm::genome_exclusion::*;
use coverm::genomes_and_contigs::GenomesAndContigs;
use coverm::mapping_index_maintenance::check_reference_existence;
use coverm::mapping_parameters::*;
use coverm::mosdepth_genome_coverage_estimators::*;
use coverm::shard_bam_reader::*;
use coverm::strobealign_aemb::strobealign_aemb_coverage;
use coverm::FlagFilter;
use coverm::OutputWriter;
use coverm::CONCATENATED_FASTA_FILE_SEPARATOR;
extern crate galah;
use galah::ClusterDistanceFinder;
extern crate rust_htslib;
use rust_htslib::bam;
use rust_htslib::bam::Read;
use std::collections::HashSet;
use std::env;
use std::process;
use std::str;
extern crate clap;
use clap::*;
extern crate clap_complete;
use clap_complete::{generate, Shell};
#[macro_use]
extern crate log;
extern crate tempfile;
use tempfile::NamedTempFile;
extern crate bird_tool_utils;
use bird_tool_utils::clap_utils::set_log_level as set_log_level_bird_tool_utils;
const CONCATENATED_REFERENCE_CACHE_STEM: &str = "coverm-genome";
const DEFAULT_MAPPING_SOFTWARE_ENUM: MappingProgram = MappingProgram::MINIMAP2_SR;
fn main() {
let mut app = build_cli();
let matches = app.clone().get_matches();
set_log_level(&matches, false);
let mut print_stream;
match matches.subcommand_name() {
Some("genome") => {
let m = matches.subcommand_matches("genome").unwrap();
bird_tool_utils::clap_utils::print_full_help_if_needed(m, genome_full_help());
set_log_level(m, true);
manually_check_args_at_runtime(m);
print_stream = OutputWriter::generate(m.get_one::<String>("output-file").map(|x| &**x));
let genome_names_content: Vec<u8>;
let mut estimators_and_taker =
EstimatorsAndTaker::generate_from_clap(m, print_stream.clone());
estimators_and_taker = estimators_and_taker.print_headers(
if m.contains_id("gff") {
"Gene\tContig\tGenome"
} else {
"Genome"
},
print_stream.clone(),
);
let filter_params = FilterParameters::generate_from_clap(m);
let separator = parse_separator(m);
let genomes_and_contigs_option_predereplication = if m.get_flag("sharded")
&& !m.contains_id("separator")
&& !m.get_flag("dereplicate")
&& !m.get_flag("single-genome")
{
parse_all_genome_definitions(m)
} else {
None
};
let mut genome_exclusion_filter_separator_type: Option<SeparatorGenomeExclusionFilter> =
None;
let mut genome_exclusion_filter_non_type: Option<NoExclusionGenomeFilter> = None;
let mut genome_exclusion_genomes_and_contigs: Option<GenomesAndContigsExclusionFilter> =
None;
enum GenomeExclusionTypes {
Separator,
None,
GenomesAndContigs,
}
let genome_exclusion_type = {
if m.get_flag("sharded") {
if m.contains_id("exclude-genomes-from-deshard") {
let filename = m.get_one::<String>("exclude-genomes-from-deshard").unwrap();
genome_names_content = std::fs::read(filename).unwrap_or_else(|_| {
panic!(
"Failed to open file '{}' containing list of excluded genomes",
filename
)
});
let mut genome_names_hash: HashSet<&[u8]> = HashSet::new();
for n in genome_names_content.split(|s| *s == b"\n"[0]) {
if n != b"" {
genome_names_hash.insert(n);
}
}
if genome_names_hash.is_empty() {
warn!("No genomes read in that are to be excluded from desharding process");
genome_exclusion_filter_non_type = Some(NoExclusionGenomeFilter {});
GenomeExclusionTypes::None
} else {
info!("Read in {} distinct genomes to exclude from desharding process e.g. '{}'",
genome_names_hash.len(),
std::str::from_utf8(genome_names_hash.iter().next().unwrap())
.unwrap());
if let Some(s) = separator {
genome_exclusion_filter_separator_type =
Some(SeparatorGenomeExclusionFilter {
split_char: s,
excluded_genomes: genome_names_hash,
});
GenomeExclusionTypes::Separator
} else {
match genomes_and_contigs_option_predereplication {
Some(ref gc) => {
genome_exclusion_genomes_and_contigs =
Some(GenomesAndContigsExclusionFilter {
genomes_and_contigs: gc,
excluded_genomes: genome_names_hash,
});
GenomeExclusionTypes::GenomesAndContigs
}
None => unreachable!(),
}
}
}
} else {
debug!("Not excluding any genomes during the deshard process");
genome_exclusion_filter_non_type = Some(NoExclusionGenomeFilter {});
GenomeExclusionTypes::None
}
} else {
genome_exclusion_filter_non_type = Some(NoExclusionGenomeFilter {});
GenomeExclusionTypes::None
}
};
if m.contains_id("bam-files") {
let bam_files: Vec<&str> = m
.get_many::<String>("bam-files")
.unwrap()
.map(|s| &**s)
.collect();
let genomes_and_contigs_option = parse_all_genome_definitions(m);
if filter_params.doing_filtering() {
run_genome(
coverm::bam_generator::generate_filtered_bam_readers_from_bam_files(
bam_files,
filter_params.flag_filters,
filter_params.min_aligned_length_single,
filter_params.min_percent_identity_single,
filter_params.min_aligned_percent_single,
filter_params.min_mapq,
filter_params.min_aligned_length_pair,
filter_params.min_percent_identity_pair,
filter_params.min_aligned_percent_pair,
),
m,
&mut estimators_and_taker,
separator,
&genomes_and_contigs_option,
&mut print_stream,
);
} else if m.get_flag("sharded") {
external_command_checker::check_for_samtools();
let sort_threads = *m.get_one::<u16>("threads").unwrap();
match genome_exclusion_type {
GenomeExclusionTypes::None => {
run_genome(
coverm::shard_bam_reader::generate_sharded_bam_reader_from_bam_files(
bam_files,
sort_threads,
&genome_exclusion_filter_non_type.unwrap()),
m,
&mut estimators_and_taker,
separator,
&genomes_and_contigs_option,
&mut print_stream,);
}
GenomeExclusionTypes::Separator => {
run_genome(
coverm::shard_bam_reader::generate_sharded_bam_reader_from_bam_files(
bam_files,
sort_threads,
&genome_exclusion_filter_separator_type.unwrap()),
m,
&mut estimators_and_taker,
separator,
&genomes_and_contigs_option,
&mut print_stream,);
}
GenomeExclusionTypes::GenomesAndContigs => {
run_genome(
coverm::shard_bam_reader::generate_sharded_bam_reader_from_bam_files(
bam_files,
sort_threads,
&genome_exclusion_genomes_and_contigs.unwrap()),
m,
&mut estimators_and_taker,
separator,
&genomes_and_contigs_option,
&mut print_stream,);
}
}
} else {
run_genome(
coverm::bam_generator::generate_named_bam_readers_from_bam_files(bam_files),
m,
&mut estimators_and_taker,
separator,
&genomes_and_contigs_option,
&mut print_stream,
);
}
} else {
let mapping_program = parse_mapping_program(m);
external_command_checker::check_for_samtools();
let genome_fasta_files_opt = resolve_and_checkm_filter_genomes(
m,
"All genomes were filtered out, so none remain to be mapped to",
);
let (concatenated_genomes, genomes_and_contigs_option) =
match m.contains_id("reference") {
true => {
check_reference_existence(
m.get_one::<String>("reference").unwrap(),
&mapping_program,
);
match genome_fasta_files_opt {
Some(genome_paths) => (
None,
extract_genomes_and_contigs_option(
m,
&genome_paths.iter().map(|s| s.as_str()).collect(),
),
),
None => match m.get_one::<String>("genome-definition") {
Some(definition_path) => (
None,
Some(coverm::genome_parsing::read_genome_definition_file(
definition_path,
)),
),
None => (None, None),
},
}
}
false => {
let dereplicated_genomes: Vec<String> = if m.get_flag("dereplicate") {
dereplicate(m, &genome_fasta_files_opt.unwrap())
} else {
genome_fasta_files_opt.unwrap()
};
info!("Profiling {} genomes", dereplicated_genomes.len());
debug!(
"Dereplicated genomes: {:?}",
dereplicated_genomes
.iter()
.map(|s| s.as_str())
.collect::<Vec<_>>()
);
let list_of_genome_fasta_files = &dereplicated_genomes;
info!(
"Generating concatenated reference FASTA file of {} genomes ..",
list_of_genome_fasta_files.len()
);
(
Some(
coverm::mapping_index_maintenance::generate_concatenated_fasta_file(
list_of_genome_fasta_files,
),
),
None,
)
}
};
if filter_params.doing_filtering() {
debug!("Mapping and filtering..");
let generator_sets = get_streamed_filtered_bam_readers(
m,
mapping_program,
&concatenated_genomes,
&filter_params,
);
let mut all_generators = vec![];
let mut indices = vec![]; for set in generator_sets {
indices.push(set.index);
for g in set.generators {
all_generators.push(g)
}
}
run_genome(
all_generators,
m,
&mut estimators_and_taker,
separator,
&genomes_and_contigs_option,
&mut print_stream,
);
} else if m.get_flag("sharded") {
match genome_exclusion_type {
GenomeExclusionTypes::None => {
run_genome(
get_sharded_bam_readers(
m,
mapping_program,
&concatenated_genomes,
&genome_exclusion_filter_non_type.unwrap(),
),
m,
&mut estimators_and_taker,
separator,
&genomes_and_contigs_option,
&mut print_stream,
);
}
GenomeExclusionTypes::Separator => {
run_genome(
get_sharded_bam_readers(
m,
mapping_program,
&concatenated_genomes,
&genome_exclusion_filter_separator_type.unwrap(),
),
m,
&mut estimators_and_taker,
separator,
&genomes_and_contigs_option,
&mut print_stream,
);
}
GenomeExclusionTypes::GenomesAndContigs => {
run_genome(
get_sharded_bam_readers(
m,
mapping_program,
&concatenated_genomes,
&genome_exclusion_genomes_and_contigs.unwrap(),
),
m,
&mut estimators_and_taker,
separator,
&genomes_and_contigs_option,
&mut print_stream,
);
}
}
} else {
let generator_sets =
get_streamed_bam_readers(m, mapping_program, &concatenated_genomes);
let mut all_generators = vec![];
let mut indices = vec![]; for set in generator_sets {
indices.push(set.index);
for g in set.generators {
all_generators.push(g)
}
}
run_genome(
all_generators,
m,
&mut estimators_and_taker,
separator,
&genomes_and_contigs_option,
&mut print_stream,
);
};
}
}
Some("filter") => {
let m = matches.subcommand_matches("filter").unwrap();
bird_tool_utils::clap_utils::print_full_help_if_needed(m, filter_full_help());
set_log_level(m, true);
manually_check_args_at_runtime(m);
let bam_files: Vec<&str> = m
.get_many::<String>("bam-files")
.unwrap()
.map(|x| &**x)
.collect();
let output_bam_files: Vec<&str> = m
.get_many::<String>("output-bam-files")
.unwrap()
.map(|x| &**x)
.collect();
if bam_files.len() != output_bam_files.len() {
error!("The number of input BAM files must be the same as the number output");
process::exit(1);
}
let filter_params = FilterParameters::generate_from_clap(m);
let num_threads: u16 = *m.get_one::<u16>("threads").unwrap();
for (bam, output) in bam_files.iter().zip(output_bam_files.iter()) {
let reader = bam::Reader::from_path(bam)
.unwrap_or_else(|_| panic!("Unable to find BAM file {}", bam));
let header = bam::header::Header::from_template(reader.header());
let mut writer =
bam::Writer::from_path(output, &header, rust_htslib::bam::Format::Bam)
.unwrap_or_else(|_| panic!("Failed to write BAM file {}", output));
writer
.set_threads(num_threads as usize)
.expect("Failed to set num threads in writer");
let mut filtered = filter::ReferenceSortedBamFilter::new(
reader,
filter_params.flag_filters.clone(),
filter_params.min_aligned_length_single,
filter_params.min_percent_identity_single,
filter_params.min_aligned_percent_single,
filter_params.min_mapq,
filter_params.min_aligned_length_pair,
filter_params.min_percent_identity_pair,
filter_params.min_aligned_percent_pair,
!m.get_flag("inverse"),
);
let mut record = bam::record::Record::new();
loop {
match filtered.read(&mut record) {
None => {
break;
}
Some(Ok(())) => {}
Some(e) => {
panic!("Failure to read filtered BAM record: {:?}", e)
}
}
debug!("Writing.. {:?}", record.qname());
writer.write(&record).expect("Failed to write BAM record");
}
}
}
Some("contig") => {
let m = matches.subcommand_matches("contig").unwrap();
bird_tool_utils::clap_utils::print_full_help_if_needed(m, contig_full_help());
set_log_level(m, true);
manually_check_args_at_runtime(m);
let print_zeros = !m.get_flag("no-zeros");
let mut filter_params1 = FilterParameters::generate_from_clap(m);
filter_params1.add_metabat_filtering_if_required(m);
let filter_params = filter_params1;
let threads = *m.get_one::<u16>("threads").unwrap();
print_stream = OutputWriter::generate(m.get_one::<String>("output-file").map(|x| &**x));
let gene_definitions: Option<coverm::genes::GeneDefinitions> =
m.get_one::<String>("gff").map(|gff_path| {
if doing_metabat(m) {
error!("The metabat method cannot be used with --gff");
process::exit(1);
}
let methods: Vec<&str> = m
.get_many::<String>("methods")
.unwrap()
.map(|x| x.as_str())
.collect();
if methods.contains(&"strobealign-aemb") {
error!("The strobealign-aemb method cannot be used with --gff");
process::exit(1);
}
coverm::genes::GeneDefinitions::read_gff(
gff_path,
m.get_one::<String>("gff-feature-type").map(|x| x.as_str()),
)
});
let mut estimators_and_taker =
EstimatorsAndTaker::generate_from_clap(m, print_stream.clone());
estimators_and_taker = estimators_and_taker.print_headers(
if gene_definitions.is_some() {
"Gene\tContig"
} else {
"Contig"
},
print_stream.clone(),
);
if let CoverageEstimator::StrobealignAembEstimator {} =
estimators_and_taker.estimators[0]
{
let mapping_params =
MappingParameters::generate_from_clap(m, MappingProgram::STROBEALIGN, &None);
debug!(
"Running strobealign-aemb coverage with mapping parameters: {mapping_params:?}"
);
strobealign_aemb_coverage(
mapping_params,
&mut estimators_and_taker.taker,
threads,
&mut estimators_and_taker.printer,
&mut print_stream,
);
} else if m.contains_id("bam-files") {
let bam_files: Vec<&str> = m
.get_many::<String>("bam-files")
.unwrap()
.map(|x| &**x)
.collect();
if filter_params.doing_filtering() {
let bam_readers =
coverm::bam_generator::generate_filtered_bam_readers_from_bam_files(
bam_files,
filter_params.flag_filters.clone(),
filter_params.min_aligned_length_single,
filter_params.min_percent_identity_single,
filter_params.min_aligned_percent_single,
filter_params.min_mapq,
filter_params.min_aligned_length_pair,
filter_params.min_percent_identity_pair,
filter_params.min_aligned_percent_pair,
);
run_contig(
&mut estimators_and_taker,
bam_readers,
print_zeros,
filter_params.flag_filters,
threads,
&gene_definitions,
&mut print_stream,
);
} else if m.get_flag("sharded") {
external_command_checker::check_for_samtools();
let bam_readers =
coverm::shard_bam_reader::generate_sharded_bam_reader_from_bam_files(
bam_files,
threads,
&NoExclusionGenomeFilter {},
);
run_contig(
&mut estimators_and_taker,
bam_readers,
print_zeros,
filter_params.flag_filters,
threads,
&gene_definitions,
&mut print_stream,
);
} else {
let bam_readers =
coverm::bam_generator::generate_named_bam_readers_from_bam_files(bam_files);
run_contig(
&mut estimators_and_taker,
bam_readers,
print_zeros,
filter_params.flag_filters,
threads,
&gene_definitions,
&mut print_stream,
);
}
} else {
let mapping_program = parse_mapping_program(m);
external_command_checker::check_for_samtools();
if filter_params.doing_filtering() {
debug!("Filtering..");
let generator_sets = get_streamed_filtered_bam_readers(
m,
mapping_program,
&None,
&filter_params,
);
let mut all_generators = vec![];
let mut indices = vec![]; for set in generator_sets {
indices.push(set.index);
for g in set.generators {
all_generators.push(g)
}
}
debug!("Finished collecting generators.");
run_contig(
&mut estimators_and_taker,
all_generators,
print_zeros,
filter_params.flag_filters,
threads,
&gene_definitions,
&mut print_stream,
);
} else if m.get_flag("sharded") {
let generator_sets = get_sharded_bam_readers(
m,
mapping_program,
&None,
&NoExclusionGenomeFilter {},
);
run_contig(
&mut estimators_and_taker,
generator_sets,
print_zeros,
filter_params.flag_filters,
threads,
&gene_definitions,
&mut print_stream,
);
} else {
debug!("Not filtering..");
let generator_sets = get_streamed_bam_readers(m, mapping_program, &None);
let mut all_generators = vec![];
let mut indices = vec![]; for set in generator_sets {
indices.push(set.index);
for g in set.generators {
all_generators.push(g)
}
}
run_contig(
&mut estimators_and_taker,
all_generators,
print_zeros,
filter_params.flag_filters,
threads,
&gene_definitions,
&mut print_stream,
);
}
}
}
Some("make") => {
let m = matches.subcommand_matches("make").unwrap();
bird_tool_utils::clap_utils::print_full_help_if_needed(m, make_full_help());
set_log_level(m, true);
manually_check_args_at_runtime(m);
let mapping_program = parse_mapping_program(m);
external_command_checker::check_for_samtools();
let output_directory = m.get_one::<String>("output-directory").unwrap();
setup_bam_cache_directory(output_directory);
let params = MappingParameters::generate_from_clap(m, mapping_program, &None);
let mut generator_sets = vec![];
let discard_unmapped_reads = m.get_flag("discard-unmapped");
for reference_wise_params in params {
let mut bam_readers = vec![];
let index = setup_mapping_index(&reference_wise_params, m, mapping_program);
let mut unique_names = HashSet::new();
for p in reference_wise_params {
let name =
generate_cached_bam_file_name(output_directory, p.reference, p.read1);
bam_readers.push(
coverm::bam_generator::generate_bam_maker_generator_from_reads(
mapping_program,
index.as_ref(),
p.read1,
p.read2,
p.read_format.clone(),
p.threads,
&name.clone(),
discard_unmapped_reads,
p.mapping_options,
),
);
if !unique_names.insert(name.clone()) {
error!("Duplicate output file name: {name}");
std::process::exit(1);
}
}
debug!("Finished BAM setup");
let to_return = BamGeneratorSet {
generators: bam_readers,
index,
};
generator_sets.push(to_return);
}
let mut i = 1;
for generator_set in generator_sets {
for generator in generator_set.generators {
info!("Running mapping number {i} ..");
generator.start().finish();
i += 1;
}
}
}
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::<String>("output-directory").unwrap();
setup_bam_cache_directory(output_directory);
let concatenated_reference: Option<String> = match resolve_and_checkm_filter_genomes(
m,
"All genomes were filtered out, so none remain to make a database from",
) {
Some(genome_fasta_files) => {
let genome_fasta_files = if m.get_flag("dereplicate") {
dereplicate(m, &genome_fasta_files)
} else {
genome_fasta_files
};
info!(
"Generating concatenated reference FASTA file of {} genomes ..",
genome_fasta_files.len()
);
let concatenated_path = std::path::Path::new(output_directory)
.join("coverm_concatenated_genomes.fna");
coverm::mapping_index_maintenance::generate_concatenated_fasta_file_at(
&genome_fasta_files,
&concatenated_path,
);
info!(
"Wrote concatenated reference to {}",
concatenated_path.display()
);
Some(concatenated_path.to_string_lossy().to_string())
}
None => None,
};
let built_from_genomes = concatenated_reference.is_some();
let references: Vec<String> = match concatenated_reference {
Some(path) => vec![path],
None => {
let references: Vec<String> = m
.get_many::<String>("reference")
.expect("No reference provided")
.map(|s| s.to_string())
.collect();
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);
}
}
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);
}
}
}
references
}
};
let mut mapping_programs = vec![];
let mut seen_mappers = HashSet::new();
for mapper_name in m.get_many::<String>("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::<u16>("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::<String>("bwa-params")
}
MappingProgram::MINIBWA
| MappingProgram::RAMMAP_SR
| MappingProgram::RAMMAP_ONT
| MappingProgram::RAMMAP_PB
| MappingProgram::RAMMAP_HIFI
| MappingProgram::RAMMAP_LR_HQ
| MappingProgram::RAMMAP_NO_PRESET => None,
MappingProgram::MINIMAP2_SR
| MappingProgram::MINIMAP2_ONT
| MappingProgram::MINIMAP2_PB
| MappingProgram::MINIMAP2_HIFI
| MappingProgram::MINIMAP2_LR_HQ
| MappingProgram::MINIMAP2_NO_PRESET => m.get_one::<String>("minimap2-params"),
MappingProgram::STROBEALIGN => m.get_one::<String>("strobealign-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 {
use coverm::mapping_index_maintenance::{makedb_usage_command, MakedbUsageMode};
let db_name =
coverm::mapping_index_maintenance::mapping_program_db_name(*mapping_program);
info!(
"To use the {db_name} database, run e.g.:\n {}",
makedb_usage_command(*mapping_program, db_path, MakedbUsageMode::Contig)
);
if built_from_genomes {
info!(
"To recover per-genome coverage from the {db_name} database, run e.g.:\n {}",
makedb_usage_command(*mapping_program, db_path, MakedbUsageMode::Genome)
);
}
}
}
Some("shell-completion") => {
let m = matches.subcommand_matches("shell-completion").unwrap();
set_log_level(m, true);
manually_check_args_at_runtime(m);
let mut file = std::fs::File::create(m.get_one::<String>("output-file").unwrap())
.expect("failed to open output file");
if let Some(generator) = m.get_one::<Shell>("shell").copied() {
let mut cmd = build_cli();
info!("Generating completion script for shell {generator}");
let name = cmd.get_name().to_string();
generate(generator, &mut cmd, name, &mut file);
}
}
Some("cluster") => {
galah::cluster_argument_parsing::run_cluster_subcommand(
&matches,
"CoverM",
crate_version!(),
);
}
_ => {
app.print_help().unwrap();
println!();
}
}
}
fn manually_check_args_at_runtime(m: &clap::ArgMatches) {
let run_checkm2 = m
.try_get_one::<bool>("run-checkm2")
.ok()
.flatten()
.copied()
.unwrap_or(false);
if ((m.try_get_one::<f32>("min-completeness").is_ok() && m.contains_id("min-completeness"))
|| (m.try_get_one::<f32>("max-contamination").is_ok()
&& m.contains_id("max-contamination")))
&& !m.contains_id("checkm-tab-table")
&& !m.contains_id("checkm2-quality-report")
&& !m.contains_id("genome-info")
&& !run_checkm2
{
error!("You must provide a CheckM tab table, CheckM2 quality report, genome info file, or use --run-checkm2 to use --min-completeness or --max-contamination");
}
}
fn setup_mapping_index(
reference_wise_params: &SingleReferenceMappingParameters,
m: &clap::ArgMatches,
mapping_program: MappingProgram,
) -> Box<dyn coverm::mapping_index_maintenance::MappingIndex> {
match mapping_program {
MappingProgram::BWA_MEM | MappingProgram::BWA_MEM2 => {
coverm::mapping_index_maintenance::generate_bwa_index(
reference_wise_params.reference,
None,
mapping_program,
)
}
MappingProgram::MINIBWA => coverm::mapping_index_maintenance::generate_minibwa_index(
reference_wise_params.reference,
Some(*m.get_one::<u16>("threads").unwrap()),
None,
),
MappingProgram::MINIMAP2_SR
| MappingProgram::MINIMAP2_ONT
| MappingProgram::MINIMAP2_HIFI
| MappingProgram::MINIMAP2_PB
| MappingProgram::MINIMAP2_LR_HQ
| MappingProgram::MINIMAP2_NO_PRESET => {
if m.get_flag("minimap2-reference-is-index") || reference_wise_params.len() == 1 {
info!("Not pre-generating minimap2 index");
if m.get_flag("minimap2-reference-is-index") {
warn!(
"Minimap2 uses mapping parameters defined when the index was created, \
not parameters defined when mapping. Proceeding on the assumption that you \
passed the correct parameters when creating the minimap2 index."
);
}
Box::new(coverm::mapping_index_maintenance::VanillaIndexStruct::new(
reference_wise_params.reference,
))
} else {
coverm::mapping_index_maintenance::generate_minimap2_index(
reference_wise_params.reference,
Some(*m.get_one::<u16>("threads").unwrap()),
Some(
m.get_one::<String>("minimap2-params")
.unwrap_or(&"".to_string()),
),
mapping_program,
)
}
}
MappingProgram::RAMMAP_SR
| MappingProgram::RAMMAP_ONT
| MappingProgram::RAMMAP_PB
| MappingProgram::RAMMAP_HIFI
| MappingProgram::RAMMAP_LR_HQ
| MappingProgram::RAMMAP_NO_PRESET => {
info!("Not pre-generating rammap index");
Box::new(coverm::mapping_index_maintenance::VanillaIndexStruct::new(
reference_wise_params.reference,
))
}
MappingProgram::STROBEALIGN => {
info!("Not pre-generating strobealign index");
if m.get_flag("strobealign-use-index") {
warn!(
"Strobealign uses mapping parameters defined when the index was created, \
not parameters defined when mapping. Proceeding on the assumption that you \
passed the correct parameters when creating the strobealign index."
);
Box::new(
coverm::mapping_index_maintenance::PregeneratedStrobealignIndexStruct::new(
reference_wise_params.reference,
),
)
} else {
Box::new(coverm::mapping_index_maintenance::VanillaIndexStruct::new(
reference_wise_params.reference,
))
}
}
}
}
fn dereplicate(m: &clap::ArgMatches, genome_fasta_files: &[String]) -> Vec<String> {
info!(
"Found {} genomes specified before dereplication",
genome_fasta_files.len()
);
let reference_genomes_owned: Option<Vec<String>> = if let Some(refs) = m.get_many::<String>(
&coverm::cli::COVERM_CLUSTER_COMMAND_DEFINITION.dereplication_reference_genomes_argument,
) {
Some(refs.cloned().collect())
} else if let Some(ref_file) = m.get_one::<String>(
&coverm::cli::COVERM_CLUSTER_COMMAND_DEFINITION
.dereplication_reference_genomes_list_argument,
) {
let content = std::fs::read_to_string(ref_file)
.unwrap_or_else(|_| panic!("Failed to read reference genomes list file: {}", ref_file));
Some(
content
.lines()
.filter(|line| !line.trim().is_empty())
.map(|s| s.split('\t').next().unwrap_or("").to_string())
.collect(),
)
} else {
None
};
let reference_genomes: Option<Vec<&str>> = reference_genomes_owned
.as_ref()
.map(|refs| refs.iter().map(String::as_str).collect());
let (combined_genomes, ref_genomes_for_clusterer): (Vec<String>, Option<Vec<&str>>) =
if let Some(ref_genomes) = &reference_genomes {
let mut combined = ref_genomes
.iter()
.map(|s| s.to_string())
.collect::<Vec<String>>();
combined.extend(genome_fasta_files.iter().cloned());
(combined, Some(ref_genomes.clone()))
} else {
(genome_fasta_files.to_vec(), None)
};
let clusterer = galah::cluster_argument_parsing::generate_galah_clusterer(
&combined_genomes,
&None, false, m,
&coverm::cli::COVERM_CLUSTER_COMMAND_DEFINITION,
ref_genomes_for_clusterer.as_deref(),
None, )
.expect("Failed to parse galah clustering arguments correctly");
let cluster_outputs = galah::cluster_argument_parsing::setup_galah_outputs(
m,
&coverm::cli::COVERM_CLUSTER_COMMAND_DEFINITION,
);
info!(
"Dereplicating genomes at {}% ANI ..",
clusterer.clusterer.get_ani_threshold()
);
let cluster_indices = clusterer.cluster();
info!(
"Finished dereplication, finding {} representative genomes.",
cluster_indices.len()
);
debug!("Found cluster indices: {cluster_indices:?}");
let reps = cluster_indices
.iter()
.map(|cluster| combined_genomes[cluster[0]].clone())
.collect::<Vec<_>>();
debug!("Found cluster representatives: {reps:?}");
galah::cluster_argument_parsing::write_galah_outputs(
cluster_outputs,
&cluster_indices,
&clusterer.genome_fasta_paths,
None, );
reps
}
fn resolve_and_checkm_filter_genomes(
m: &clap::ArgMatches,
no_survivors_error: &str,
) -> Option<Vec<String>> {
match bird_tool_utils::clap_utils::parse_list_of_genome_fasta_files(m, false) {
Ok(paths) => {
if paths.is_empty() {
error!("Genome paths were described, but ultimately none were found");
process::exit(1);
}
if m.contains_id("checkm-tab-table")
|| m.contains_id("checkm2-quality-report")
|| m.contains_id("genome-info")
|| m.get_flag(
&coverm::cli::COVERM_CLUSTER_COMMAND_DEFINITION
.dereplication_run_checkm2_argument,
)
{
let genomes_after_filtering =
galah::cluster_argument_parsing::filter_genomes_through_checkm(
&paths,
m,
&coverm::cli::COVERM_CLUSTER_COMMAND_DEFINITION,
None, )
.expect("Error parsing CheckM-related options");
info!(
"After filtering by CheckM, {} genomes remained",
genomes_after_filtering.len()
);
if genomes_after_filtering.is_empty() {
error!("{no_survivors_error}");
process::exit(1);
}
Some(
genomes_after_filtering
.iter()
.map(|s| s.to_string())
.collect(),
)
} else {
Some(paths)
}
}
Err(_) => None,
}
}
fn parse_mapping_program(m: &clap::ArgMatches) -> MappingProgram {
let mapping_program = mapping_program_from_name(m.get_one::<String>("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,
Some("minimap2-ont") => MappingProgram::MINIMAP2_ONT,
Some("minimap2-pb") => MappingProgram::MINIMAP2_PB,
Some("minimap2-hifi") => MappingProgram::MINIMAP2_HIFI,
Some("minimap2-lr-hq") => MappingProgram::MINIMAP2_LR_HQ,
Some("minimap2-no-preset") => MappingProgram::MINIMAP2_NO_PRESET,
Some("strobealign") => MappingProgram::STROBEALIGN,
Some("minibwa") => MappingProgram::MINIBWA,
Some("rammap-sr") => MappingProgram::RAMMAP_SR,
Some("rammap-ont") => MappingProgram::RAMMAP_ONT,
Some("rammap-pb") => MappingProgram::RAMMAP_PB,
Some("rammap-hifi") => MappingProgram::RAMMAP_HIFI,
Some("rammap-lr-hq") => MappingProgram::RAMMAP_LR_HQ,
Some("rammap-no-preset") => MappingProgram::RAMMAP_NO_PRESET,
None => DEFAULT_MAPPING_SOFTWARE_ENUM,
_ => 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();
}
MappingProgram::BWA_MEM2 => {
external_command_checker::check_for_bwa_mem2();
}
MappingProgram::MINIMAP2_SR
| MappingProgram::MINIMAP2_ONT
| MappingProgram::MINIMAP2_HIFI
| MappingProgram::MINIMAP2_PB
| MappingProgram::MINIMAP2_LR_HQ
| MappingProgram::MINIMAP2_NO_PRESET => {
external_command_checker::check_for_minimap2();
}
MappingProgram::STROBEALIGN => {
external_command_checker::check_for_strobealign();
}
MappingProgram::MINIBWA => {
external_command_checker::check_for_minibwa();
}
MappingProgram::RAMMAP_SR
| MappingProgram::RAMMAP_ONT
| MappingProgram::RAMMAP_PB
| MappingProgram::RAMMAP_HIFI
| MappingProgram::RAMMAP_LR_HQ
| MappingProgram::RAMMAP_NO_PRESET => {
external_command_checker::check_for_rammap();
}
}
}
struct EstimatorsAndTaker {
estimators: Vec<CoverageEstimator>,
taker: CoverageTakerType,
columns_to_normalise: Vec<usize>,
rpkm_column: Option<usize>,
tpm_column: Option<usize>,
printer: CoveragePrinter,
}
fn extract_genomes_and_contigs_option(
m: &clap::ArgMatches,
genome_fasta_files: &Vec<&str>,
) -> Option<GenomesAndContigs> {
match m.contains_id("genome-definition") {
true => Some(coverm::genome_parsing::read_genome_definition_file(
m.get_one::<String>("genome-definition").unwrap(),
)),
false => Some(coverm::genome_parsing::read_genome_fasta_files(
genome_fasta_files,
m.get_flag("use-full-contig-names"),
)),
}
}
fn parse_all_genome_definitions(m: &clap::ArgMatches) -> Option<GenomesAndContigs> {
if m.get_flag("single-genome") || m.contains_id("separator") {
None
} else if m.contains_id("genome-definition") {
Some(coverm::genome_parsing::read_genome_definition_file(
m.get_one::<String>("genome-definition").unwrap(),
))
} else {
extract_genomes_and_contigs_option(
m,
&bird_tool_utils::clap_utils::parse_list_of_genome_fasta_files(m, true)
.expect("Failed to parse genome paths")
.iter()
.map(|s| s.as_str())
.collect(),
)
}
}
fn parse_percentage(m: &clap::ArgMatches, parameter: &str) -> f32 {
if m.contains_id(parameter) {
let mut percentage: f32 = *m.get_one::<f32>(parameter).unwrap_or(&0.0);
if (1.0..=100.0).contains(&percentage) {
percentage /= 100.0;
} else if !(0.0..=100.0).contains(&percentage) {
error!("Invalid alignment percentage: '{percentage}'");
process::exit(1);
}
if m.value_source(parameter) == Some(clap::parser::ValueSource::CommandLine) {
info!("Using {} {}%", parameter, percentage * 100.0);
}
percentage
} else {
0.0
}
}
impl EstimatorsAndTaker {
pub fn generate_from_clap(m: &clap::ArgMatches, stream: OutputWriter) -> EstimatorsAndTaker {
let mut estimators = vec![];
let min_fraction_covered = parse_percentage(m, "min-covered-fraction");
let contig_end_exclusion = *m.get_one::<u64>("contig-end-exclusion").unwrap();
let methods: Vec<&str> = m
.get_many::<String>("methods")
.unwrap()
.map(|x| &**x)
.collect();
let mut columns_to_normalise: Vec<usize> = vec![];
let taker;
let output_format = m.get_one::<String>("output-format").unwrap().as_str();
let printer;
let mut rpkm_column = None;
let mut tpm_column = None;
if doing_metabat(m) {
estimators.push(CoverageEstimator::new_estimator_length());
estimators.push(CoverageEstimator::new_estimator_mean(
min_fraction_covered,
contig_end_exclusion,
false,
));
estimators.push(CoverageEstimator::new_estimator_variance(
min_fraction_covered,
contig_end_exclusion,
));
debug!("Cached regular coverage taker for metabat mode being used");
taker = CoverageTakerType::new_cached_single_float_coverage_taker(estimators.len());
printer = CoveragePrinter::MetabatAdjustedCoveragePrinter;
} else {
for (i, method) in methods.iter().enumerate() {
match *method {
"mean" => {
estimators.push(CoverageEstimator::new_estimator_mean(
min_fraction_covered,
contig_end_exclusion,
false,
)); }
"coverage_histogram" => {
estimators.push(CoverageEstimator::new_estimator_pileup_counts(
min_fraction_covered,
contig_end_exclusion,
));
}
"trimmed_mean" => {
let min = parse_percentage(m, "trim-min");
let max = parse_percentage(m, "trim-max");
estimators.push(CoverageEstimator::new_estimator_trimmed_mean(
min,
max,
min_fraction_covered,
contig_end_exclusion,
));
}
"covered_fraction" => {
estimators.push(CoverageEstimator::new_estimator_covered_fraction(
min_fraction_covered,
));
}
"covered_bases" => {
estimators.push(CoverageEstimator::new_estimator_covered_bases(
min_fraction_covered,
));
}
"rpkm" => {
if rpkm_column.is_some() {
error!("The RPKM column cannot be specified more than once");
process::exit(1);
}
rpkm_column = Some(i);
estimators.push(CoverageEstimator::new_estimator_rpkm(min_fraction_covered))
}
"tpm" => {
if tpm_column.is_some() {
error!("The TPM column cannot be specified more than once");
process::exit(1);
}
tpm_column = Some(i);
estimators.push(CoverageEstimator::new_estimator_tpm(min_fraction_covered))
}
"variance" => {
estimators.push(CoverageEstimator::new_estimator_variance(
min_fraction_covered,
contig_end_exclusion,
));
}
"length" => {
estimators.push(CoverageEstimator::new_estimator_length());
}
"relative_abundance" => {
columns_to_normalise.push(i);
estimators.push(CoverageEstimator::new_estimator_mean(
min_fraction_covered,
contig_end_exclusion,
false,
));
}
"count" => {
estimators.push(CoverageEstimator::new_estimator_read_count());
}
"reads_per_base" => {
estimators.push(CoverageEstimator::new_estimator_reads_per_base());
}
"anir" => {
estimators.push(CoverageEstimator::new_estimator_anir());
}
"strobealign-aemb" => {
if methods.len() > 1 {
error!("Cannot (currently) specify the strobealign-aemb method with any other coverage methods");
process::exit(1);
}
estimators.push(CoverageEstimator::new_estimator_strobealign_aemb());
}
_ => unreachable!(),
};
}
if methods.contains(&"coverage_histogram") {
if methods.len() > 1 {
error!("Cannot specify the coverage_histogram method with any other coverage methods");
process::exit(1);
} else {
debug!("Coverage histogram type coverage taker being used");
taker = CoverageTakerType::new_pileup_coverage_coverage_printer(stream);
printer = CoveragePrinter::StreamedCoveragePrinter;
}
} else if columns_to_normalise.is_empty()
&& rpkm_column.is_none()
&& tpm_column.is_none()
&& output_format == "sparse"
{
debug!("Streaming regular coverage output");
taker =
CoverageTakerType::new_single_float_coverage_streaming_coverage_printer(stream);
printer = CoveragePrinter::StreamedCoveragePrinter;
} else {
debug!(
"Cached regular coverage taker with columns to normlise: {columns_to_normalise:?} and rpkm_column: {rpkm_column:?} and tpm_column: {tpm_column:?}"
);
taker = CoverageTakerType::new_cached_single_float_coverage_taker(estimators.len());
printer = match output_format {
"sparse" => CoveragePrinter::SparseCachedCoveragePrinter,
"dense" => CoveragePrinter::DenseCachedCoveragePrinter {
entry_type: None,
estimator_headers: None,
},
_ => unreachable!(),
}
}
}
if min_fraction_covered != 0.0 {
let die = |estimator_name| {
error!(
"The '{estimator_name}' coverage estimator cannot be used when \
--min-covered-fraction is > 0 as it does not calculate \
the covered fraction. You may wish to set the \
--min-covered-fraction to 0 and/or run this estimator \
separately."
);
process::exit(1)
};
for e in &estimators {
match e {
CoverageEstimator::ReadCountCalculator { .. } => die("counts"),
CoverageEstimator::ReferenceLengthCalculator { .. } => die("length"),
CoverageEstimator::ReadsPerBaseCalculator { .. } => die("reads_per_base"),
CoverageEstimator::AverageIdentityEstimator { .. } => die("anir"),
CoverageEstimator::StrobealignAembEstimator { .. } => die("strobealign-aemb"),
_ => {}
}
}
}
EstimatorsAndTaker {
estimators,
taker,
columns_to_normalise,
rpkm_column,
tpm_column,
printer,
}
}
pub fn print_headers(mut self, entry_type: &str, print_stream: OutputWriter) -> Self {
let mut headers: Vec<String> = vec![];
for e in self.estimators.iter() {
for h in e.column_headers() {
headers.push(h.to_string())
}
}
for i in self.columns_to_normalise.iter() {
headers[*i] = "Relative Abundance (%)".to_string();
}
self.printer
.print_headers(entry_type, headers, print_stream);
self
}
}
fn parse_separator(m: &clap::ArgMatches) -> Option<u8> {
let single_genome = m.get_flag("single-genome");
if single_genome {
Some("0".as_bytes()[0])
} else if m.contains_id("separator") {
m.get_one::<char>("separator").map(|c| *c as u8)
} else if m.contains_id("bam-files") || m.contains_id("reference") {
None
} else {
Some(CONCATENATED_FASTA_FILE_SEPARATOR.as_bytes()[0])
}
}
fn run_genome<
R: coverm::bam_generator::NamedBamReader,
T: coverm::bam_generator::NamedBamReaderGenerator<R>,
>(
bam_generators: Vec<T>,
m: &clap::ArgMatches,
estimators_and_taker: &mut EstimatorsAndTaker,
separator: Option<u8>,
genomes_and_contigs_option: &Option<GenomesAndContigs>,
print_stream: &mut OutputWriter,
) {
let print_zeros = !m.get_flag("no-zeros");
let flag_filter = FilterParameters::generate_from_clap(m).flag_filters;
let single_genome = m.get_flag("single-genome");
let threads = *m.get_one::<u16>("threads").unwrap();
let reads_mapped = if let Some(gff_path) = m.get_one::<String>("gff") {
let gene_definitions = coverm::genes::GeneDefinitions::read_gff(
gff_path,
m.get_one::<String>("gff-feature-type").map(|x| x.as_str()),
);
let genome_namer: Box<coverm::genes::GenomeNamer> = if single_genome {
Box::new(|_: &str| Some("genome1".to_string()))
} else if let Some(sep) = separator {
Box::new(move |contig: &str| {
contig
.split_once(sep as char)
.map(|(genome, _)| genome.to_string())
})
} else if let Some(gc) = genomes_and_contigs_option {
Box::new(move |contig: &str| gc.genome_of_contig(&contig.to_string()).cloned())
} else {
unreachable!("A genome definition is required when using --gff in genome mode")
};
coverm::genes::gene_coverage(
bam_generators,
&mut estimators_and_taker.taker,
&mut estimators_and_taker.estimators,
&gene_definitions,
Some(genome_namer.as_ref()),
print_zeros,
&flag_filter,
threads,
)
} else {
match separator.is_some() || single_genome {
true => coverm::genome::mosdepth_genome_coverage(
bam_generators,
separator.unwrap(),
&mut estimators_and_taker.taker,
print_zeros,
&mut estimators_and_taker.estimators,
&flag_filter,
single_genome,
threads,
),
false => match genomes_and_contigs_option {
Some(gc) => coverm::genome::mosdepth_genome_coverage_with_contig_names(
bam_generators,
gc,
&mut estimators_and_taker.taker,
print_zeros,
&flag_filter,
&mut estimators_and_taker.estimators,
threads,
),
None => unreachable!(),
},
}
};
debug!("Finalising printing ..");
estimators_and_taker.printer.finalise_printing(
&estimators_and_taker.taker,
print_stream,
Some(&reads_mapped),
&estimators_and_taker.columns_to_normalise,
estimators_and_taker.rpkm_column,
estimators_and_taker.tpm_column,
);
}
fn doing_metabat(m: &clap::ArgMatches) -> bool {
let methods: Vec<&str> = m
.get_many::<String>("methods")
.unwrap()
.map(|x| &**x)
.collect();
if methods.contains(&"metabat") {
if methods.len() > 1 {
error!("Cannot specify the metabat method with any other coverage methods");
process::exit(1);
} else {
return true;
}
}
false
}
#[derive(Debug)]
struct FilterParameters {
flag_filters: FlagFilter,
min_aligned_length_single: u32,
min_percent_identity_single: f32,
min_aligned_percent_single: f32,
min_mapq: u8, min_aligned_length_pair: u32,
min_percent_identity_pair: f32,
min_aligned_percent_pair: f32,
}
impl FilterParameters {
pub fn generate_from_clap(m: &clap::ArgMatches) -> FilterParameters {
let f = FilterParameters {
flag_filters: FlagFilter {
include_improper_pairs: !m.get_flag("proper-pairs-only"),
include_secondary: m.get_flag("include-secondary"),
include_supplementary: !m.get_flag("exclude-supplementary"),
},
min_aligned_length_single: *m.get_one::<u32>("min-read-aligned-length").unwrap_or(&0),
min_percent_identity_single: parse_percentage(m, "min-read-percent-identity"),
min_aligned_percent_single: parse_percentage(m, "min-read-aligned-percent"),
min_mapq: *m.get_one::<u8>("min-mapq").unwrap_or(&255),
min_aligned_length_pair: *m
.get_one::<u32>("min-read-aligned-length-pair")
.unwrap_or(&0),
min_percent_identity_pair: parse_percentage(m, "min-read-percent-identity-pair"),
min_aligned_percent_pair: parse_percentage(m, "min-read-aligned-percent-pair"),
};
debug!("Filter parameters set as {f:?}");
f
}
pub fn add_metabat_filtering_if_required(&mut self, m: &clap::ArgMatches) {
if doing_metabat(m) {
info!(
"Setting single read percent identity threshold at 0.97 for \
MetaBAT adjusted coverage, and not filtering out supplementary, \
secondary and improper pair alignments"
);
self.min_percent_identity_single = 0.97001;
self.flag_filters.include_improper_pairs = true;
self.flag_filters.include_supplementary = true;
self.flag_filters.include_secondary = true;
}
}
pub fn doing_filtering(&self) -> bool {
self.min_percent_identity_single > 0.0
|| self.min_percent_identity_pair > 0.0
|| self.min_aligned_percent_single > 0.0
|| self.min_mapq < 255
|| self.min_aligned_percent_pair > 0.0
|| self.min_aligned_length_single > 0
|| self.min_aligned_length_pair > 0
}
}
fn get_sharded_bam_readers<'a, 'b, T>(
m: &'a clap::ArgMatches,
mapping_program: MappingProgram,
reference_tempfile: &'a Option<NamedTempFile>,
genome_exclusion: &'b T,
) -> Vec<ShardedBamReaderGenerator<'b, T>>
where
T: GenomeExclusion,
{
if m.contains_id("cache-unfiltered-bam-directory") {
setup_bam_cache_directory(
m.get_one::<String>("cache-unfiltered-bam-directory")
.unwrap(),
);
}
let discard_unmapped = m.get_flag("discard-unmapped");
let sort_threads = *m.get_one::<u16>("threads").unwrap();
let params = MappingParameters::generate_from_clap(m, mapping_program, reference_tempfile);
let mut bam_readers = vec![];
let mut concatenated_reference_name: Option<String> = None;
let mut concatenated_read_names: Option<String> = None;
for reference_wise_params in params {
let index = setup_mapping_index(&reference_wise_params, m, mapping_program);
let reference = reference_wise_params.reference;
let reference_name = std::path::Path::new(reference)
.file_name()
.expect("Unable to convert reference to file name")
.to_str()
.expect("Unable to covert file name into str")
.to_string();
concatenated_reference_name = match concatenated_reference_name {
Some(prev) => Some(format!("{prev}|{reference_name}")),
None => Some(reference_name),
};
let mut bam_file_cache = build_bam_file_cache_fn(m, reference_tempfile, reference);
for p in reference_wise_params {
bam_readers.push(
coverm::shard_bam_reader::generate_named_sharded_bam_readers_from_reads(
mapping_program,
index.as_ref(),
p.read1,
p.read2,
p.read_format.clone(),
p.threads,
bam_file_cache(p.read1, p.read2)
.as_ref()
.map(String::as_ref),
discard_unmapped,
p.mapping_options,
),
);
let name = &std::path::Path::new(p.read1)
.file_name()
.expect("Unable to convert read1 name to file name")
.to_str()
.expect("Unable to covert file name into str")
.to_string();
concatenated_read_names = match concatenated_read_names {
Some(prev) => Some(format!("{prev}|{name}")),
None => Some(name.to_string()),
};
}
debug!("Finished BAM setup");
}
let gen = ShardedBamReaderGenerator {
stoit_name: format!(
"{}/{}",
concatenated_reference_name.unwrap(),
concatenated_read_names.unwrap()
),
read_sorted_bam_readers: bam_readers,
sort_threads,
genome_exclusion,
};
vec![gen]
}
fn get_streamed_bam_readers(
m: &clap::ArgMatches,
mapping_program: MappingProgram,
reference_tempfile: &Option<NamedTempFile>,
) -> Vec<BamGeneratorSet<StreamingNamedBamReaderGenerator>> {
if m.contains_id("cache-unfiltered-bam-directory") {
setup_bam_cache_directory(
m.get_one::<String>("cache-unfiltered-bam-directory")
.unwrap(),
);
}
let discard_unmapped = m.get_flag("discard-unmapped");
let params = MappingParameters::generate_from_clap(m, mapping_program, reference_tempfile);
let mut generator_set = vec![];
for reference_wise_params in params {
let mut bam_readers = vec![];
let index = setup_mapping_index(&reference_wise_params, m, mapping_program);
let reference = reference_wise_params.reference;
let mut bam_file_cache = build_bam_file_cache_fn(m, reference_tempfile, reference);
for p in reference_wise_params {
bam_readers.push(
coverm::bam_generator::generate_named_bam_readers_from_reads(
mapping_program,
index.as_ref(),
p.read1,
p.read2,
p.read_format.clone(),
p.threads,
bam_file_cache(p.read1, p.read2)
.as_ref()
.map(String::as_ref),
discard_unmapped,
p.mapping_options,
reference_tempfile.is_none(),
),
);
}
debug!("Finished BAM setup");
let to_return = BamGeneratorSet {
generators: bam_readers,
index,
};
generator_set.push(to_return);
}
generator_set
}
fn generate_cached_bam_file_name(directory: &str, reference: &str, read1_path: &str) -> String {
debug!(
"Constructing BAM file cache name in directory {directory}, reference {reference}, read1_path {read1_path}"
);
std::path::Path::new(directory)
.to_str()
.expect("Unable to convert cache-unfiltered-bam-directory name into str")
.to_string()
+ "/"
+ std::path::Path::new(reference)
.file_name()
.expect("Unable to convert reference to file name")
.to_str()
.expect("Unable to covert file name into str")
+ "."
+ std::path::Path::new(read1_path)
.file_name()
.expect("Unable to convert read1 name to file name")
.to_str()
.expect("Unable to covert file name into str")
+ ".bam"
}
fn setup_bam_cache_directory(cache_directory: &str) {
let path = std::path::Path::new(cache_directory);
if path.is_dir() {
if path
.metadata()
.expect("Unable to read metadata for cache directory")
.permissions()
.readonly()
{
error!(
"Cache directory {cache_directory} does not appear to be writeable, not continuing"
);
process::exit(1);
} else {
info!("Writing output files to already existing directory {cache_directory}")
}
} else {
match path.parent() {
Some(parent) => {
let parent2 = match parent == std::path::Path::new("") {
true => std::path::Path::new("."),
false => parent,
};
if parent2
.canonicalize()
.unwrap_or_else(|_| {
panic!(
"Unable to canonicalize parent of cache directory {}",
cache_directory
)
})
.is_dir()
{
if parent2
.metadata()
.unwrap_or_else(|_| {
panic!(
"Unable to get metadata for parent of cache directory {}",
cache_directory
)
})
.permissions()
.readonly()
{
error!(
"The parent directory of the (currently non-existent) \
cache directory {cache_directory} is not writeable, not continuing"
);
process::exit(1);
} else {
info!("Creating cache directory {cache_directory}");
std::fs::create_dir(path).expect("Unable to create cache directory");
}
} else {
error!(
"The parent directory of the cache directory {cache_directory} does not \
yet exist, so not creating that cache directory, and not continuing."
);
process::exit(1);
}
}
None => {
error!("Cannot create root directory {cache_directory}");
process::exit(1);
}
}
}
let tf_result = tempfile::tempfile_in(path);
if tf_result.is_err() {
error!(
"Failed to create test file in bam cache directory: {}",
tf_result.err().unwrap()
);
process::exit(1);
}
}
fn build_cache_name_iter(m: &clap::ArgMatches) -> Option<std::vec::IntoIter<String>> {
if m.contains_id("cache-unfiltered-bam-files") {
let names: Vec<String> = m
.get_many::<String>("cache-unfiltered-bam-files")
.unwrap()
.map(|s| s.to_string())
.collect();
let single_count = m.get_many::<String>("single").map(|v| v.len()).unwrap_or(0);
let read1_count = m.get_many::<String>("read1").map(|v| v.len()).unwrap_or(0);
let coupled_count = m
.get_many::<String>("coupled")
.map(|v| v.len() / 2)
.unwrap_or(0);
let interleaved_count = m
.get_many::<String>("interleaved")
.map(|v| v.len())
.unwrap_or(0);
let expected = single_count + read1_count + coupled_count + interleaved_count;
if names.len() != expected {
error!(
"--cache-unfiltered-bam-files specified {names_len} names but {expected} read sets were provided",
names_len = names.len()
);
process::exit(1);
}
let mut idx = 0;
let single_names = names[idx..idx + single_count].to_vec();
idx += single_count;
let read1_names = names[idx..idx + read1_count].to_vec();
idx += read1_count;
let coupled_names = names[idx..idx + coupled_count].to_vec();
idx += coupled_count;
let interleaved_names = names[idx..idx + interleaved_count].to_vec();
let mut iter_order = Vec::new();
iter_order.extend(read1_names);
iter_order.extend(coupled_names);
iter_order.extend(interleaved_names);
iter_order.extend(single_names);
Some(iter_order.into_iter())
} else {
None
}
}
fn build_bam_file_cache_fn<'a>(
m: &'a clap::ArgMatches,
reference_tempfile: &'a Option<NamedTempFile>,
reference: &'a str,
) -> impl FnMut(&str, Option<&str>) -> Option<String> + 'a {
let mut bam_cache_name_iter = build_cache_name_iter(m);
move |read1: &str, read2: Option<&str>| -> Option<String> {
if let Some(iter) = bam_cache_name_iter.as_mut() {
let name = iter.next().unwrap_or_else(|| {
error!("Not enough BAM file cache names specified");
process::exit(1);
});
match read2 {
Some(r2) => info!("Caching BAM file to {name} for readset {read1} {r2}"),
None => info!("Caching BAM file to {name} for readset {read1}"),
}
Some(name)
} else if m.contains_id("cache-unfiltered-bam-directory") {
let path = generate_cached_bam_file_name(
m.get_one::<String>("cache-unfiltered-bam-directory")
.unwrap(),
match reference_tempfile {
Some(_) => CONCATENATED_REFERENCE_CACHE_STEM,
None => reference,
},
read1,
);
match read2 {
Some(r2) => info!("Caching BAM file to {path} for readset {read1} {r2}"),
None => info!("Caching BAM file to {path} for readset {read1}"),
}
Some(path)
} else {
None
}
}
}
fn get_streamed_filtered_bam_readers(
m: &clap::ArgMatches,
mapping_program: MappingProgram,
reference_tempfile: &Option<NamedTempFile>,
filter_params: &FilterParameters,
) -> Vec<BamGeneratorSet<StreamingFilteredNamedBamReaderGenerator>> {
if m.contains_id("cache-unfiltered-bam-directory") {
setup_bam_cache_directory(
m.get_one::<String>("cache-unfiltered-bam-directory")
.unwrap(),
);
}
let discard_unmapped = m.get_flag("discard-unmapped");
let params = MappingParameters::generate_from_clap(m, mapping_program, reference_tempfile);
let mut generator_set = vec![];
for reference_wise_params in params {
let mut bam_readers = vec![];
let index = setup_mapping_index(&reference_wise_params, m, mapping_program);
let reference = reference_wise_params.reference;
let mut bam_file_cache = build_bam_file_cache_fn(m, reference_tempfile, reference);
for p in reference_wise_params {
bam_readers.push(
coverm::bam_generator::generate_filtered_named_bam_readers_from_reads(
mapping_program,
index.as_ref(),
p.read1,
p.read2,
p.read_format.clone(),
p.threads,
bam_file_cache(p.read1, p.read2)
.as_ref()
.map(String::as_ref),
filter_params.flag_filters.clone(),
filter_params.min_aligned_length_single,
filter_params.min_percent_identity_single,
filter_params.min_aligned_percent_single,
filter_params.min_mapq,
filter_params.min_aligned_length_pair,
filter_params.min_percent_identity_pair,
filter_params.min_aligned_percent_pair,
p.mapping_options,
discard_unmapped,
reference_tempfile.is_none(),
),
);
}
debug!("Finished BAM setup");
let to_return = BamGeneratorSet {
generators: bam_readers,
index,
};
generator_set.push(to_return);
}
generator_set
}
fn run_contig<
R: coverm::bam_generator::NamedBamReader,
T: coverm::bam_generator::NamedBamReaderGenerator<R>,
>(
estimators_and_taker: &mut EstimatorsAndTaker,
bam_readers: Vec<T>,
print_zeros: bool,
flag_filters: FlagFilter,
threads: u16,
gene_definitions: &Option<coverm::genes::GeneDefinitions>,
print_stream: &mut OutputWriter,
) {
let reads_mapped = match gene_definitions {
Some(gene_definitions) => coverm::genes::gene_coverage(
bam_readers,
&mut estimators_and_taker.taker,
&mut estimators_and_taker.estimators,
gene_definitions,
None, print_zeros,
&flag_filters,
threads,
),
None => coverm::contig::contig_coverage(
bam_readers,
&mut estimators_and_taker.taker,
&mut estimators_and_taker.estimators,
print_zeros,
&flag_filters,
threads,
),
};
debug!("Finalising printing ..");
estimators_and_taker.printer.finalise_printing(
&estimators_and_taker.taker,
print_stream,
Some(&reads_mapped),
&estimators_and_taker.columns_to_normalise,
estimators_and_taker.rpkm_column,
estimators_and_taker.tpm_column,
);
}
fn set_log_level(matches: &clap::ArgMatches, is_last: bool) {
set_log_level_bird_tool_utils(matches, is_last, "CoverM", crate_version!());
}