use std::{collections::HashSet, process::Command};
use anyhow::Result;
use log::{info, debug};
use crate::coverage::{coverage_calculator::ReadCollection, coverage_table::CoverageTable};
pub struct CovermEngine<'a> {
assembly: &'a str,
threads: usize,
settings: CovermSettings,
}
impl<'a> CovermEngine<'a> {
pub fn new(m: &'a clap::ArgMatches) -> Result<Self> {
let assembly = m.get_one::<String>("assembly").unwrap();
let output_directory = m.get_one::<String>("output-directory").unwrap();
std::fs::create_dir_all(&output_directory)?;
let settings = CovermSettings::new(m)?;
Ok(
Self {
assembly,
threads: *m.get_one::<usize>("threads").unwrap(),
settings
}
)
}
pub fn run(&self, samples_names_to_run: HashSet<&str>, read_collection: &ReadCollection) -> Result<CoverageTable> {
debug!("Sample names to run: {:?}", samples_names_to_run);
let mut coverage_tables = Vec::with_capacity(4);
let short_reads_to_map = read_collection.subset_short_reads(&samples_names_to_run);
if short_reads_to_map.len() > 0 {
info!("Mapping {} short reads.", short_reads_to_map.len());
coverage_tables.push(self.run_coverm(short_reads_to_map, MappingMode::ShortRead)?);
}
let long_reads_to_map = read_collection.subset_long_reads(&samples_names_to_run);
if long_reads_to_map.len() > 0 {
info!("Mapping {} long reads.", long_reads_to_map.len());
coverage_tables.push(self.run_coverm(long_reads_to_map, MappingMode::LongRead)?);
}
let short_bams_to_use = read_collection.subset_short_read_bams(&samples_names_to_run);
if short_bams_to_use.len() > 0 {
info!("Calculating coverage for {} short read bams.", short_bams_to_use.len());
coverage_tables.push(self.run_coverm(short_bams_to_use, MappingMode::ShortBam)?);
}
let long_bams_to_use = read_collection.subset_long_read_bams(&samples_names_to_run);
if long_bams_to_use.len() > 0 {
info!("Calculating coverage for {} long read bams.", long_bams_to_use.len());
coverage_tables.push(self.run_coverm(long_bams_to_use, MappingMode::LongBam)?);
}
CoverageTable::merge_many(coverage_tables)
}
fn run_coverm(&self, read_collection: ReadCollection, mode: MappingMode) -> Result<CoverageTable> {
let mut coverm_command = Command::new("coverm");
coverm_command
.arg("contig")
.arg("--threads").arg(&format!("{}", self.threads))
.arg("--min-covered-fraction").arg(&format!("{}", self.settings.min_covered_fraction));
match mode {
MappingMode::ShortRead => {
coverm_command
.arg("--mapper")
.arg(&self.settings.short_read_mapper)
.arg("--reference")
.arg(&self.assembly);
}
MappingMode::LongRead => {
coverm_command
.arg("--mapper")
.arg(&self.settings.long_read_mapper)
.arg("--reference")
.arg(&self.assembly);
}
_ => {
}
};
if let Some(minimap_params) = &self.settings.minimap_params {
coverm_command
.arg("--minimap2-params")
.arg(minimap_params);
}
if let Some(bwa_params) = &self.settings.bwa_params {
coverm_command
.arg("--bwa-params")
.arg(bwa_params);
}
if let Some(min_read_aligned_length) = &self.settings.min_read_aligned_length {
coverm_command
.arg("--min-read-aligned-length")
.arg(&format!("{}", min_read_aligned_length));
}
if let Some(min_read_percent_identity) = &self.settings.min_read_percent_identity {
coverm_command
.arg("--min-read-percent-identity")
.arg(&format!("{}", min_read_percent_identity));
}
if let Some(min_read_aligned_percent) = &self.settings.min_read_aligned_percent {
coverm_command
.arg("--min-read-aligned-percent")
.arg(&format!("{}", min_read_aligned_percent));
}
if self.settings.include_secondary {
coverm_command
.arg("--include-secondary");
}
if self.settings.exclude_supplementary {
coverm_command
.arg("--exclude-supplementary");
}
coverm_command
.arg("--contig-end-exclusion")
.arg(&format!("{}", self.settings.contig_end_exclusion))
.arg("--trim-min")
.arg(&format!("{}", self.settings.trim_min))
.arg("--trim-max")
.arg(&format!("{}", self.settings.trim_max));
match mode {
MappingMode::ShortRead | MappingMode::ShortBam => {
coverm_command
.arg("--methods")
.arg("metabat");
if let Some(min_read_aligned_length_pair) = &self.settings.min_read_aligned_length_pair {
coverm_command
.arg("--min-read-aligned-length-pair")
.arg(&format!("{}", min_read_aligned_length_pair));
}
if let Some(min_read_percent_identity_pair) = &self.settings.min_read_percent_identity_pair {
coverm_command
.arg("--min-read-percent-identity-pair")
.arg(&format!("{}", min_read_percent_identity_pair));
}
if let Some(min_read_aligned_percent_pair) = &self.settings.min_read_aligned_percent_pair {
coverm_command
.arg("--min-read-aligned-percent-pair")
.arg(&format!("{}", min_read_aligned_percent_pair));
}
if self.settings.proper_pairs_only {
coverm_command
.arg("--proper-pairs-only");
}
}
MappingMode::LongRead | MappingMode::LongBam => {
coverm_command
.arg("--methods")
.arg("length")
.arg("trimmed_mean")
.arg("variance");
}
};
read_collection.add_to_coverm_command(&mut coverm_command);
let temp_file = tempfile::NamedTempFile::new()?;
let temp_file_path = temp_file.path().to_str().unwrap();
coverm_command
.arg("--output-file")
.arg(temp_file_path);
coverm_command
.stdout(std::process::Stdio::piped())
.stderr(std::process::Stdio::piped());
match coverm_command.output() {
Ok(output) => {
if output.status.success() {
let coverage_table = CoverageTable::from_file(temp_file_path, mode)?;
Ok(coverage_table)
} else {
Err(anyhow::anyhow!("Coverm failed with exit code: {} {} {}",
output.status,
std::str::from_utf8(output.stdout.as_slice())?,
std::str::from_utf8(output.stderr.as_slice())?
))
}
}
Err(e) => {
Err(anyhow::anyhow!("Coverm failed with error: {}", e))
}
}
}
}
#[derive(Debug, Clone, Copy)]
pub enum MappingMode {
ShortRead,
LongRead,
ShortBam,
LongBam,
}
struct CovermSettings {
short_read_mapper: String,
long_read_mapper: String,
minimap_params: Option<String>,
bwa_params: Option<String>,
min_read_aligned_length: Option<u32>,
min_read_percent_identity: Option<f32>,
min_read_aligned_percent: Option<f32>,
min_read_aligned_length_pair: Option<u32>,
min_read_percent_identity_pair: Option<f32>,
min_read_aligned_percent_pair: Option<f32>,
proper_pairs_only: bool,
min_covered_fraction: f32,
include_secondary: bool,
exclude_supplementary: bool,
contig_end_exclusion: usize,
trim_min: f32,
trim_max: f32,
}
impl CovermSettings {
pub fn new(m: &clap::ArgMatches) -> Result<Self> {
let minimap_params = m.get_one::<String>("minimap2-params").cloned();
let bwa_params = m.get_one::<String>("bwa-params").cloned();
let min_read_aligned_length = m.get_one::<u32>("min-read-aligned-length").cloned();
let min_read_percent_identity = m.get_one::<f32>("min-read-percent-identity").cloned();
let min_read_aligned_percent = m.get_one::<f32>("min-read-aligned-percent").cloned();
let min_read_aligned_length_pair = m.get_one::<u32>("min-read-aligned-length-pair").cloned();
let min_read_percent_identity_pair = m.get_one::<f32>("min-read-percent-identity-pair").cloned();
let min_read_aligned_percent_pair = m.get_one::<f32>("min-read-aligned-percent-pair").cloned();
let proper_pairs_only = m.get_flag("proper-pairs-only");
let min_covered_fraction = m.get_one::<f32>("min-covered-fraction").unwrap().clone();
let include_secondary = m.get_flag("include-secondary");
let exclude_supplementary = m.get_flag("exclude-supplementary");
let contig_end_exclusion = m.get_one::<usize>("contig-end-exclusion").unwrap().clone();
let trim_min = m.get_one::<f32>("trim-min").unwrap().clone();
let trim_max = m.get_one::<f32>("trim-max").unwrap().clone();
Ok(
Self {
short_read_mapper: m.get_one::<String>("mapper").unwrap().clone(),
long_read_mapper: m.get_one::<String>("longread-mapper").unwrap().clone(),
minimap_params,
bwa_params,
min_read_aligned_length,
min_read_percent_identity,
min_read_aligned_percent,
min_read_aligned_length_pair,
min_read_percent_identity_pair,
min_read_aligned_percent_pair,
proper_pairs_only,
min_covered_fraction,
include_secondary,
exclude_supplementary,
contig_end_exclusion,
trim_min,
trim_max,
}
)
}
}