diff --git a/.claude/worktrees/agent-aa8bc29ee5b6ad4d5 b/.claude/worktrees/agent-aa8bc29ee5b6ad4d5 new file mode 160000 index 0000000..40d9286 --- /dev/null +++ b/.claude/worktrees/agent-aa8bc29ee5b6ad4d5 @@ -0,0 +1 @@ +Subproject commit 40d92861cf3a1c2125d68a589655c0277fa519cb diff --git a/CHANGELOG.md b/CHANGELOG.md index 7c20e6c..1a64eaa 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/src/junction/sj_output.rs b/src/junction/sj_output.rs index 76b6527..05c3082 100644 --- a/src/junction/sj_output.rs +++ b/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 = ¶ms.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; diff --git a/src/ruSTAR.code-workspace b/src/ruSTAR.code-workspace new file mode 100644 index 0000000..9e68e72 --- /dev/null +++ b/src/ruSTAR.code-workspace @@ -0,0 +1,7 @@ +{ + "folders": [ + { + "path": ".." + } + ] +} \ No newline at end of file diff --git a/test/__pycache__/bench_report.cpython-312.pyc b/test/__pycache__/bench_report.cpython-312.pyc new file mode 100644 index 0000000..d9ad800 Binary files /dev/null and b/test/__pycache__/bench_report.cpython-312.pyc differ diff --git a/test/__pycache__/nfcore_diff.cpython-312.pyc b/test/__pycache__/nfcore_diff.cpython-312.pyc new file mode 100644 index 0000000..605b024 Binary files /dev/null and b/test/__pycache__/nfcore_diff.cpython-312.pyc differ diff --git a/test/__pycache__/speed_bench.cpython-312.pyc b/test/__pycache__/speed_bench.cpython-312.pyc new file mode 100644 index 0000000..ca6d35d Binary files /dev/null and b/test/__pycache__/speed_bench.cpython-312.pyc differ diff --git a/tests/determinism.rs b/tests/determinism.rs new file mode 100644 index 0000000..a778748 --- /dev/null +++ b/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 { + 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 { + 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)> { + 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 = if *name == "Aligned.out.sam" { + String::from_utf8_lossy(&bytes) + .lines() + .filter(|l| !l.starts_with("@PG")) + .collect::>() + .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" + ); +}