Add configurable --display-thresh parameter for .prob file filtering

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "oarfish"
-version = "0.9.2"
+version = "0.9.3"
 edition = "2024"
 rust-version = "1.91"
 authors = [
@@ -31,7 +31,7 @@ seqcol_rs = { git = "https://github.com/COMBINE-lab/seqcol_rs", branch = "dev",
 anyhow = "1.0.100"
 bstr = "1.12.1"
 bio-types = { version = "1.0.4", features = ["serde"] }
-clap = { version = "4.5.54", features = ["derive"] }
+clap = { version = "4.5.57", features = ["derive"] }
 noodles-bam = "0.85.0"
 noodles-sam = "0.81.0"
 num-format = "0.4.4"
@@ -66,8 +66,8 @@ crossbeam = { version = "0.8.4", features = [
 sprs = "0.11.4"
 
 # alternative sources for dev
-minimap2-sys = { version = "0.1.24" } #, git = "https://github.com/rob-p/minimap2-rs.git", branch = "main" }
-minimap2 = { version = "0.1.28" }     # git = "https://github.com/rob-p/minimap2-rs.git", branch = "main" }
+minimap2-sys = { version = "0.1.30" } #, git = "https://github.com/rob-p/minimap2-rs.git", branch = "main" }
+minimap2 = { version = "0.1.30" }     # git = "https://github.com/rob-p/minimap2-rs.git", branch = "main" }
 
 needletail = "0.6.3"
 indicatif = "0.18.3"

src/bulk.rs

@@ -200,6 +200,7 @@ fn perform_inference_and_write_output(
             &counts,
             name_vec,
             txps_name,
+            args.display_thresh,
         )?;
     }
 

src/prog_opts.rs

@@ -43,6 +43,20 @@ fn parse_assign_prob_out_value(s: &str) -> anyhow::Result<ReadAssignmentProbOut>
     }
 }
 
+fn parse_display_thresh(s: &str) -> anyhow::Result<f64> {
+    match s.to_lowercase().as_str() {
+        "none" => Ok(f64::MIN_POSITIVE),
+        _ => {
+            let val = s.parse::<f64>()?;
+            if val >= 0.0 && val <= 1.0 {
+                Ok(val)
+            } else {
+                anyhow::bail!("display-thresh must be between 0.0 and 1.0, got {}", val)
+            }
+        }
+    }
+}
+
 #[derive(Debug, Clone, clap::ValueEnum, Serialize)]
 pub enum SequencingTech {
     OntCDNA,
@@ -383,6 +397,17 @@ pub struct Args {
     )]
     pub write_assignment_probs: Option<ReadAssignmentProbOut>,
 
+    /// minimum posterior probability threshold for a read-transcript assignment to be written
+    /// to the .prob file. Accepts a float value between 0.0 and 1.0, or the sentinel value
+    /// 'none' to use the minimum machine precision (f64::MIN_POSITIVE).
+    #[arg(
+        long,
+        help_heading = "output read-txps probabilities",
+        default_value_t = 1e-6,
+        value_parser = parse_display_thresh
+    )]
+    pub display_thresh: f64,
+
     /// maximum number of iterations for which to run the EM algorithm
     #[arg(long, help_heading = "EM", default_value_t = 1000)]
     pub max_em_iter: u32,

src/util/write_function.rs

@@ -214,6 +214,7 @@ pub fn write_out_prob(
     counts: &[f64],
     names_vec: SwapVec<String>,
     txps_name: &[String],
+    display_thresh: f64,
 ) -> anyhow::Result<()> {
     if let Some(p) = output.parent() {
         // unless this was a relative path with one component,
@@ -278,15 +279,15 @@ pub fn write_out_prob(
         txps.clear();
         txp_probs.clear();
 
-        const DISPLAY_THRESH: f64 = 0.001;
         let mut denom2 = 0.0_f64;
 
         for (a, p, cp) in izip!(alns, probs, coverage_probs) {
             let target_id = a.ref_id as usize;
             let prob = *p as f64;
             let cov_prob = if model_coverage { *cp } else { 1.0 };
             let nprob = ((counts[target_id] * prob * cov_prob) / denom).clamp(0.0, 1.0);
-            if nprob >= DISPLAY_THRESH {
+
+            if nprob >= display_thresh {
                 txps.push(target_id);
                 txp_probs.push(nprob);
                 denom2 += nprob;