From d5a1d953b80bed867c686cd01aa4d237a7dc8439 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 3 Nov 2026 10:41:36 +0100 Subject: [PATCH] CellReads.stats rows in whitelist order, not thread order The cell rows were written in reverse order of first appearance, taken from whichever thread folded a read first, so the file changed with --runThreadN. Write them in whitelist order (sorted by barcode) from the BTreeMap that already holds them. STAR's own order is a libstdc++ unordered_map walk; reproducing it would mean porting libstdc++'s hash table, so DIVERGENCE.md 3.2 documents the difference for sign-off. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 7 ++++ DIVERGENCE.md | 12 ++++--- src/solo/cell_reads.rs | 72 +++++++++++++++++------------------------- 3 files changed, 43 insertions(+), 48 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1ab8f0e6..7a164e1b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -232,6 +232,13 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Bug fixes +- **`CellReads.stats` no longer changes with `--runThreadN`.** The cell rows + were written in reverse order of first appearance, taken from whichever + thread finished a read first, so the same input gave the same rows in a + different order at 1, 4 and 8 threads. They are now written in whitelist + order (sorted by barcode). STAR's order is a libstdc++ hash-table walk that + is not reproduced; DIVERGENCE.md 3.2 says why. + - **STARsolo's stats files as STAR writes them.** `Solo.out/Barcodes.stats`, `/Features.stats` and `/UMIperCellSorted.txt` were not written, SJ and Velocyto had no `Summary.csv`, and the Gene/GeneFull diff --git a/DIVERGENCE.md b/DIVERGENCE.md index f0890fbd..43bfc110 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -225,15 +225,17 @@ There are none at present. The four yeast PE reads once listed here now match ST ### 3.2 `CellReads.stats` row order -**What STAR does.** `--soloCellReadStats CB` emits its rows by iterating a libc++ `std::unordered_map`, so the order is a hash-table walk rather than a sort. At the map sizes this produces, libc++ chains new entries at the head of their bucket and walks buckets in order, which comes out as the reverse of each barcode's first appearance in read order. +**Needs maintainer sign-off.** -**What rustar-aligner does.** Emits that same reverse-first-appearance order, including across threads: the per-read accumulator merges in read order, so a threaded run writes the same file as a serial one. +**What STAR does.** `--soloCellReadStats Standard` writes the cell rows by walking a `std::unordered_map` (`SoloCommon.h:47`, `SoloFeature_statsOutput.cpp:102`). The map is reserved for `nCB*3/2` entries (`SoloFeature_countCBgeneUMI.cpp:42`) and filled from each thread's temporary file in thread order (`:45-48`, `SoloReadFeature_inputRecords.cpp:141-152`). So the order is the standard library's bucket layout after its rehashes, not a sort. The Linux STAR build uses libstdc++; feeding its `unordered_map` the same keys in read order, reserved the same way, gives STAR's order. -**Why.** Reproducing the order where it is reproducible costs nothing and keeps a byte-comparison against STAR usable on the sizes where it can work at all. +**What rustar-aligner does.** Writes the cell rows in whitelist order, which is sorted by barcode. The `CBnotInPasslist` row stays first, as in STAR. -**Impact.** Past libc++'s load factor the map rehashes, and the order then depends on the bucket count, which depends on how many distinct barcodes were seen; beyond that size the order diverges. The **values never do** — only which line they appear on. Reading the file by barcode rather than by position is unaffected either way. +**Why.** Before this, rustar-aligner wrote the rows in reverse order of first appearance, taken from the order threads happened to finish reads, so the file changed with `--runThreadN`. That order did not match STAR's either. Matching STAR would mean re-implementing libstdc++'s hash table (its prime bucket table, rehash policy and node chaining). libstdc++ is GPL with a runtime exception, a sensitive licence under CONTRIBUTING.md, so that is not done without a maintainer decision. -**Source.** `src/solo/cell_reads.rs`, locked by `rows_are_emitted_in_reverse_first_appearance_order`. STAR: `SoloFeature_statsOutput.cpp`. +**Impact.** The **values** of every row are unchanged; only which line a barcode is on differs from STAR. A reader that looks rows up by barcode is unaffected. A line-by-line diff against STAR's file needs both sides sorted first. + +**Source.** `src/solo/cell_reads.rs`, locked by `rows_follow_whitelist_order_not_arrival_order`. STAR: `SoloFeature_statsOutput.cpp:102`. --- diff --git a/src/solo/cell_reads.rs b/src/solo/cell_reads.rs index 5d140a51..20d33e21 100644 --- a/src/solo/cell_reads.rs +++ b/src/solo/cell_reads.rs @@ -9,25 +9,15 @@ //! `CBnotInPasslist` row, so the columns account for every read rather than //! only the ones that succeeded. //! -//! # D24: row order +//! # Row order (DIVERGENCE.md 3.2) //! -//! STAR iterates a libc++ `std::unordered_map` to emit these rows, so the order -//! is a hash-table walk, not a sort. For the small maps this produces, libc++ -//! chains new entries at the head of their bucket and walks buckets in order, -//! which comes out as the reverse of first appearance in read order. That is -//! what this reproduces. -//! -//! It is not reproducible in general: past the load factor libc++ rehashes, and -//! the order after a rehash depends on the bucket count, which depends on how -//! many distinct barcodes were seen. At that size the order diverges. The -//! **values never do** — only which line they appear on. A consumer that reads -//! this file by barcode rather than by position is unaffected either way, and -//! sorting by barcode is the only stable thing to do with it. -//! -//! Reads are folded in under a mutex held by `SoloContext`, so the order is the -//! order reads were processed in regardless of thread count. +//! STAR walks a `std::unordered_map` to emit these rows, so its order is the +//! standard library's hash-table layout, not a sort. rustar-aligner writes the +//! cell rows in whitelist order (sorted by barcode) instead. That does not +//! depend on the thread count or on which thread saw a barcode first. The +//! **values** are the same either way; only which line a barcode is on differs. -use std::collections::{BTreeMap, HashSet}; +use std::collections::BTreeMap; use std::fmt::Write as _; /// What happened to one read, as the fourteen optional flags STAR tracks. @@ -73,13 +63,10 @@ pub struct CellReadFlag { /// The accumulator behind `CellReads.stats`. #[derive(Debug, Default, Clone)] pub struct CellReadStats { - /// Whitelist index to its fifteen counters. + /// Whitelist index to its fifteen counters; also the output row order. cells: BTreeMap, /// The single bucket for reads whose barcode did not resolve. no_cb: [u64; 15], - /// Whitelist indices in first-appearance order; emitted reversed (D24). - order: Vec, - seen: HashSet, } impl CellReadStats { @@ -116,9 +103,6 @@ impl CellReadStats { /// Record a read that resolved to whitelist cell `cb`. pub fn add_cell(&mut self, cb: u32, flag: &CellReadFlag) { - if self.seen.insert(cb) { - self.order.push(cb); - } Self::fold(self.cells.entry(cb).or_insert([0; 15]), flag); } @@ -146,10 +130,10 @@ impl CellReadStats { s.push_str(&v.to_string()); } s.push_str("\t0\t0\t0\t0\n"); - // Reverse first-appearance order — see D24 in the module docs. - for &cb in self.order.iter().rev() { + // Whitelist order, whatever order the threads folded reads in. + for (&cb, counts) in &self.cells { s.push_str(&barcode_of(cb)); - for v in &self.cells[&cb] { + for v in counts { s.push('\t'); s.push_str(&v.to_string()); } @@ -245,22 +229,25 @@ mod tests { assert_eq!(&row[16..], ["0", "0", "0", "0"], "no UMI columns for it"); } - /// D24: rows come out in reverse first-appearance order, which is what - /// STAR's libc++ hash-map walk produces at these sizes. + /// Rows come out in whitelist order whatever order the reads arrived in, + /// so a threaded run writes the same file as a serial one. #[test] - fn rows_are_emitted_in_reverse_first_appearance_order() { - let mut st = CellReadStats::new(); - for cb in [4u32, 1, 9] { - st.add_cell(cb, &flag_perfect_counted()); - } - st.add_cell(1, &flag_perfect_counted()); // seen again: order unchanged - let out = st.render(|i| format!("CB{i}"), &BTreeMap::new()); + fn rows_follow_whitelist_order_not_arrival_order() { + let render = |arrival: &[u32]| { + let mut st = CellReadStats::new(); + for &cb in arrival { + st.add_cell(cb, &flag_perfect_counted()); + } + st.render(|i| format!("CB{i}"), &BTreeMap::new()) + }; + let out = render(&[4, 1, 9, 1]); let cbs: Vec<&str> = out .lines() .skip(2) .map(|l| l.split('\t').next().unwrap()) .collect(); - assert_eq!(cbs, ["CB9", "CB1", "CB4"]); + assert_eq!(cbs, ["CB1", "CB4", "CB9"]); + assert_eq!(out, render(&[9, 1, 4, 1])); } /// The UMI and gene totals come from the final matrix, not from the read @@ -269,8 +256,8 @@ mod tests { #[test] fn umi_and_gene_totals_come_from_the_matrix() { let mut st = CellReadStats::new(); - st.add_cell(2, &flag_perfect_counted()); st.add_cell(3, &flag_perfect_counted()); + st.add_cell(2, &flag_perfect_counted()); let mut umi_gene = BTreeMap::new(); umi_gene.insert(2u32, (17u32, 5u32)); let out = st.render(|i| format!("CB{i}"), &umi_gene); @@ -279,10 +266,9 @@ mod tests { .skip(2) .map(|l| l.split('\t').collect()) .collect(); - // Reverse order: CB3 first, then CB2. - assert_eq!(rows[0][0], "CB3"); - assert_eq!(&rows[0][16..18], ["0", "0"]); - assert_eq!(rows[1][0], "CB2"); - assert_eq!(&rows[1][16..18], ["17", "5"]); + assert_eq!(rows[0][0], "CB2"); + assert_eq!(&rows[0][16..18], ["17", "5"]); + assert_eq!(rows[1][0], "CB3"); + assert_eq!(&rows[1][16..18], ["0", "0"]); } }