feat: add metagenomic mode support and update gene scoring logic

Author: Ffloriel

orphos-cli/benches/prodigal_benchmarks.rs

@@ -12,6 +12,16 @@ fn run_original_prodigal(
     input_file: &str,
     output_file: &str,
     threads: Option<usize>,
+) -> Result<Duration, Box<dyn std::error::Error>> {
+    run_original_prodigal_with_mode(input_file, output_file, threads, "single")
+}
+
+// Helper function to run original prodigal with specified mode
+fn run_original_prodigal_with_mode(
+    input_file: &str,
+    output_file: &str,
+    threads: Option<usize>,
+    mode: &str,
 ) -> Result<Duration, Box<dyn std::error::Error>> {
     let start = std::time::Instant::now();
 
@@ -23,7 +33,7 @@ fn run_original_prodigal(
         .arg("-f")
         .arg("gff")
         .arg("-p")
-        .arg("single");
+        .arg(mode);
 
     // Original prodigal doesn't have threading, but we include this for consistency
     if threads.is_some() {
@@ -49,6 +59,16 @@ fn run_pyrodigal(
     input_file: &str,
     output_file: &str,
     threads: Option<usize>,
+) -> Result<Duration, Box<dyn std::error::Error>> {
+    run_pyrodigal_with_mode(input_file, output_file, threads, "single")
+}
+
+// Helper function to run pyrodigal with specified mode
+fn run_pyrodigal_with_mode(
+    input_file: &str,
+    output_file: &str,
+    threads: Option<usize>,
+    mode: &str,
 ) -> Result<Duration, Box<dyn std::error::Error>> {
     let start = std::time::Instant::now();
 
@@ -60,7 +80,7 @@ fn run_pyrodigal(
         .arg("-f")
         .arg("gff")
         .arg("-p")
-        .arg("single");
+        .arg(mode);
 
     if let Some(t) = threads {
         cmd.arg("-j").arg(t.to_string());
@@ -85,6 +105,16 @@ fn run_orphos_cli(
     input_file: &str,
     output_file: &str,
     threads: Option<usize>,
+) -> Result<Duration, Box<dyn std::error::Error>> {
+    run_orphos_cli_with_mode(input_file, output_file, threads, "single")
+}
+
+// Helper function to run orphos with specified mode
+fn run_orphos_cli_with_mode(
+    input_file: &str,
+    output_file: &str,
+    threads: Option<usize>,
+    mode: &str,
 ) -> Result<Duration, Box<dyn std::error::Error>> {
     let start = std::time::Instant::now();
 
@@ -106,7 +136,7 @@ fn run_orphos_cli(
         .arg("-f")
         .arg("gff")
         .arg("-p")
-        .arg("single");
+        .arg(mode);
 
     let output = cmd.output()?;
     let duration = start.elapsed();
@@ -426,11 +456,109 @@ fn benchmark_scaling_analysis(c: &mut Criterion) {
     group.finish();
 }
 
+/// Benchmark meta mode (metagenomic) - single threaded comparison
+fn benchmark_meta_mode_single_threaded(c: &mut Criterion) {
+    let test_files = [
+        ("ecoli", "tests/data/ecoli.fasta"),
+        ("salmonella", "tests/data/salmonella.fasta"),
+    ];
+
+    for (name, input_file) in test_files {
+        if !Path::new(input_file).exists() {
+            eprintln!("Warning: {} not found, skipping benchmark", input_file);
+            continue;
+        }
+
+        let file_size = get_file_size(input_file);
+        let mut group = c.benchmark_group(format!("{}_meta_single_threaded", name));
+        group.throughput(Throughput::Bytes(file_size));
+
+        // Check if original prodigal is available
+        if Command::new("prodigal").arg("--help").output().is_ok() {
+            group.bench_function("original_prodigal_meta", |b| {
+                b.iter_custom(|iters| {
+                    let mut total_duration = Duration::new(0, 0);
+                    for _ in 0..iters {
+                        let output_file = NamedTempFile::new().unwrap();
+                        let duration = run_original_prodigal_with_mode(
+                            black_box(input_file),
+                            output_file.path().to_str().unwrap(),
+                            None,
+                            "meta",
+                        )
+                        .unwrap_or_else(|e| {
+                            eprintln!("Original prodigal meta benchmark failed: {}", e);
+                            Duration::from_secs(0)
+                        });
+                        total_duration += duration;
+                    }
+                    total_duration
+                });
+            });
+        } else {
+            eprintln!(
+                "Warning: Original prodigal not found, skipping original_prodigal_meta benchmark"
+            );
+        }
+
+        // Check if pyrodigal is available
+        if Command::new("pyrodigal").arg("--help").output().is_ok() {
+            group.bench_function("pyrodigal_meta", |b| {
+                b.iter_custom(|iters| {
+                    let mut total_duration = Duration::new(0, 0);
+                    for _ in 0..iters {
+                        let output_file = NamedTempFile::new().unwrap();
+                        let duration = run_pyrodigal_with_mode(
+                            black_box(input_file),
+                            output_file.path().to_str().unwrap(),
+                            Some(1),
+                            "meta",
+                        )
+                        .unwrap_or_else(|e| {
+                            eprintln!("Pyrodigal meta benchmark failed: {}", e);
+                            Duration::from_secs(0)
+                        });
+                        total_duration += duration;
+                    }
+                    total_duration
+                });
+            });
+        } else {
+            eprintln!("Warning: Pyrodigal not found, skipping pyrodigal_meta benchmark");
+        }
+
+        // Benchmark orphos meta mode (single-threaded)
+        group.bench_function("orphos_meta", |b| {
+            b.iter_custom(|iters| {
+                let mut total_duration = Duration::new(0, 0);
+                for _ in 0..iters {
+                    let output_file = NamedTempFile::new().unwrap();
+                    let duration = run_orphos_cli_with_mode(
+                        black_box(input_file),
+                        output_file.path().to_str().unwrap(),
+                        Some(1),
+                        "meta",
+                    )
+                    .unwrap_or_else(|e| {
+                        eprintln!("orphos-cli meta benchmark failed: {}", e);
+                        Duration::from_secs(0)
+                    });
+                    total_duration += duration;
+                }
+                total_duration
+            });
+        });
+
+        group.finish();
+    }
+}
+
 criterion_group!(
     name = benches;
     config = configure_criterion();
     targets = benchmark_single_threaded,
     // benchmark_multi_threaded,
-    benchmark_scaling_analysis
+    benchmark_scaling_analysis,
+    benchmark_meta_mode_single_threaded
 );
 criterion_main!(benches);

orphos-cli/tests/snapshots/meta_genomes_equivalence__meta_genomes_gff_summary.snap

@@ -3,16 +3,16 @@ source: orphos-cli/tests/meta_genomes_equivalence.rs
 expression: summary
 ---
 file	similarity(%)	coord_match(%)	exact	orig_hash	rust_hash
-NC_000913.3.fasta	5.957	99.652	no	9995c9861491	624c54e087bd
-NC_000915.1.fasta	0.829	99.171	no	661adac5bfad	3f661f42ee96
-NC_000962.3.fasta	9.641	99.799	no	5a26d7cdd39c	77a5b452ab9b
-NC_000964.3.fasta	1.017	99.410	no	80bd80b02c2e	dbf14d16e018
-NC_001318.1.fasta	15.599	100.000	no	fd3c884611ff	0966ac145b56
-NC_002516.2.fasta	16.726	99.930	no	38c34e144238	dac532d8e321
-NC_002528.1.fasta	38.384	100.000	no	f7f46f0a7a71	22cc1aff9eae
-NC_002745.2.fasta	3.599	99.846	no	797f59fb9866	4cf1d1ca1aad
-NC_003888.3.fasta	4.550	99.923	no	082e59efe934	0651fa453cde
-NC_005072.1.fasta	13.788	99.736	no	f1398c92c767	94f273ed0d9a
-NC_008512.1.fasta	24.599	100.000	no	927216b98ab4	fcbf39edf768
-fragment.fasta	14.761	99.708	no	17bf514adad9	f1da8e43ec38
-fragment2.fasta	25.926	100.000	no	df7200909f9a	2cb02a1b7e96
+NC_000913.3.fasta	56.452	99.930	no	9995c9861491	41c6dad6f214
+NC_000915.1.fasta	61.900	99.873	no	661adac5bfad	718323dfeee3
+NC_000962.3.fasta	28.998	99.950	no	5a26d7cdd39c	107315ed518a
+NC_000964.3.fasta	9.312	100.000	no	80bd80b02c2e	76f7d8cda2d5
+NC_001318.1.fasta	57.557	100.000	no	fd3c884611ff	333684f6b76e
+NC_002516.2.fasta	20.098	99.965	no	38c34e144238	b5d18dc81ee8
+NC_002528.1.fasta	100.000	100.000	yes	f7f46f0a7a71	f7f46f0a7a71
+NC_002745.2.fasta	12.859	99.962	no	797f59fb9866	2b72e6297f92
+NC_003888.3.fasta	10.750	99.923	no	082e59efe934	f68bb4b2a257
+NC_005072.1.fasta	13.664	100.000	no	f1398c92c767	11103121e304
+NC_008512.1.fasta	100.000	100.000	yes	927216b98ab4	927216b98ab4
+fragment.fasta	99.422	100.000	no	17bf514adad9	18f398b80b81
+fragment2.fasta	100.000	100.000	yes	df7200909f9a	df7200909f9a

orphos-cli/tests/snapshots/meta_genomes_equivalence__meta_mode_individual_files_summary.snap

@@ -0,0 +1,7 @@
+---
+source: orphos-cli/tests/meta_genomes_equivalence.rs
+expression: summary
+---
+file	similarity(%)	coord_match(%)	exact	orig_hash	rust_hash
+ecoli.fasta	56.452	99.930	no	00a4ae70baa3	73394cf5247e
+salmonella.fasta	17.403	99.957	no	1f477c1793a5	a25cdf8bc731

orphos-core/src/algorithms/gene_finding/annotation.rs

@@ -3,7 +3,7 @@ use bio::bio_types::strand::Strand;
 use crate::{
     constants::{MAX_CONFIDENCE_SCORE, RBS_DESCRIPTIONS},
     sequence::mer_text,
-    types::{Gene, GeneAnnotation, GeneScore, Node, StartType, Training},
+    types::{CodonType, Gene, GeneAnnotation, GeneScore, Node, StartType, Training},
 };
 
 /// Record gene annotation and score.
@@ -125,6 +125,48 @@ fn calculate_confidence(score: f64, start_weight: f64) -> f64 {
     }
 }
 
+/// Update gene scores from re-scored nodes.
+///
+/// This function is used in metagenomic mode to update gene scores after the final
+/// re-scoring with the best model. It finds the corresponding start node by matching
+/// the gene's begin/end positions and strand, then updates the score fields.
+///
+/// # Arguments
+/// * `genes` - Mutable slice of genes to update
+/// * `nodes` - Slice of re-scored nodes
+/// * `training` - Training parameters for confidence calculation
+pub fn update_gene_scores_from_nodes(genes: &mut [Gene], nodes: &[Node], training: &Training) {
+    for gene in genes.iter_mut() {
+        // Find the start node by matching position and strand
+        // For forward strand: start_node.index + 1 == gene.begin
+        // For reverse strand: start_node.index + 1 == gene.end
+        let target_pos = if gene.coordinates.strand == Strand::Forward {
+            gene.coordinates.begin.saturating_sub(1)
+        } else {
+            gene.coordinates.end.saturating_sub(1)
+        };
+
+        if let Some(start_node) = nodes.iter().find(|n| {
+            n.position.index == target_pos
+                && n.position.strand == gene.coordinates.strand
+                && n.position.codon_type != CodonType::Stop
+        }) {
+            let total_score = start_node.scores.coding_score + start_node.scores.start_score;
+            let confidence = calculate_confidence(total_score, training.start_weight_factor);
+
+            gene.score = GeneScore {
+                confidence,
+                total_score,
+                coding_score: start_node.scores.coding_score,
+                start_score: start_node.scores.start_score,
+                ribosome_binding_score: start_node.scores.ribosome_binding_score,
+                upstream_score: start_node.scores.upstream_score,
+                type_score: start_node.scores.type_score,
+            };
+        }
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -186,6 +228,7 @@ mod tests {
             },
             score: GeneScore::default(),
             annotation: GeneAnnotation::default(),
+            display_score: None,
         }
     }
 

orphos-core/src/algorithms/gene_finding/mod.rs

@@ -3,4 +3,5 @@ pub mod creation;
 pub mod optimization;
 pub mod scoring;
 
+pub use annotation::update_gene_scores_from_nodes;
 pub use optimization::GeneBuilder;