fix: make dupRadar and Qualimap outputs reproducible run-to-run

src/rna/dupradar/plots.rs

@@ -15,6 +15,12 @@ use plotters::prelude::*;
 use plotters_svg::SVGBackend;
 use std::collections::HashMap;
 
+/// One deduplicated scatter point: `(pixel coordinate, (data x, data y, density))`.
+///
+/// The pixel coordinate is kept alongside the data so the draw order can be
+/// tie-broken deterministically when two points share a density.
+type PixelPoint = ((i32, i32), (f64, f64, f64));
+
 // ---------------------------------------------------------------------------
 // Constants
 // ---------------------------------------------------------------------------
@@ -597,13 +603,22 @@ where
             .or_insert((x, y, d));
     }
 
-    // Sort by density ascending so high-density points draw on top
-    let mut deduped: Vec<(f64, f64, f64)> = pixel_map.into_values().collect();
-    deduped.sort_by(|a, b| a.2.partial_cmp(&b.2).unwrap_or(std::cmp::Ordering::Equal));
-
-    let max_dens = deduped.iter().map(|d| d.2).fold(0.0f64, f64::max);
-
-    for (x, y, d) in deduped {
+    // Sort by density ascending so high-density points draw on top.
+    // Points of equal density would otherwise keep `pixel_map`'s iteration
+    // order (stable sort over a HashMap), which varies per process and makes
+    // the SVG/PNG differ between runs on identical input. Tie-break on the
+    // pixel coordinate so the draw order is reproducible.
+    let mut deduped: Vec<PixelPoint> = pixel_map.into_iter().collect();
+    deduped.sort_by(|a, b| {
+        a.1 .2
+            .partial_cmp(&b.1 .2)
+            .unwrap_or(std::cmp::Ordering::Equal)
+            .then_with(|| a.0.cmp(&b.0))
+    });
+
+    let max_dens = deduped.iter().map(|d| d.1 .2).fold(0.0f64, f64::max);
+
+    for (_, (x, y, d)) in deduped {
         let t = if max_dens > 0.0 { d / max_dens } else { 0.0 };
         let c = density_color(t);
         chart.draw_series(std::iter::once(Circle::new(

src/rna/qualimap/output.rs

@@ -105,8 +105,10 @@ struct TranscriptCoverageEntry {
     strand: char,
     /// Gene index for best-per-gene selection.
     gene_idx: u32,
-    /// Flat transcript index (for diagnostics).
-    #[allow(dead_code)]
+    /// Flat transcript index, assigned in GTF order.
+    ///
+    /// Used as the deterministic tie-breaker when ranking transcripts by mean
+    /// coverage in [`compute_bias`].
     flat_idx: u32,
 }
 
@@ -202,8 +204,17 @@ fn compute_bias(entries: &[TranscriptCoverageEntry]) -> (f64, f64, f64) {
         return (f64::NAN, f64::NAN, f64::NAN);
     }
 
-    // Sort by mean coverage descending, take top N
-    qualifying.sort_by(|a, b| b.0.partial_cmp(&a.0).unwrap_or(std::cmp::Ordering::Equal));
+    // Sort by mean coverage descending, take top N.
+    // `best_per_gene` is a HashMap and `sort_by` is stable, so without a
+    // tie-breaker the transcripts that survive `truncate` at the 1000-entry
+    // boundary would depend on HashMap iteration order — and with them the
+    // reported bias values. Tie-break on the transcript's flat index, which is
+    // assigned in GTF order and is therefore stable across runs.
+    qualifying.sort_by(|a, b| {
+        b.0.partial_cmp(&a.0)
+            .unwrap_or(std::cmp::Ordering::Equal)
+            .then_with(|| a.1.flat_idx.cmp(&b.1.flat_idx))
+    });
     qualifying.truncate(NUM_TRANSCRIPTS_FOR_BIAS);
 
     let mut five_prime_biases = Vec::with_capacity(qualifying.len());
@@ -755,7 +766,14 @@ fn write_results_file(
             .iter()
             .map(|(motif, &count)| (motif.clone(), count as f64 * 100.0 / total_junctions))
             .collect();
-        motif_pcts.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal));
+        // `junction_motifs` is a HashMap and many 4-mers tie at the same
+        // percentage, so tie-break on the motif to keep the top-11 list (and
+        // its order) identical across runs.
+        motif_pcts.sort_by(|a, b| {
+            b.1.partial_cmp(&a.1)
+                .unwrap_or(std::cmp::Ordering::Equal)
+                .then_with(|| a.0.cmp(&b.0))
+        });
 
         // Qualimap shows top 11 motifs (count <= 10 in their loop)
         for (motif, pct) in motif_pcts.iter().take(11) {

src/rna/qualimap/report.rs

@@ -491,7 +491,10 @@ fn write_summary_section(html: &mut String, data: &ReportData) {
     // Top junction motifs sorted by count descending — Qualimap shows motif / percentage
     // Use reads_at_junctions as denominator (matches text file output and upstream Qualimap)
     let mut motifs: Vec<(&String, &u64)> = data.junction_motifs.iter().collect();
-    motifs.sort_by(|a, b| b.1.cmp(a.1));
+    // Tie-break on the motif: `junction_motifs` is a HashMap and `sort_by` is
+    // stable, so equal counts would otherwise be ordered by HashMap iteration
+    // order and the report would differ between runs.
+    motifs.sort_by(|a, b| b.1.cmp(a.1).then_with(|| a.0.cmp(b.0)));
     let total_junctions = data.reads_at_junctions;
     for (motif, &count) in motifs.iter().take(11) {
         let pct = if total_junctions > 0 {