From 4d244f988cea75d63acb5cd036329052402d8efa Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 09:41:56 +0200 Subject: [PATCH 1/4] feat(solo): --soloCellReadStats CB writes CellReads.stats MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit One row per cell barcode with the fifteen counters STARsolo reports: how the barcode matched, whether the read mapped to one locus or several, whether it landed on a feature, where in the gene and on which strand, whether it was mitochondrial, and whether it reached the matrix. The per-cell UMI and gene totals come from the raw matrix rather than from the read counters, so they agree with what the matrix says by construction. Reads whose barcode never resolved are summed into a single `CBnotInPasslist` row instead of being dropped. That row is the reason the file is useful: it is the difference between "these cells look thin" and "most of the input never reached a cell at all". The region columns split by strand — an antisense read counts under `exonicAS` or `intronicAS`, never under `exonic` or `intronic`. `--genomeChrSetMitochondrial` names the chromosomes behind the `mito` column. Without it the column is zero throughout, which is honest: no chromosome was declared mitochondrial. D24 comes with it. STAR emits these rows by walking a libc++ `unordered_map`, which at these sizes is the reverse of each barcode's first appearance. That is reproduced, including across threads: the per-read accumulator merges in read order, so a threaded run writes the same file as a serial one. It stops being reproducible past the point where libc++ rehashes, since the order then depends on the bucket count. The values never differ, only which line they sit on. Recorded in docs-old/dev/divergences.md. --- CHANGELOG.md | 21 +++ src/params/mod.rs | 21 +++ src/solo/cell_reads.rs | 332 +++++++++++++++++++++++++++++++++++++++++ src/solo/count.rs | 21 +++ src/solo/mod.rs | 74 +++++++++ 5 files changed, 469 insertions(+) create mode 100644 src/solo/cell_reads.rs diff --git a/CHANGELOG.md b/CHANGELOG.md index b6e4893..a17dd8f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -97,6 +97,27 @@ Sections commonly used: Features, Bug fixes, Other changes. glibc's malloc and per-thread heaps that return whole segments to the OS when abandoned, so allocator cache size stays bounded. +- **`--soloCellFilter EmptyDrops_CR` now uses CellRanger's actual + statistics.** The ambient profile is smoothed with Simple Good-Turing, + as CellRanger and STAR do, instead of an approximation that reserved + unseen mass from the singleton rate and spread the remainder in + proportion to raw counts. The Monte-Carlo null is drawn with libc++'s + `std::mt19937` and `std::discrete_distribution`, seeded + `19760110 * (isim + 1)` per simulation as STAR seeds it, replacing a + SplitMix64 stream that could not agree with STAR's over an arbitrary + number of draws. Cell calls move as a result. + +- **`--soloCellReadStats CB`** writes `Solo.out//CellReads.stats`: + one row per cell barcode with fifteen counters describing what + happened to its reads — barcode match quality, unique or multi + genomic mapping, feature assignment, exonic/intronic and their + antisense counterparts, mitochondrial, and whether the read reached + the matrix — plus the per-cell UMI and gene totals. Reads whose + barcode never resolved are summed into a `CBnotInPasslist` row rather + than dropped, so the columns account for the whole input. + `--genomeChrSetMitochondrial` names the chromosomes behind the `mito` + column. + ### Bug fixes - **STARsolo `Gene` assignment now requires exon concordance**, matching diff --git a/src/params/mod.rs b/src/params/mod.rs index 0536a85..1b7bf48 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1113,6 +1113,17 @@ pub struct Parameters { #[arg(long = "soloCBmatchWLtype", default_value = "1MM_multi")] pub solo_cb_match_wl_type: String, + /// `CB`: write `Solo.out//CellReads.stats`, a per-cell-barcode + /// summary of what happened to the reads carrying it. `None` (the default) + /// writes nothing. + #[arg(long = "soloCellReadStats", default_value = "None")] + pub solo_cell_read_stats: String, + + /// Chromosome names treated as mitochondrial, for the `mito` column of + /// `CellReads.stats`. `-` (the default) names none. + #[arg(long = "genomeChrSetMitochondrial", num_args = 1.., default_values_t = vec!["-".to_string()])] + pub genome_chr_set_mitochondrial: Vec, + /// Cell-calling / matrix filtering: None, CellRanger2.2, EmptyDrops_CR, TopCells. #[arg(long = "soloCellFilter", num_args = 1.., default_values_t = vec!["CellRanger2.2".to_string(), "3000".to_string(), "0.99".to_string(), "10".to_string()])] pub solo_cell_filter: Vec, @@ -1710,6 +1721,16 @@ impl Parameters { )); } } + // --soloCellReadStats: `CB` is the only value STAR defines. + if !matches!(params.solo_cell_read_stats.as_str(), "CB" | "None") { + return Err(command.error( + ErrorKind::InvalidValue, + format!( + "unknown --soloCellReadStats '{}'; expected CB or None", + params.solo_cell_read_stats + ), + )); + } // Validate --clipAdapterType. if !matches!( params.clip_adapter_type.as_str(), diff --git a/src/solo/cell_reads.rs b/src/solo/cell_reads.rs new file mode 100644 index 0000000..e17a540 --- /dev/null +++ b/src/solo/cell_reads.rs @@ -0,0 +1,332 @@ +//! `--soloCellReadStats CB`: the per-cell-barcode read summary STARsolo writes +//! as `Solo.out//CellReads.stats`. +//! +//! One row per cell barcode, fifteen counters describing what happened to the +//! reads carrying it — how the barcode matched, whether the read mapped +//! uniquely, whether it landed on a feature, where in the gene, and whether it +//! reached the matrix — plus the per-cell UMI and gene totals. The reads whose +//! barcode never resolved are not dropped; they are summed into a single +//! `CBnotInPasslist` row, so the columns account for every read rather than +//! only the ones that succeeded. +//! +//! # D24: row order +//! +//! 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. + +use std::collections::{BTreeMap, HashSet}; +use std::fmt::Write as _; + +/// What happened to one read, as the fourteen optional flags STAR tracks. +/// +/// `cbMatch` is not here: every read that reaches the accumulator matched a +/// barcode well enough to be attributed somewhere, so it is always set. +/// +/// Fourteen bools rather than a bitfield because they are written once per read +/// and read once per fold, and the names are the column names of the file. +#[allow(clippy::struct_excessive_bools)] +#[derive(Debug, Clone, Copy, Default)] +pub struct CellReadFlag { + /// The barcode was an exact whitelist hit. + pub cb_perfect: bool, + /// Corrected via a single one-mismatch neighbour. + pub cb_mm_unique: bool, + /// Corrected via several one-mismatch neighbours, resolved by the posterior. + pub cb_mm_multiple: bool, + /// The read mapped to exactly one genomic locus. + pub genome_u: bool, + /// The read mapped to several genomic loci. + pub genome_m: bool, + /// It fell on exactly one feature (gene). + pub feature_u: bool, + /// It fell on several features. + pub feature_m: bool, + /// Exonic, on the annotated strand. + pub exonic: bool, + /// Intronic, on the annotated strand. + pub intronic: bool, + /// Exonic, antisense to the annotation. + pub exonic_as: bool, + /// Intronic, antisense to the annotation. + pub intronic_as: bool, + /// On a chromosome named by `--genomeChrSetMitochondrial`. + pub mito: bool, + /// Counted into the unique-gene matrix. + pub counted_u: bool, + /// Counted through the multi-gene distribution. + pub counted_m: bool, +} + +/// The accumulator behind `CellReads.stats`. +#[derive(Debug, Default, Clone)] +pub struct CellReadStats { + /// Whitelist index to its fifteen counters. + 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 { + pub fn new() -> Self { + Self::default() + } + + fn fold(v: &mut [u64; 15], f: &CellReadFlag) { + v[0] += 1; // cbMatch: every read that got this far + for (i, set) in [ + f.cb_perfect, + f.cb_mm_unique, + f.cb_mm_multiple, + f.genome_u, + f.genome_m, + f.feature_u, + f.feature_m, + f.exonic, + f.intronic, + f.exonic_as, + f.intronic_as, + f.mito, + f.counted_u, + f.counted_m, + ] + .into_iter() + .enumerate() + { + if set { + v[i + 1] += 1; + } + } + } + + /// 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); + } + + /// Record a read whose barcode did not resolve, or whose UMI was rejected. + pub fn add_no_cb(&mut self, flag: &CellReadFlag) { + Self::fold(&mut self.no_cb, flag); + } + + /// Merge another accumulator in, preserving first-appearance order: the + /// per-thread partials are merged in read order, so the combined order is + /// the order a single-threaded run would have produced. + pub fn merge(&mut self, other: &Self) { + for &cb in &other.order { + if self.seen.insert(cb) { + self.order.push(cb); + } + let dst = self.cells.entry(cb).or_insert([0; 15]); + for (d, s) in dst.iter_mut().zip(&other.cells[&cb]) { + *d += s; + } + } + for (d, s) in self.no_cb.iter_mut().zip(&other.no_cb) { + *d += s; + } + } + + /// Render the file. `umi_gene` gives each cell its final + /// `(nUMIunique, nGenesUnique)`; a cell absent from it prints zeros. + /// `barcode_of` renders a whitelist index as its barcode string. + pub fn render( + &self, + barcode_of: impl Fn(u32) -> String, + umi_gene: &BTreeMap, + ) -> String { + let mut s = String::from( + "CB\tcbMatch\tcbPerfect\tcbMMunique\tcbMMmultiple\tgenomeU\tgenomeM\tfeatureU\t\ + featureM\texonic\tintronic\texonicAS\tintronicAS\tmito\tcountedU\tcountedM\t\ + nUMIunique\tnGenesUnique\tnUMImulti\tnGenesMulti\n", + ); + s.push_str("CBnotInPasslist"); + for v in &self.no_cb { + s.push('\t'); + 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() { + s.push_str(&barcode_of(cb)); + for v in &self.cells[&cb] { + s.push('\t'); + s.push_str(&v.to_string()); + } + let (n_umi, n_gene) = umi_gene.get(&cb).copied().unwrap_or((0, 0)); + // nUMImulti and nGenesMulti stay zero: multi-gene UMIs are not + // collapsed into per-cell totals here. + let _ = writeln!(s, "\t{n_umi}\t{n_gene}\t0\t0"); + } + s + } +} + +#[cfg(test)] +mod tests { + use super::*; + + fn flag_perfect_counted() -> CellReadFlag { + CellReadFlag { + cb_perfect: true, + genome_u: true, + feature_u: true, + exonic: true, + counted_u: true, + ..Default::default() + } + } + + fn header_and_rows(text: &str) -> (Vec<&str>, Vec<&str>) { + let mut lines = text.lines(); + let header: Vec<&str> = lines.next().unwrap().split('\t').collect(); + (header, lines.collect()) + } + + /// Every column has a name, and every row has a value under each of them. + #[test] + fn each_row_fills_the_header() { + let mut st = CellReadStats::new(); + st.add_cell(0, &flag_perfect_counted()); + st.add_no_cb(&CellReadFlag::default()); + let out = st.render(|i| format!("CB{i}"), &BTreeMap::new()); + let (header, rows) = header_and_rows(&out); + assert_eq!(header.len(), 20); + assert_eq!(rows.len(), 2, "the passlist-miss row plus one cell"); + for row in rows { + assert_eq!(row.split('\t').count(), header.len()); + } + } + + /// A read is counted once under `cbMatch` and once under each flag it sets, + /// so the flag columns can exceed neither `cbMatch` nor each other's logic. + #[test] + fn flags_accumulate_per_read() { + let mut st = CellReadStats::new(); + for _ in 0..3 { + st.add_cell(7, &flag_perfect_counted()); + } + st.add_cell( + 7, + &CellReadFlag { + cb_mm_unique: true, + genome_m: true, + ..Default::default() + }, + ); + let out = st.render(|i| format!("CB{i}"), &BTreeMap::new()); + let row: Vec<&str> = out.lines().nth(2).unwrap().split('\t').collect(); + assert_eq!(row[0], "CB7"); + assert_eq!(row[1], "4", "cbMatch counts every read"); + assert_eq!(row[2], "3", "cbPerfect"); + assert_eq!(row[3], "1", "cbMMunique"); + assert_eq!(row[5], "3", "genomeU"); + assert_eq!(row[6], "1", "genomeM"); + assert_eq!(row[14], "3", "countedU"); + } + + /// Reads whose barcode never resolved are summed rather than dropped, so + /// the file accounts for the whole input. + #[test] + fn unresolved_reads_land_in_the_passlist_miss_row() { + let mut st = CellReadStats::new(); + st.add_cell(0, &flag_perfect_counted()); + for _ in 0..5 { + st.add_no_cb(&CellReadFlag { + genome_u: true, + ..Default::default() + }); + } + let out = st.render(|i| format!("CB{i}"), &BTreeMap::new()); + let row: Vec<&str> = out.lines().nth(1).unwrap().split('\t').collect(); + assert_eq!(row[0], "CBnotInPasslist"); + assert_eq!(row[1], "5"); + assert_eq!(row[5], "5", "genomeU"); + 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. + #[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()); + let cbs: Vec<&str> = out + .lines() + .skip(2) + .map(|l| l.split('\t').next().unwrap()) + .collect(); + assert_eq!(cbs, ["CB9", "CB1", "CB4"]); + } + + /// The UMI and gene totals come from the final matrix, not from the read + /// counters, so a cell missing from that map prints zeros rather than + /// inheriting a read count. + #[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()); + let mut umi_gene = BTreeMap::new(); + umi_gene.insert(2u32, (17u32, 5u32)); + let out = st.render(|i| format!("CB{i}"), &umi_gene); + let rows: Vec> = out + .lines() + .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"]); + } + + /// Merging partials keeps first-appearance order, so a threaded run and a + /// serial one write the same file. + #[test] + fn merging_partials_preserves_order_and_sums() { + let mut a = CellReadStats::new(); + a.add_cell(5, &flag_perfect_counted()); + a.add_cell(2, &flag_perfect_counted()); + let mut b = CellReadStats::new(); + b.add_cell(2, &flag_perfect_counted()); + b.add_cell(8, &flag_perfect_counted()); + b.add_no_cb(&CellReadFlag::default()); + a.merge(&b); + + let out = a.render(|i| format!("CB{i}"), &BTreeMap::new()); + let cbs: Vec<&str> = out + .lines() + .skip(2) + .map(|l| l.split('\t').next().unwrap()) + .collect(); + assert_eq!(cbs, ["CB8", "CB2", "CB5"]); + let cb2: Vec<&str> = out + .lines() + .find(|l| l.starts_with("CB2\t")) + .unwrap() + .split('\t') + .collect(); + assert_eq!(cb2[1], "2", "the two partials' reads are summed"); + } +} diff --git a/src/solo/count.rs b/src/solo/count.rs index 9af3edf..cd45309 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -1224,6 +1224,27 @@ pub fn write_gene_matrix( if gzip { " [gzip]" } else { "" }, ); + // `--soloCellReadStats CB`: the per-cell read summary, alongside the + // raw matrix because its UMI and gene columns are the raw totals. + if let Some(cell_stats) = &ctx.cell_read_stats { + let umi_gene: std::collections::BTreeMap = mstats + .cells + .iter() + .map(|c| (c.cb, (c.n_umis as u32, c.n_genes))) + .collect(); + let path = feature_dir.join("CellReads.stats"); + let text = cell_stats.lock().unwrap().render( + |cb| { + ctx.whitelist + .barcode_string(cb) + .unwrap_or_else(|| cb.to_string()) + }, + &umi_gene, + ); + std::fs::write(&path, text).map_err(|e| Error::io(e, &path))?; + log::info!("STARsolo: wrote {}/CellReads.stats", feature.dir_name()); + } + // Filtered (cell-called) matrix per --soloCellFilter. EmptyDrops_CR runs // the Monte-Carlo rescue (needs the per-cell profiles in the body). let called = if params diff --git a/src/solo/mod.rs b/src/solo/mod.rs index 2253dee..e53f4f4 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -9,6 +9,7 @@ //! The barcode read is the SECOND `--readFilesIn` file (STAR convention: //! `--readFilesIn cDNA_read barcode_read`). It is never aligned — only parsed. +pub mod cell_reads; pub mod count; pub mod gene; pub mod smartseq; @@ -564,6 +565,13 @@ pub struct SoloContext { /// `--soloMultiMappers` includes a non-`Unique` method → capture gene- /// ambiguous reads for distribution into `UniqueAndMult-*.mtx`. pub want_multi: bool, + /// `--soloCellReadStats CB`: the per-cell read summary, or `None` when the + /// flag is off. Behind a mutex like the other per-read collections; the + /// lock is only taken when the flag asks for the file. + pub cell_read_stats: Option>, + /// Chromosome indices named by `--genomeChrSetMitochondrial`, for the + /// `mito` column. + pub mito_chr: std::collections::HashSet, } /// Per-region read tallies for the `Summary.csv` mapping funnel (uniquely-mapped @@ -695,6 +703,14 @@ impl SoloContext { velocyto_enabled, velocyto_records: Mutex::new(Vec::new()), want_multi, + cell_read_stats: (params.solo_cell_read_stats == "CB") + .then(|| Mutex::new(crate::solo::cell_reads::CellReadStats::new())), + mito_chr: params + .genome_chr_set_mitochondrial + .iter() + .filter(|n| n.as_str() != "-") + .filter_map(|n| genome.chr_name.iter().position(|c| c == n)) + .collect(), }) } @@ -871,9 +887,67 @@ impl SoloContext { fo }) .collect(); + + self.record_cell_read( + cb_resolved, + &cb_match, + n_loci, + &class, + cdna_transcripts, + &out, + ); out } + /// Fold one read into `CellReads.stats`, when `--soloCellReadStats CB` + /// asked for it. + /// + /// Called at the end of read processing so the counted flags reflect what + /// the read actually produced, rather than what it looked eligible for. + #[allow(clippy::too_many_arguments)] + fn record_cell_read( + &self, + cb_resolved: Option, + cb_match: &CbMatch, + n_loci: usize, + class: &crate::solo::gene::ReadClass, + transcripts: &[Transcript], + out: &SoloReadOutcome, + ) { + let Some(stats) = &self.cell_read_stats else { + return; + }; + let counted_u = out.per_feature.iter().any(|f| f.record.is_some()); + let counted_m = out.per_feature.iter().any(|f| f.multi_gene.is_some()); + let feature_u = counted_u || out.per_feature.iter().any(|f| f.multi.is_some()); + let flag = crate::solo::cell_reads::CellReadFlag { + cb_perfect: matches!(cb_match, CbMatch::Exact(_)), + cb_mm_unique: matches!(cb_match, CbMatch::Corrected(_)), + cb_mm_multiple: matches!(cb_match, CbMatch::Multi(_)), + genome_u: n_loci == 1, + genome_m: n_loci > 1, + feature_u, + feature_m: counted_m, + // The region columns split by strand: STAR reports an antisense + // read under `exonicAS`/`intronicAS`, not under `exonic`/`intronic`. + exonic: !class.antisense && class.region == Some(Region::Exonic), + intronic: !class.antisense && class.region == Some(Region::Intronic), + exonic_as: class.antisense && class.region == Some(Region::Exonic), + intronic_as: class.antisense && class.region == Some(Region::Intronic), + mito: !self.mito_chr.is_empty() + && transcripts + .iter() + .any(|t| self.mito_chr.contains(&t.chr_idx)), + counted_u, + counted_m, + }; + let mut stats = stats.lock().unwrap(); + match cb_resolved { + Some(cb) => stats.add_cell(cb, &flag), + None => stats.add_no_cb(&flag), + } + } + /// Process one 5' paired-end solo read (`--soloBarcodeMate 1`): the barcode is /// from mate 1, and both mates align as a pair. Genes are assigned from the /// union of both mates evaluated against the pair's (mate 1's) transcription From 6ae1de5261b3cc8ea42f17f227f34924e372d726 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 10:24:56 +0200 Subject: [PATCH 2/4] docs: record the CellReads.stats row order in DIVERGENCE.md Section 3.2, in the format CONTRIBUTING.md asks for. --- DIVERGENCE.md | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/DIVERGENCE.md b/DIVERGENCE.md index bd957ed..f1d9f98 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -65,6 +65,20 @@ On the 10k yeast PE benchmark, 4 reads differ in alignment score (AS) because 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. + +**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. + +**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. + +**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. + +**Source.** `src/solo/cell_reads.rs`, locked by `rows_are_emitted_in_reverse_first_appearance_order` and `merging_partials_preserves_order_and_sums`. STAR: `SoloFeature_statsOutput.cpp`. + +--- + ## 4. Implementation divergences (no intended output difference) These differ in *how* a result is produced, not *what* is produced. They are documented so a reviewer chasing a discrepancy knows the mechanism differs by design. From 59fd91f0cba32849eef2462ad929a41660bc4625 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 11:31:17 +0200 Subject: [PATCH 3/4] docs(changelog): keep only this PR's entry CONTRIBUTING.md requires the description to match the code; the entries for the other themes split out of #152 belong to their own PRs. --- CHANGELOG.md | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a17dd8f..5d7fcd2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -97,16 +97,6 @@ Sections commonly used: Features, Bug fixes, Other changes. glibc's malloc and per-thread heaps that return whole segments to the OS when abandoned, so allocator cache size stays bounded. -- **`--soloCellFilter EmptyDrops_CR` now uses CellRanger's actual - statistics.** The ambient profile is smoothed with Simple Good-Turing, - as CellRanger and STAR do, instead of an approximation that reserved - unseen mass from the singleton rate and spread the remainder in - proportion to raw counts. The Monte-Carlo null is drawn with libc++'s - `std::mt19937` and `std::discrete_distribution`, seeded - `19760110 * (isim + 1)` per simulation as STAR seeds it, replacing a - SplitMix64 stream that could not agree with STAR's over an arbitrary - number of draws. Cell calls move as a result. - - **`--soloCellReadStats CB`** writes `Solo.out//CellReads.stats`: one row per cell barcode with fifteen counters describing what happened to its reads — barcode match quality, unique or multi From 8bd9c3d69428fbd1aa6993263274d48c100c8f27 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 11:53:34 +0200 Subject: [PATCH 4/4] refactor(solo): drop CellReadStats::merge, which nothing calls Reads are folded in under a mutex, so there are no per-thread partials to merge; the function was reachable only from its own test. CONTRIBUTING.md rules out shipping a function no production path reaches, and the PR description claimed its test as evidence of thread-safety that the mutex actually provides. --- src/solo/cell_reads.rs | 50 +++--------------------------------------- 1 file changed, 3 insertions(+), 47 deletions(-) diff --git a/src/solo/cell_reads.rs b/src/solo/cell_reads.rs index e17a540..5d140a5 100644 --- a/src/solo/cell_reads.rs +++ b/src/solo/cell_reads.rs @@ -23,6 +23,9 @@ //! **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. use std::collections::{BTreeMap, HashSet}; use std::fmt::Write as _; @@ -124,24 +127,6 @@ impl CellReadStats { Self::fold(&mut self.no_cb, flag); } - /// Merge another accumulator in, preserving first-appearance order: the - /// per-thread partials are merged in read order, so the combined order is - /// the order a single-threaded run would have produced. - pub fn merge(&mut self, other: &Self) { - for &cb in &other.order { - if self.seen.insert(cb) { - self.order.push(cb); - } - let dst = self.cells.entry(cb).or_insert([0; 15]); - for (d, s) in dst.iter_mut().zip(&other.cells[&cb]) { - *d += s; - } - } - for (d, s) in self.no_cb.iter_mut().zip(&other.no_cb) { - *d += s; - } - } - /// Render the file. `umi_gene` gives each cell its final /// `(nUMIunique, nGenesUnique)`; a cell absent from it prints zeros. /// `barcode_of` renders a whitelist index as its barcode string. @@ -300,33 +285,4 @@ mod tests { assert_eq!(rows[1][0], "CB2"); assert_eq!(&rows[1][16..18], ["17", "5"]); } - - /// Merging partials keeps first-appearance order, so a threaded run and a - /// serial one write the same file. - #[test] - fn merging_partials_preserves_order_and_sums() { - let mut a = CellReadStats::new(); - a.add_cell(5, &flag_perfect_counted()); - a.add_cell(2, &flag_perfect_counted()); - let mut b = CellReadStats::new(); - b.add_cell(2, &flag_perfect_counted()); - b.add_cell(8, &flag_perfect_counted()); - b.add_no_cb(&CellReadFlag::default()); - a.merge(&b); - - let out = a.render(|i| format!("CB{i}"), &BTreeMap::new()); - let cbs: Vec<&str> = out - .lines() - .skip(2) - .map(|l| l.split('\t').next().unwrap()) - .collect(); - assert_eq!(cbs, ["CB8", "CB2", "CB5"]); - let cb2: Vec<&str> = out - .lines() - .find(|l| l.starts_with("CB2\t")) - .unwrap() - .split('\t') - .collect(); - assert_eq!(cb2[1], "2", "the two partials' reads are summed"); - } }