add process subcommand that runs analyse and cluster

Author: AroneyS

src/analyse.rs

@@ -25,6 +25,7 @@ pub fn analyse<Q: QualityFinder, R: RrnaFinder, T: TrnaFinder>(
     rrna_finder: &R,
     trna_finder: &T,
     checkm2_quality_report: &Option<String>,
+    output_quality_report_path: &Option<String>,
     checkm_tab_table: &Option<String>,
     barrnap_gff_list: &Option<String>,
     trnascan_out_list: &Option<String>,
@@ -119,6 +120,27 @@ pub fn analyse<Q: QualityFinder, R: RrnaFinder, T: TrnaFinder>(
         cache
     } else {
         quality_finder.prepare_comp_cont(genomes, threads, tmp_path);
+
+        // If requested, copy the CheckM2 quality report out of the temporary directory now
+        if let Some(dest) = output_quality_report_path {
+            let src = tmp_path.join("checkm2").join("quality_report.tsv");
+            if let Some(parent) = std::path::Path::new(dest).parent() {
+                if !parent.as_os_str().is_empty() {
+                    std::fs::create_dir_all(parent).map_err(|e| {
+                        format!("Failed to create parent directory for quality report output: {e}")
+                    })?;
+                }
+            }
+            std::fs::copy(&src, dest).map_err(|e| {
+                format!(
+                    "Failed to copy CheckM2 quality report from {} to {}: {}",
+                    src.display(),
+                    dest,
+                    e
+                )
+            })?;
+        }
+
         genomes
             .iter()
             .map(|g| (g.clone(), quality_finder.find_comp_cont(g)))
@@ -219,7 +241,7 @@ fn parse_barrnap_gff_list(
     Ok(rrna_cache)
 }
 
-/// Parse two-column TSV mapping genome names to tRNAscan-SE output files  
+/// Parse two-column TSV mapping genome names to tRNAscan-SE output files
 fn parse_trnascan_out_list(list_path: &str) -> Result<HashMap<String, usize>, String> {
     let mut trna_cache = HashMap::new();
     let content = std::fs::read_to_string(list_path)

src/analyse_argument_parsing.rs

@@ -88,14 +88,18 @@ pub struct GalahAnalyser<'a> {
 }
 
 impl GalahAnalyser<'_> {
-    pub fn analyse(&mut self) -> Result<std::collections::HashMap<String, GenomeOutput>, String> {
+    pub fn analyse(
+        &mut self,
+        output_quality_report_path: &Option<String>,
+    ) -> Result<std::collections::HashMap<String, GenomeOutput>, String> {
         crate::analyse::analyse(
             self.genome_fasta_files,
             self.threads,
             &mut self.quality_analyser,
             &self.rrna_analyser,
             &self.trna_analyser,
             &self.checkm2_quality_report,
+            output_quality_report_path,
             &self.checkm_tab_table,
             &self.barrnap_gff_list,
             &self.trnascan_out_list,
@@ -108,6 +112,7 @@ pub struct GalahAnalyserCommandDefinition {
     pub rrna_method_argument: String,
     pub trna_method_argument: String,
     pub output_mimag_summary_argument: String,
+    pub output_quality_report_argument: String,
     pub checkm2_db_path_argument: String,
     pub checkm2_quality_report_argument: String,
     pub checkm_tab_table_argument: String,
@@ -116,12 +121,13 @@ pub struct GalahAnalyserCommandDefinition {
 }
 
 lazy_static! {
-    static ref ANALYSE_COMMAND_DEFINITION: GalahAnalyserCommandDefinition = {
+    pub static ref ANALYSE_COMMAND_DEFINITION: GalahAnalyserCommandDefinition = {
         GalahAnalyserCommandDefinition {
             quality_method_argument: "quality-method".to_string(),
             rrna_method_argument: "rrna-method".to_string(),
             trna_method_argument: "trna-method".to_string(),
             output_mimag_summary_argument: "output-mimag-summary".to_string(),
+            output_quality_report_argument: "output-quality-report".to_string(),
             checkm2_db_path_argument: "checkm2-db-path".to_string(),
             checkm2_quality_report_argument: "checkm2-quality-report".to_string(),
             checkm_tab_table_argument: "checkm-tab-table".to_string(),
@@ -314,7 +320,21 @@ pub fn add_analyse_subcommand(app: clap::Command) -> clap::Command {
             Arg::new(&*ANALYSE_COMMAND_DEFINITION.output_mimag_summary_argument)
                 .long("output-mimag-summary")
                 .value_name("SUMMARY")
-                .help("Path to output MIMAG summary file"),
+                .help("Path to output MIMAG summary file")
+            .required_unless_present_any([
+                "output-quality-report",
+                "full-help",
+                "full-help-roff",]),
+        )
+        .arg(
+            Arg::new(&*ANALYSE_COMMAND_DEFINITION.output_quality_report_argument)
+                .long("output-quality-report")
+                .value_name("REPORT")
+                .help("Path to output CheckM2-format quality report")
+            .required_unless_present_any([
+                "output-mimag-summary",
+                "full-help",
+                "full-help-roff",]),
         )
         .arg(
             Arg::new(&*ANALYSE_COMMAND_DEFINITION.rrna_method_argument)
@@ -460,15 +480,22 @@ pub fn add_analyse_output_parameters_to_section(
     section: Section,
     definition: &GalahAnalyserCommandDefinition,
 ) -> Section {
-    section.option(
-        Opt::new("PATH")
-            .long(&format!("--{}", definition.output_mimag_summary_argument))
-            .help("Output a tsv file summarising the MIMAG status for each genome."),
-    )
+    section
+        .option(
+            Opt::new("PATH")
+                .long(&format!("--{}", definition.output_mimag_summary_argument))
+                .help("Output a tsv file summarising the MIMAG status for each genome."),
+        )
+        .option(
+            Opt::new("PATH")
+                .long(&format!("--{}", definition.output_quality_report_argument))
+                .help("Output a CheckM2-format quality report TSV file."),
+        )
 }
 
 pub struct AnalyseOutput {
     pub output_mimag_summary: Option<std::fs::File>,
+    pub output_quality_report_path: Option<String>,
 }
 
 pub fn setup_analyse_outputs(
@@ -479,8 +506,13 @@ pub fn setup_analyse_outputs(
         .get_one::<String>(&command_definition.output_mimag_summary_argument)
         .map(|o| std::fs::File::create(o).expect("Failed to open output MIMAG summary file"));
 
+    let output_quality_report_path = m
+        .get_one::<String>(&command_definition.output_quality_report_argument)
+        .map(|s| s.to_string());
+
     AnalyseOutput {
         output_mimag_summary,
+        output_quality_report_path,
     }
 }
 
@@ -511,9 +543,11 @@ pub fn run_analyse_subcommand(
     let output_definitions = setup_analyse_outputs(m, &ANALYSE_COMMAND_DEFINITION);
 
     info!("Analysing {} genomes ..", genome_fasta_files.len());
-    let analysis = galah.analyse().expect("Failed to analyse genomes");
+    let analysis = galah
+        .analyse(&output_definitions.output_quality_report_path)
+        .expect("Failed to analyse genomes");
 
-    write_analyse_outputs(output_definitions, &analysis, genome_fasta_files);
+    write_analyse_outputs(output_definitions, &analysis, &genome_fasta_files);
     info!("Finished printing genome analysis");
 }
 
@@ -580,10 +614,10 @@ fn generate_galah_analyser<'a>(
     })
 }
 
-fn write_analyse_outputs(
+pub fn write_analyse_outputs(
     output_definitions: AnalyseOutput,
     analysis: &HashMap<String, GenomeOutput>,
-    genome_fasta_files: Vec<String>,
+    genome_fasta_files: &Vec<String>,
 ) {
     if let Some(mut f) = output_definitions.output_mimag_summary {
         writeln!(
@@ -592,7 +626,7 @@ fn write_analyse_outputs(
         )
         .unwrap();
         for genome in genome_fasta_files {
-            if let Some(output_data) = analysis.get(&*genome) {
+            if let Some(output_data) = analysis.get(&**genome) {
                 writeln!(
                     f,
                     "{genome}\t{completeness:.2}\t{contamination:.2}\t{r5s}\t{r16s}\t{r23s}\t{trnas}\t{mimag_quality}",

src/checkm2.rs

@@ -12,20 +12,36 @@ pub struct CheckM2Analyser {
     // Cache for completeness and contamination results
     pub comp_cont_cache: HashMap<String, (f64, f64)>,
     pub database_path: String,
+    pub quality_report_source_path: Option<std::path::PathBuf>,
 }
 
 impl CheckM2Analyser {
     pub fn new(database_path: String) -> Self {
         Self {
             comp_cont_cache: HashMap::new(),
             database_path,
+            quality_report_source_path: None,
+        }
+    }
+
+    pub fn copy_quality_report(&self, dest_path: &str) -> Result<(), String> {
+        if let Some(src_path) = &self.quality_report_source_path {
+            std::fs::copy(src_path, dest_path).map_err(|e| {
+                format!("Failed to copy quality report from {src_path:?} to {dest_path}: {e}")
+            })?;
+            Ok(())
+        } else {
+            Err("No quality report available to copy (CheckM2 may not have been run)".to_string())
         }
     }
 }
 
 impl QualityFinder for CheckM2Analyser {
     fn prepare_comp_cont(&mut self, genome_paths: &[String], threads: usize, tmp_path: &Path) {
-        self.comp_cont_cache = get_comp_cont(genome_paths, threads, tmp_path, &self.database_path);
+        let (cache, quality_report_path) =
+            get_comp_cont(genome_paths, threads, tmp_path, &self.database_path);
+        self.comp_cont_cache = cache;
+        self.quality_report_source_path = Some(quality_report_path);
     }
 
     fn find_comp_cont(&self, genome_path: &str) -> (f64, f64) {
@@ -45,7 +61,7 @@ fn get_comp_cont(
     threads: usize,
     tmp_path: &Path,
     database_path: &str,
-) -> HashMap<String, (f64, f64)> {
+) -> (HashMap<String, (f64, f64)>, std::path::PathBuf) {
     let mut comp_cont_cache = HashMap::new();
     let checkm2_path = tmp_path.join("checkm2");
 
@@ -136,5 +152,5 @@ fn get_comp_cont(
             );
         }
     }
-    comp_cont_cache
+    (comp_cont_cache, quality_report_path)
 }

src/cluster_argument_parsing.rs

@@ -182,13 +182,13 @@ lazy_static! {
 
 {}
 
-  {} cluster --genome-fasta-directory input_genomes/ 
+  {} cluster --genome-fasta-directory input_genomes/
     --output-representative-fasta-directory output_directory/
 
 {}
 
   {} cluster --ani 95 --precluster-ani 90 --precluster-method finch
-    --genome-fasta-list genomes.txt 
+    --genome-fasta-list genomes.txt
     --output-cluster-definition clusters.tsv
 
 {}
@@ -199,7 +199,7 @@ lazy_static! {
 
 {}
 
-  {} cluster --cluster-contigs --small-contigs --genome-fasta-files contigs.fasta 
+  {} cluster --cluster-contigs --small-contigs --genome-fasta-files contigs.fasta
     --output-cluster-definition contig_clusters.tsv
 
 See {} cluster --full-help for further options and further detail.
@@ -308,7 +308,7 @@ pub fn add_dereplication_clustering_parameters_to_section(
         .option(
             Opt::new("FLOAT")
                 .long(&format!("--{}", definition.dereplication_ani_argument))
-                .help(&format!("Overall ANI level to dereplicate at with the primary clusterer. {}", 
+                .help(&format!("Overall ANI level to dereplicate at with the primary clusterer. {}",
                     &default_roff(crate::DEFAULT_ANI))),
         )
         .option(
@@ -579,9 +579,11 @@ pub fn run_cluster_subcommand(
     }
 
     let contig_names_owned = if cluster_contigs {
-        if m.contains_id("output-representative-fasta-directory")
-            || m.contains_id("output-representative-fasta-directory-copy")
-        {
+        if m.contains_id(
+            &GALAH_COMMAND_DEFINITION.dereplication_output_representative_fasta_directory,
+        ) || m.contains_id(
+            &GALAH_COMMAND_DEFINITION.dereplication_output_representative_fasta_directory_copy,
+        ) {
             panic!("Cannot specify --cluster-contigs with --output-representative-fasta-directory or --output-representative-fasta-directory-copy");
         }
 
@@ -674,6 +676,7 @@ pub fn run_cluster_subcommand(
         m,
         &GALAH_COMMAND_DEFINITION,
         ref_genomes_for_clusterer.as_deref(),
+        None,
     )
     .expect("Failed to parse galah clustering arguments correctly");
 
@@ -844,6 +847,7 @@ pub fn filter_genomes_through_checkm<'a>(
     genome_fasta_files: &'a Vec<String>,
     clap_matches: &clap::ArgMatches,
     argument_definition: &GalahClustererCommandDefinition,
+    injected_quality_report: Option<String>,
 ) -> std::result::Result<Vec<&'a str>, String> {
     if clap_matches.get_flag(&argument_definition.dereplication_cluster_contigs_argument) {
         return Ok(genome_fasta_files.iter().map(|s| &**s).collect());
@@ -852,7 +856,9 @@ pub fn filter_genomes_through_checkm<'a>(
     match clap_matches.contains_id("checkm-tab-table")
         || clap_matches.contains_id("genome-info")
         || clap_matches.contains_id("checkm2-quality-report")
-        || (clap_matches.contains_id("run-checkm2") && clap_matches.get_flag("run-checkm2"))
+        || injected_quality_report.is_some()
+        || (clap_matches.contains_id(&argument_definition.dereplication_run_checkm2_argument)
+            && clap_matches.get_flag(&argument_definition.dereplication_run_checkm2_argument))
     {
         false => {
             warn!("Since CheckM input has not been provided and CheckM2 has been disabled, genomes are not being ordered by quality. Instead the order of their input is being used");
@@ -897,11 +903,20 @@ pub fn filter_genomes_through_checkm<'a>(
                     )
                     .expect("Error parsing genomeInfo file"),
                 }
-            } else if clap_matches.contains_id("run-checkm2")
-                && clap_matches.get_flag("run-checkm2")
+            } else if injected_quality_report.is_some() {
+                info!("Reading injected CheckM2 Quality report ..");
+                CheckMResultEnum::CheckM2Result {
+                    result: checkm::CheckM2QualityReport::read_file_path(
+                        injected_quality_report.as_deref().unwrap(),
+                    )
+                    .unwrap(),
+                }
+            } else if clap_matches
+                .contains_id(&argument_definition.dereplication_run_checkm2_argument)
+                && clap_matches.get_flag(&argument_definition.dereplication_run_checkm2_argument)
             {
                 // Run CheckM2 as in analyse
-                let db_path = clap_matches.get_one::<String>("checkm2-db-path")
+                let db_path = clap_matches.get_one::<String>(&argument_definition.dereplication_checkm2_db_path_argument)
                     .map(|s| s.to_string())
                     .or_else(|| std::env::var("CHECKM2DB").ok())
                     .expect("CheckM2 database path must be provided via --checkm2-db-path or CHECKM2DB env var");
@@ -1190,6 +1205,7 @@ pub fn generate_galah_clusterer<'a>(
     clap_matches: &clap::ArgMatches,
     argument_definition: &GalahClustererCommandDefinition,
     reference_genomes: Option<&[&str]>,
+    injected_quality_report: Option<String>,
 ) -> std::result::Result<GalahClusterer<'a>, String> {
     crate::external_command_checker::check_for_fastani();
 
@@ -1204,7 +1220,12 @@ pub fn generate_galah_clusterer<'a>(
                 .as_str()
     };
 
-    match filter_genomes_through_checkm(genome_fasta_paths, clap_matches, argument_definition) {
+    match filter_genomes_through_checkm(
+        genome_fasta_paths,
+        clap_matches,
+        argument_definition,
+        injected_quality_report,
+    ) {
         Err(e) => std::result::Result::Err(e),
 
         Ok(v2) => {

src/lib.rs

@@ -10,6 +10,8 @@ pub mod fastani;
 pub mod finch;
 pub mod genome_info_file;
 pub mod genome_stats;
+pub mod process;
+pub mod process_argument_parsing;
 pub mod skani;
 pub mod sorted_pair_genome_distance_cache;
 pub mod trnascan;

src/main.rs

@@ -53,6 +53,13 @@ fn main() {
                 crate_version!(),
             );
         }
+        Some("process") => {
+            galah::process_argument_parsing::run_process_subcommand(
+                &matches,
+                "galah",
+                crate_version!(),
+            );
+        }
         _ => panic!("Programming error"),
     }
 }
@@ -122,5 +129,6 @@ fn build_cli() -> Command {
 
     app = galah::cluster_argument_parsing::add_cluster_subcommand(app);
     app = galah::analyse_argument_parsing::add_analyse_subcommand(app);
+    app = galah::process_argument_parsing::add_process_subcommand(app);
     app
 }