From 046412767d640ae68dccf56d39ea099f197fd08f Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 14 Aug 2026 02:27:34 +0200 Subject: [PATCH] fix: make dupRadar and Qualimap outputs reproducible MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Fixes #141. Three outputs differed between two runs of the same binary on the same input. In each case a HashMap was collected into a Vec and sorted with a comparator that has ties; `sort_by` is stable, so tied elements kept HashMap iteration order, which `RandomState` randomises per process. - `qualimap/output.rs`: junction motifs sorted by percentage, then truncated to the top 11. Many 4-mers tie, so which ones appear and in what order varied. Tie-break on the motif. - `qualimap/report.rs`: the same motif list in the HTML report. Same fix. - `dupradar/plots.rs`: the density scatter deduplicates points per pixel into a HashMap, then sorts by density so dense points draw last. Points of equal density came out in a different sequence each run, changing the SVG element order and the PNG bytes. Keep the pixel coordinate alongside the point and tie-break on it. Also fixes the same pattern in `compute_bias()`, which had a worse consequence: the transcript list is sorted by mean coverage and truncated to the top 1000, so ties at that boundary decided which transcripts fed the 5'/3' bias calculation. Tie-break on the transcript's flat index, assigned in GTF order. `flat_idx` is now used rather than dead code. Upstream Qualimap orders these by Java HashMap iteration, i.e. arbitrarily, so a deterministic tie-break is no less faithful than the current behaviour. Verified with three consecutive runs on a 4M-read BAM: `rnaseq_qc_results.txt`, the scatter `.svg` and `.png` are now byte-identical across runs. The contents are unchanged from before — sorting the old and new files gives identical multisets, only the order of tied entries moved. `qualimapReport.html` still differs between runs by exactly two lines: the "Analysis date" timestamp and its copy in the footer. That is intentional report metadata and matches upstream Qualimap. `cargo test --release` passes (200 + 12 + 18 + 2). Co-Authored-By: Claude Opus 5 (1M context) --- src/rna/dupradar/plots.rs | 29 ++++++++++++++++++++++------- src/rna/qualimap/output.rs | 28 +++++++++++++++++++++++----- src/rna/qualimap/report.rs | 5 ++++- 3 files changed, 49 insertions(+), 13 deletions(-) diff --git a/src/rna/dupradar/plots.rs b/src/rna/dupradar/plots.rs index 802fc37d..6166fa8a 100644 --- a/src/rna/dupradar/plots.rs +++ b/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 = 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( diff --git a/src/rna/qualimap/output.rs b/src/rna/qualimap/output.rs index 5d53082e..bd980118 100644 --- a/src/rna/qualimap/output.rs +++ b/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) { diff --git a/src/rna/qualimap/report.rs b/src/rna/qualimap/report.rs index bf0b4882..2e31f859 100644 --- a/src/rna/qualimap/report.rs +++ b/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 {