test(determinism): lock output against map iteration order (#210)

CHANGELOG.md

@@ -20,7 +20,14 @@ Sections commonly used: Features, Bug fixes, Other changes.
   implementation built the SA in 172.953 s versus 267.592 s for its original
   0.7 baseline: 35.4% faster, with peak RSS reduced from 10,512,408 to
   9,169,892 KiB. The complete output hash was unchanged.
-
+- The splice-junction sorts that produce `SJ.out.tab`, the `SJ` solo-feature
+  rows and the `BySJout` survivor set now order on the whole key (chromosome,
+  start, end, strand, motif) rather than on coordinates alone. The counts come
+  from a `DashMap`, whose iteration order varies with hashing and with
+  concurrent insertion, so a tie left to that order would have been a file
+  that differs between runs or thread counts. `tests/determinism.rs` locks it:
+  the same reads at one and at eight threads, and two runs at eight threads,
+  produce byte-identical output in single-pass and two-pass mode. Answers #210.
 - `cluster_seeds` reuses its window-bin map across reads on a thread instead
   of rebuilding it per read. Merging two windows re-keys every bin in the
   merged span, so the per-read pre-sizing was only a floor and the map

src/junction/sj_output.rs

@@ -165,12 +165,17 @@ impl SpliceJunctionStats {
             })
             .collect();
 
-        // Sort by chromosome, start, end (for distance calculation)
+        // Sort by chromosome, start, end (for distance calculation). Strand and
+        // motif join the key because the source is a `DashMap`, whose iteration
+        // order is not stable across runs or thread counts: a tie left to that
+        // order would carry it forward (#210).
         junctions.sort_by(|a, b| {
             a.0.chr_idx
                 .cmp(&b.0.chr_idx)
                 .then(a.0.intron_start.cmp(&b.0.intron_start))
                 .then(a.0.intron_end.cmp(&b.0.intron_end))
+                .then(a.0.strand.cmp(&b.0.strand))
+                .then(a.0.motif.cmp(&b.0.motif))
         });
 
         let overhang_min = &params.out_sj_filter_overhang_min;
@@ -274,17 +279,18 @@ impl SpliceJunctionStats {
     /// keys so the SJ recorder can be mapped to matrix rows.
     pub(crate) fn sj_feature_order(&self, params: &Parameters) -> Vec<(u64, u64)> {
         let surviving = self.compute_surviving_junctions(params);
-        let mut keys: Vec<(usize, u64, u64)> = self
+        let mut keys: Vec<(usize, u64, u64, u8, u8)> = self
             .junctions
             .iter()
             .filter(|e| surviving.contains(e.key()))
             .map(|e| {
                 let k = e.key();
-                (k.chr_idx, k.intron_start, k.intron_end)
+                (k.chr_idx, k.intron_start, k.intron_end, k.strand, k.motif)
             })
             .collect();
+        // Total key again: these are the SJ matrix row positions (#210).
         keys.sort_unstable();
-        keys.into_iter().map(|(_, s, e)| (s, e)).collect()
+        keys.into_iter().map(|(_, s, e, _, _)| (s, e)).collect()
     }
 
     /// Write the 9-column `SJ.out.tab` lines (sorted) to `writer`; returns the
@@ -315,11 +321,16 @@ impl SpliceJunctionStats {
             })
             .collect();
 
+        // Total order, for the same reason as `compute_surviving_junctions`:
+        // these rows are the bytes of SJ.out.tab, so a tie broken by `DashMap`
+        // iteration order would be a file that differs between runs (#210).
         output_junctions.sort_by(|a, b| {
             a.0.chr_idx
                 .cmp(&b.0.chr_idx)
                 .then(a.0.intron_start.cmp(&b.0.intron_start))
                 .then(a.0.intron_end.cmp(&b.0.intron_end))
+                .then(a.0.strand.cmp(&b.0.strand))
+                .then(a.0.motif.cmp(&b.0.motif))
         });
 
         let mut written = 0u32;

src/ruSTAR.code-workspace

@@ -0,0 +1,7 @@
+{
+    "folders": [
+        {
+            "path": ".."
+        }
+    ]
+}

tests/determinism.rs

@@ -0,0 +1,276 @@
+//! Output must not depend on how many threads produced it.
+//!
+//! Junction counts live in a `DashMap`, whose iteration order varies with
+//! hashing and with concurrent insertion. Every path that emits an order sorts
+//! first, and these tests are what keeps that true: they align the same reads
+//! at one thread and at eight and compare the output files byte for byte
+//! (issue #210).
+
+use assert_cmd::cargo::cargo_bin_cmd;
+use std::fs;
+use std::io::Write;
+use std::path::{Path, PathBuf};
+use tempfile::TempDir;
+
+/// Two exons with a GT-AG intron between them, so reads spanning the junction
+/// produce SJ.out.tab rows rather than an empty file.
+const EXON1: (usize, usize) = (2_000, 2_300);
+const INTRON: (usize, usize) = (2_300, 2_800);
+const EXON2: (usize, usize) = (2_800, 3_100);
+/// A second junction close to the first, to give the neighbour-distance filter
+/// something to compute and the sort something to order.
+const EXON3: (usize, usize) = (6_000, 6_200);
+const INTRON2: (usize, usize) = (6_200, 6_600);
+const EXON4: (usize, usize) = (6_600, 6_900);
+
+fn lcg_seq(seed: u32, length: usize) -> Vec<u8> {
+    let bases: [u8; 4] = *b"ACGT";
+    let mut state = seed;
+    let mut seq = Vec::with_capacity(length);
+    for _ in 0..length {
+        state = state.wrapping_mul(1_103_515_245).wrapping_add(12345);
+        seq.push(bases[((state >> 16) & 3) as usize]);
+    }
+    seq
+}
+
+fn build_genome() -> Vec<u8> {
+    let mut genome = lcg_seq(88888, 20_000);
+    for (start, end) in [INTRON, INTRON2] {
+        genome[start] = b'G';
+        genome[start + 1] = b'T';
+        genome[end - 2] = b'A';
+        genome[end - 1] = b'G';
+    }
+    genome
+}
+
+fn write_fasta(dir: &Path, genome: &[u8]) -> PathBuf {
+    let path = dir.join("genome.fa");
+    let mut f = fs::File::create(&path).unwrap();
+    writeln!(f, ">chr1").unwrap();
+    f.write_all(genome).unwrap();
+    writeln!(f).unwrap();
+    path
+}
+
+fn build_index(fasta: &Path, genome_dir: &Path) {
+    fs::create_dir_all(genome_dir).unwrap();
+    cargo_bin_cmd!("rustar-aligner")
+        .args([
+            "--runMode",
+            "genomeGenerate",
+            "--genomeDir",
+            genome_dir.to_str().unwrap(),
+            "--genomeFastaFiles",
+            fasta.to_str().unwrap(),
+            "--genomeSAindexNbases",
+            "7",
+        ])
+        .assert()
+        .success();
+}
+
+/// Reads spanning both junctions plus unspliced filler, enough of them that
+/// several threads each hold some.
+fn write_fastq(dir: &Path, genome: &[u8]) -> PathBuf {
+    let path = dir.join("reads.fq");
+    let mut f = fs::File::create(&path).unwrap();
+    let mut n = 0usize;
+
+    let spliced = |f: &mut fs::File, e_end: usize, i_end: usize, overhang: usize, n: &mut usize| {
+        let left = &genome[e_end - overhang..e_end];
+        let right = &genome[i_end..i_end + (100 - overhang)];
+        let mut seq = left.to_vec();
+        seq.extend_from_slice(right);
+        writeln!(f, "@sj{n}").unwrap();
+        f.write_all(&seq).unwrap();
+        writeln!(f, "\n+\n{}", "I".repeat(seq.len())).unwrap();
+        *n += 1;
+    };
+
+    for overhang in 30..70 {
+        spliced(&mut f, INTRON.0, INTRON.1, overhang, &mut n);
+        spliced(&mut f, INTRON2.0, INTRON2.1, overhang, &mut n);
+    }
+    for i in 0..200usize {
+        let start = 8_000 + i * 40;
+        writeln!(f, "@u{i}").unwrap();
+        f.write_all(&genome[start..start + 100]).unwrap();
+        writeln!(f, "\n+\n{}", "I".repeat(100)).unwrap();
+    }
+    let _ = (EXON1, EXON2, EXON3, EXON4);
+    path
+}
+
+/// Align with `threads` threads and return the bytes of every output file that
+/// carries an order.
+fn align(
+    genome_dir: &Path,
+    fastq: &Path,
+    prefix: &str,
+    threads: &str,
+    extra: &[&str],
+) -> Vec<(String, Vec<u8>)> {
+    let mut cmd = cargo_bin_cmd!("rustar-aligner");
+    cmd.args([
+        "--runMode",
+        "alignReads",
+        "--genomeDir",
+        genome_dir.to_str().unwrap(),
+        "--readFilesIn",
+        fastq.to_str().unwrap(),
+        "--runThreadN",
+        threads,
+        "--outFileNamePrefix",
+        prefix,
+    ]);
+    cmd.args(extra);
+    cmd.assert().success();
+
+    ["SJ.out.tab", "Aligned.out.sam"]
+        .iter()
+        .map(|name| {
+            let bytes = fs::read(format!("{prefix}{name}")).unwrap_or_default();
+            // Drop the @PG header line: it records the command line, which
+            // differs by the thread count itself.
+            let filtered: Vec<u8> = if *name == "Aligned.out.sam" {
+                String::from_utf8_lossy(&bytes)
+                    .lines()
+                    .filter(|l| !l.starts_with("@PG"))
+                    .collect::<Vec<_>>()
+                    .join("\n")
+                    .into_bytes()
+            } else {
+                bytes
+            };
+            ((*name).to_string(), filtered)
+        })
+        .collect()
+}
+
+#[test]
+fn output_is_identical_at_one_and_eight_threads() {
+    let dir = TempDir::new().unwrap();
+    let root = dir.path();
+    let genome = build_genome();
+    let genome_dir = root.join("genome");
+    build_index(&write_fasta(root, &genome), &genome_dir);
+    let fq = write_fastq(root, &genome);
+
+    let one = align(
+        &genome_dir,
+        &fq,
+        &format!("{}/t1_", root.display()),
+        "1",
+        &[],
+    );
+    let eight = align(
+        &genome_dir,
+        &fq,
+        &format!("{}/t8_", root.display()),
+        "8",
+        &[],
+    );
+
+    // The fixture has to actually produce junctions, or this test would pass
+    // on two empty files.
+    let sj = &one[0].1;
+    let rows = String::from_utf8_lossy(sj).lines().count();
+    println!("fixture produced {rows} SJ.out.tab rows");
+    assert!(
+        rows >= 2,
+        "the fixture has to produce at least two junction rows for an order to \
+         exist at all, got {rows}"
+    );
+
+    for ((name, a), (_, b)) in one.iter().zip(eight.iter()) {
+        assert_eq!(
+            String::from_utf8_lossy(a),
+            String::from_utf8_lossy(b),
+            "{name} differs between 1 and 8 threads"
+        );
+    }
+}
+
+#[test]
+fn output_is_identical_across_two_runs_at_eight_threads() {
+    // Guards against the pair above agreeing only because both runs happened
+    // to hit the same map order.
+    let dir = TempDir::new().unwrap();
+    let root = dir.path();
+    let genome = build_genome();
+    let genome_dir = root.join("genome");
+    build_index(&write_fasta(root, &genome), &genome_dir);
+    let fq = write_fastq(root, &genome);
+
+    let a = align(
+        &genome_dir,
+        &fq,
+        &format!("{}/a_", root.display()),
+        "8",
+        &[],
+    );
+    let b = align(
+        &genome_dir,
+        &fq,
+        &format!("{}/b_", root.display()),
+        "8",
+        &[],
+    );
+
+    for ((name, x), (_, y)) in a.iter().zip(b.iter()) {
+        assert_eq!(
+            String::from_utf8_lossy(x),
+            String::from_utf8_lossy(y),
+            "{name} differs between two 8-thread runs"
+        );
+    }
+}
+
+#[test]
+fn two_pass_output_is_identical_at_one_and_eight_threads() {
+    // Two-pass feeds pass 1's junctions back into the alignment, so any order
+    // escaping the junction map would show up here rather than in a single
+    // pass.
+    let dir = TempDir::new().unwrap();
+    let root = dir.path();
+    let genome = build_genome();
+    let genome_dir = root.join("genome");
+    build_index(&write_fasta(root, &genome), &genome_dir);
+    let fq = write_fastq(root, &genome);
+
+    let extra = ["--twopassMode", "Basic"];
+    let one = align(
+        &genome_dir,
+        &fq,
+        &format!("{}/tp1_", root.display()),
+        "1",
+        &extra,
+    );
+    let eight = align(
+        &genome_dir,
+        &fq,
+        &format!("{}/tp8_", root.display()),
+        "8",
+        &extra,
+    );
+
+    for ((name, a), (_, b)) in one.iter().zip(eight.iter()) {
+        assert_eq!(
+            String::from_utf8_lossy(a),
+            String::from_utf8_lossy(b),
+            "two-pass {name} differs between 1 and 8 threads"
+        );
+    }
+
+    // The pass-1 junction file is written from the same map and is equally
+    // order-sensitive.
+    let p1 = fs::read(format!("{}/tp1__STARpass1/SJ.out.tab", root.display())).unwrap();
+    let p8 = fs::read(format!("{}/tp8__STARpass1/SJ.out.tab", root.display())).unwrap();
+    assert_eq!(
+        String::from_utf8_lossy(&p1),
+        String::from_utf8_lossy(&p8),
+        "pass-1 SJ.out.tab differs between 1 and 8 threads"
+    );
+}