Skip to content

Three outputs are not reproducible run-to-run (HashMap iteration order reaching output) #141

Description

@BenjaminDEMAILLE

Summary

Three output files are not reproducible run-to-run. Running the same binary on the same input twice produces files that do not hash-match:

  • dupradar/<sample>_duprateExpDens.svg and the matching .png
  • qualimap/qualimapReport.html
  • qualimap/rnaseq_qc_results.txt

No numeric value changes — the differences are pure reordering of equal-valued entries — but the files differ byte-for-byte, which breaks checksum-based pipeline caching (Nextflow storeDir/resume, md5sum in test assertions) and makes it impossible to assert byte-identical output in CI.

Reproduce

rustqc rna sample.bam --gtf annotation.gtf -o run_a --threads 4
rustqc rna sample.bam --gtf annotation.gtf -o run_b --threads 4
diff run_a/qualimap/rnaseq_qc_results.txt run_b/qualimap/rnaseq_qc_results.txt

Observed on main (bea5571), aarch64 macOS, on a 4M-read coordinate-sorted BAM with a 60k-gene GTF. Example diff:

55d54
<     GCCC : 0.41%
56a56
>     GCCC : 0.41%

GCCC and its neighbour both sit at 0.41%; which one is printed first varies per run.

For comparison, every other output is stable: featureCounts.tsv, dupMatrix.txt, tin.xls, lc_extrap.txt, samtools stats/flagstat/idxstats, and all RSeQC outputs are byte-identical across runs.

Cause

The same pattern in three places: a HashMap is collected into a Vec and then sorted with a comparator that has ties. slice::sort_by is stable, so tied elements keep their input order — which here is HashMap iteration order, randomised per process by RandomState.

  1. src/rna/qualimap/output.rs:757 — junction motifs sorted by percentage descending, then .take(11). junction_motifs is HashMap<String, u64> and many 4-mers tie at the same percentage.

    let mut motif_pcts: Vec<(String, f64)> = result.junction_motifs.iter()...collect();
    motif_pcts.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(Ordering::Equal));
  2. src/rna/qualimap/report.rs:494 — the same motif list rendered into the HTML report:

    let mut motifs: Vec<(&String, &u64)> = data.junction_motifs.iter().collect();
    motifs.sort_by(|a, b| b.1.cmp(a.1));
  3. src/rna/dupradar/plots.rs:601 — the density scatter deduplicates points per pixel into HashMap<(i32, i32), (f64, f64, f64)>, then sorts by density ascending so dense points draw last. Points with equal density are emitted in HashMap order, so the <circle> elements come out in a different sequence each run:

    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(Ordering::Equal));

    Since ties are only reordered among points of equal density, the rendered image is visually identical; only the file bytes move.

Related latent case

compute_bias() in src/rna/qualimap/output.rs:206 has the same shape but a worse consequence:

let mut qualifying: Vec<(f64, &TranscriptCoverageEntry)> = best_per_gene.values()...collect();
qualifying.sort_by(|a, b| b.0.partial_cmp(&a.0).unwrap_or(Ordering::Equal));
qualifying.truncate(NUM_TRANSCRIPTS_FOR_BIAS);

best_per_gene is a HashMap, and the sort is truncated to the top 1000. When transcripts tie at the 1000-entry boundary, which transcripts survive depends on HashMap iteration order, so the reported 5'/3'/5'-3' bias values themselves can change between runs. I have not caught this firing in practice (bias values were stable across my runs), but it is the same bug with numeric rather than cosmetic impact. The existing comment at output.rs:176-179 acknowledges a related tie-break decision.

Suggested fix

Add a deterministic tie-breaker to each comparator rather than swapping map types:

  • motifs: tie-break on the motif string
  • scatter points: tie-break on the pixel coordinate
  • compute_bias: tie-break on the transcript's flat index

Upstream Qualimap's own ordering here is Java HashMap order, i.e. arbitrary, so any deterministic tie-break is at least as faithful to upstream as the current behaviour while being reproducible.

Happy to send a PR.


🤖 Generated with Claude Code

Activity

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions