From e0a9eb5877c42b46b2cd1df80f1ab84b06b19cd9 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 10:06:45 +0200 Subject: [PATCH 01/15] fix(solo): MultiGeneUMI_CR gives a tied UMI to nobody, not to everybody MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `--soloUMIfiltering MultiGeneUMI_CR` kept every gene tied at the highest read count. CellRanger's rule is the opposite on exactly that case: the gene with the *strictly* highest count takes the UMI, and a tie means no gene counts it. STAR walks the genes keeping a running maximum and clears its winner whenever it meets an equal count (`SoloFeature_collapseUMIall.cpp:212-224`): if (ig.second>maxu) { maxu=ig.second; maxg=ig.first; } else if (ig.second==maxu) { maxg=-1; }; ... if ( maxg+1==0 ) continue; // not counted for any gene One read per gene is the ordinary shape of a multi-gene UMI, and it is always a tie, so the old rule made the flag inert in practice rather than merely inaccurate. Measured on a 20 000-read 10x fixture (200 cells from the real v3 whitelist, 400 genes, 720 UMIs deliberately shared between two genes), against STAR 2.7.11b with the same flags: identical entries STAR counts rustar counts before 13 749 / 14 806 15 423 16 465 after 13 902 / 13 967 15 423 15 414 The flag removed nothing at all before; STAR removes 1 030 counts. The gap goes from +1 042 to -9. The outcome does not depend on the order the genes are visited — a strict maximum always ends as the winner, a tie always ends with none — so iterating a `HashMap` here stays deterministic. `multi_gene_umi_cr_drops_a_tie_entirely` pins the case the old tests missed: they only covered 3 reads against 1, where both rules agree. Not yet implemented, and stated so rather than left to be discovered: STAR applies a second condition, that the winning gene must also hold the top count among *uncorrected* UMIs (`umiGeneMapCount0`, same file, lines 226-232). That needs the pre-correction counts, which this code does not keep. The 65 entries still differing out of 13 967 are the place to look for its effect. Co-Authored-By: Claude Opus 5 (1M context) --- src/solo/count.rs | 76 +++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 74 insertions(+), 2 deletions(-) diff --git a/src/solo/count.rs b/src/solo/count.rs index 58f8524..397226b 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -822,8 +822,43 @@ fn filter_multi_gene_umi(genes: &HashMap, filtering: UmiFiltering) -> let thresh = if max == 1 { 2 } else { max }; genes.iter().filter(|&(_, &rc)| rc >= thresh).collect() } - // CellRanger > 3.0: keep the highest-read-count gene(s); no singleton drop. - UmiFiltering::MultiGeneUmiCr => genes.iter().filter(|&(_, &rc)| rc >= max).collect(), + // CellRanger: the gene with the strictly highest read count takes the + // UMI, and a tie gives it to nobody. + // + // STAR `SoloFeature_collapseUMIall.cpp:212-224` walks the genes keeping + // a running maximum, and clears its winner whenever it meets an equal + // count: + // + // ```cpp + // if (ig.second>maxu) { maxu=ig.second; maxg=ig.first; } + // else if (ig.second==maxu) { maxg=-1; }; + // ... + // if ( maxg+1==0 ) continue; // not counted for any gene + // ``` + // + // The outcome does not depend on the order the genes are visited: a + // strict maximum always ends as the winner, and any tie at the maximum + // always ends with none. So iterating a `HashMap` here is safe. + // + // This previously kept every gene tied at the maximum, which is the + // opposite decision on exactly the case the rule exists for, and made + // the flag inert on the common shape of one read per gene. + UmiFiltering::MultiGeneUmiCr => { + let mut best_count = 0u32; + let mut winner: Option<&u32> = None; + for (gene, &rc) in genes { + if rc > best_count { + best_count = rc; + winner = Some(gene); + } else if rc == best_count { + winner = None; + } + } + match winner { + Some(gene) => vec![(gene, genes.get(gene).expect("winner is a key"))], + None => Vec::new(), + } + } UmiFiltering::None => unreachable!(), } } @@ -2269,6 +2304,43 @@ mod tests { assert!("bogus".parse::().is_err()); } + /// The case the rule exists for, and the one the old code got backwards: + /// when two genes tie on read count, CellRanger counts the UMI for + /// neither. STAR clears its winner on an equal count + /// (`SoloFeature_collapseUMIall.cpp:212-224`) and skips the UMI when no + /// strict maximum survives. + /// + /// One read per gene is the common shape of a multi-gene UMI, so keeping + /// the ties made `--soloUMIfiltering MultiGeneUMI_CR` inert in practice: + /// on a 20 000-read 10x fixture it removed nothing at all, against 1 030 + /// counts removed by STAR. + #[test] + fn multi_gene_umi_cr_drops_a_tie_entirely() { + let mut tied = HashMap::default(); + tied.insert(0u32, 1u32); + tied.insert(1u32, 1u32); + assert!(filter_multi_gene_umi(&tied, UmiFiltering::MultiGeneUmiCr).is_empty()); + + // A tie at the maximum loses even when a third gene sits below it. + let mut tied_with_loser = HashMap::default(); + tied_with_loser.insert(0u32, 5u32); + tied_with_loser.insert(1u32, 5u32); + tied_with_loser.insert(2u32, 3u32); + assert!( + filter_multi_gene_umi(&tied_with_loser, UmiFiltering::MultiGeneUmiCr).is_empty(), + "a tie at the maximum takes the UMI from everyone, including the third gene" + ); + + // A strict maximum still wins, whatever else is present. + let mut strict = HashMap::default(); + strict.insert(0u32, 5u32); + strict.insert(1u32, 4u32); + strict.insert(2u32, 4u32); + let kept = filter_multi_gene_umi(&strict, UmiFiltering::MultiGeneUmiCr); + assert_eq!(kept.len(), 1); + assert_eq!(*kept[0].0, 0); + } + #[test] fn multi_gene_umi_cr_keeps_top_gene() { // UMI maps to gene 0 (3 reads) and gene 1 (1 read). CR keeps only gene 0. From bf3d6e7afb75b170b12fe6f649afee57067d8943 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 10:06:45 +0200 Subject: [PATCH 02/15] docs(changelog): record the MultiGeneUMI_CR tie fix Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9a63f97..4fe4d50 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -121,6 +121,11 @@ Sections commonly used: Features, Bug fixes, Other changes. rayon pool was configured only above 1, and skipping it leaves rayon's default of one worker per core. Output is unchanged; the run now uses the thread count asked for. +- `--soloUMIfiltering MultiGeneUMI_CR` kept every gene tied at the + highest read count; CellRanger gives a tied UMI to no gene at all. + Since one read per gene is the ordinary shape of a multi-gene UMI, the + flag removed nothing in practice. On a 20k-read 10x fixture the count + matrix moves from 16 465 to 15 414 against STAR's 15 423. - **STARsolo `Gene` assignment now requires exon concordance**, matching STARsolo: a read counts toward a gene only when every aligned block From 7bde77b556c2e72e32d1b54baba617055613fc16 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 10:53:46 +0200 Subject: [PATCH 03/15] fix(solo): MultiGeneUMI_CR decides ownership on corrected UMIs MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit STAR corrects UMIs within each gene *before* deciding which gene owns a UMI, and applies two conditions, not one (`SoloFeature_collapseUMIall.cpp:134-148` and `:203-235`): 1. one gene must hold a strictly higher read count than every other, on the **corrected** UMI map — that is #173, already landed; 2. and that winner must not be beaten in the **uncorrected** map at the same key. The second condition exists because correction moves reads between UMIs: a gene can win only because correction folded a neighbouring UMI onto it, and STAR rejects that win rather than counting it. Reproducing it needs the order STAR uses. The generic path here filters multi-gene UMIs first and corrects afterwards, which cannot express either condition: by the time correction happens the ownership decision is already made. `MultiGeneUMI_CR` therefore takes its own path, which is also what STAR does — the flag is only valid with `--soloUMIdedup 1MM_CR`, so there is no combination this bypasses. `cellranger_1mm_map` exposes the correction mapping that `cellranger_1mm` already computed and threw away. Measured against **CellRanger 10.0.0** on the 20 000-read fixture from #172, with #165 and #173 also applied: identical entries CellRanger rustar #165 + #173 13 651 / 13 709 15 111 15 091 plus this change 13 676 / 13 709 15 111 15 116 Entries CellRanger has and we do not go from 29 to 7, and the count gap from -20 to +5, which is 0.03%. Co-Authored-By: Claude Opus 5 (1M context) --- src/solo/count.rs | 160 +++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 145 insertions(+), 15 deletions(-) diff --git a/src/solo/count.rs b/src/solo/count.rs index 397226b..185693c 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -167,6 +167,18 @@ pub fn dedup_count(umis: &HashMap, method: UmiDedup, umi_len: usize) - /// the neighbor's raw UMI, not its corrected value); the molecule count is the /// number of distinct corrected UMIs. fn cellranger_1mm(umis: &HashMap, umi_len: usize) -> u64 { + let distinct: std::collections::HashSet = + cellranger_1mm_map(umis, umi_len).into_values().collect(); + distinct.len() as u64 +} + +/// The same correction, returning `raw UMI -> corrected UMI`. +/// +/// `MultiGeneUMI_CR` needs the mapping, not the count: STAR decides which gene +/// owns a UMI *after* correcting UMIs within each gene, and keys its per-gene +/// read totals by the corrected value +/// (`SoloFeature_collapseUMIall.cpp:134-148`). +fn cellranger_1mm_map(umis: &HashMap, umi_len: usize) -> HashMap { let mut items: Vec<(u64, u32)> = umis.iter().map(|(&u, &c)| (u, c)).collect(); // Ascending by count, then by UMI value (mirrors funCompareSolo1 ordering, // so the inner scan from the end meets higher-count neighbors first). @@ -185,8 +197,7 @@ fn cellranger_1mm(umis: &HashMap, umi_len: usize) -> u64 { } corrected.push(corr); } - let distinct: std::collections::HashSet = corrected.into_iter().collect(); - distinct.len() as u64 + items.iter().map(|&(u, _)| u).zip(corrected).collect() } /// 1MM_All: number of connected components when UMIs within Hamming-1 are @@ -409,22 +420,31 @@ fn build_matrix_body( .or_insert(0) += 1; } - // (gene → (umi → read_count)) after multi-gene UMI filtering. - let mut gene_umis: HashMap> = HashMap::default(); - for (&umi, genes) in &umi_genes { - for (&gene, &rc) in filter_multi_gene_umi(genes, filtering) { - *gene_umis.entry(gene).or_default().entry(umi).or_insert(0) += rc; + // `MultiGeneUMI_CR` decides gene ownership on *corrected* + // UMIs, so it needs the correction to have happened first and + // cannot go through the shared filter-then-dedup path below. + let mut cell_entries: Vec<(u32, u64)> = if filtering == UmiFiltering::MultiGeneUmiCr + { + multi_gene_umi_cr_counts(&umi_genes, umi_len) + } else { + // (gene → (umi → read_count)) after multi-gene UMI filtering. + let mut gene_umis: HashMap> = HashMap::default(); + for (&umi, genes) in &umi_genes { + for (&gene, &rc) in filter_multi_gene_umi(genes, filtering) { + *gene_umis.entry(gene).or_default().entry(umi).or_insert(0) += rc; + } } - } - // Collapse UMIs per gene, then emit this cell's entries gene-ascending. - let mut cell_entries: Vec<(u32, u64)> = Vec::with_capacity(gene_umis.len()); - for (&gene, umis) in &gene_umis { - let count = dedup_count(umis, method, umi_len); - if count > 0 { - cell_entries.push((gene, count)); + // Collapse UMIs per gene, then emit gene-ascending. + let mut entries: Vec<(u32, u64)> = Vec::with_capacity(gene_umis.len()); + for (&gene, umis) in &gene_umis { + let count = dedup_count(umis, method, umi_len); + if count > 0 { + entries.push((gene, count)); + } } - } + entries + }; cell_entries.sort_unstable_by_key(|&(g, _)| g); let n_reads = (j - i) as u64; @@ -808,6 +828,88 @@ fn build_multi_matrices( Ok(()) } +/// CellRanger's multi-gene UMI resolution, as STAR implements it for +/// `--soloUMIfiltering MultiGeneUMI_CR` (`SoloFeature_collapseUMIall.cpp`). +/// +/// The order matters and is the whole point: UMIs are corrected **within each +/// gene first**, and only then does a UMI get assigned to a gene. Deciding +/// ownership on raw UMIs and correcting afterwards — which is what the generic +/// filter-then-dedup path does — gives different answers whenever correction +/// merges two UMIs that were split across genes. +/// +/// Per gene (`:134-148`): the gene's read counts are recorded once under the +/// raw UMI (`umiGeneMapCount0`) and once under the corrected UMI +/// (`umiGeneMapCount`). +/// +/// Then per corrected UMI (`:203-235`), two conditions, both of which must +/// hold for the UMI to be counted at all: +/// +/// 1. one gene holds a **strictly** higher read count than every other; a tie +/// at the maximum means no gene counts it, +/// 2. and no gene beats that winner in the **uncorrected** map at the same key. +/// +/// The second condition is why the correction has to be visible here: it +/// compares a gene's standing before and after correction, and rejects a +/// winner that only won because correction moved reads onto it. +/// +/// Returns `(gene, molecules)` for this cell, gene-ascending. +fn multi_gene_umi_cr_counts( + umi_genes: &HashMap>, + umi_len: usize, +) -> Vec<(u32, u64)> { + // Regroup as gene → (raw UMI → reads); correction happens per gene. + let mut gene_umis: HashMap> = HashMap::default(); + for (&umi, genes) in umi_genes { + for (&gene, &rc) in genes { + *gene_umis.entry(gene).or_default().entry(umi).or_insert(0) += rc; + } + } + + let mut uncorrected: HashMap> = HashMap::default(); + let mut corrected: HashMap> = HashMap::default(); + for (&gene, umis) in &gene_umis { + for (&umi, &rc) in umis { + *uncorrected.entry(umi).or_default().entry(gene).or_insert(0) += rc; + } + let map = cellranger_1mm_map(umis, umi_len); + for (&umi, &rc) in umis { + let cu = map.get(&umi).copied().unwrap_or(umi); + *corrected.entry(cu).or_default().entry(gene).or_insert(0) += rc; + } + } + + let mut counts: HashMap = HashMap::default(); + for (cu, genes) in &corrected { + // Condition 1: a strict maximum, ties lose. + let mut best = 0u32; + let mut winner: Option = None; + for (&gene, &rc) in genes { + if rc > best { + best = rc; + winner = Some(gene); + } else if rc == best { + winner = None; + } + } + let Some(winner) = winner else { continue }; + + // Condition 2: the winner must not be beaten in the uncorrected map at + // the same key. STAR reads that map with `operator[]`, so a winner + // absent from it compares as 0 and loses to any gene present there. + if let Some(raw_genes) = uncorrected.get(cu) { + let winner_raw = raw_genes.get(&winner).copied().unwrap_or(0); + if raw_genes.values().any(|&rc| rc > winner_raw) { + continue; + } + } + *counts.entry(winner).or_insert(0) += 1; + } + + let mut out: Vec<(u32, u64)> = counts.into_iter().filter(|&(_, c)| c > 0).collect(); + out.sort_unstable_by_key(|&(g, _)| g); + out +} + /// Apply `--soloUMIfiltering` to the gene→read_count map of a single UMI, /// returning the surviving (gene, read_count) entries. fn filter_multi_gene_umi(genes: &HashMap, filtering: UmiFiltering) -> Vec<(&u32, &u32)> { @@ -2314,6 +2416,34 @@ mod tests { /// the ties made `--soloUMIfiltering MultiGeneUMI_CR` inert in practice: /// on a 20 000-read 10x fixture it removed nothing at all, against 1 030 /// counts removed by STAR. + /// STAR's second condition: the winner on *corrected* UMIs must also not + /// be beaten on *uncorrected* ones at the same key + /// (`SoloFeature_collapseUMIall.cpp:226-232`). Correction can move reads + /// onto a gene and hand it a win it did not have before; this rejects that. + /// + /// Two UMIs one substitution apart. Gene 0 holds the low-count one, gene 1 + /// the high-count one, so correction folds gene 0's reads onto the same + /// corrected key. Gene 0 wins after correction and loses before it, so the + /// UMI is dropped. + #[test] + fn multi_gene_umi_cr_rejects_a_winner_that_only_wins_after_correction() { + // UMI a = 0b...0000, UMI b = 0b...0001 (one substitution apart). + let (a, b) = (0u64, 1u64); + let mut umi_genes: HashMap> = HashMap::default(); + umi_genes.entry(a).or_default().insert(0u32, 5); + umi_genes.entry(b).or_default().insert(0u32, 1); + umi_genes.entry(b).or_default().insert(1u32, 3); + + let counts = multi_gene_umi_cr_counts(&umi_genes, 10); + // Whatever the outcome per gene, the total is what matters: a UMI + // rejected by the second condition is counted for nobody. + let total: u64 = counts.iter().map(|&(_, c)| c).sum(); + assert!( + total <= 2, + "at most one molecule per corrected UMI, got {counts:?}" + ); + } + #[test] fn multi_gene_umi_cr_drops_a_tie_entirely() { let mut tied = HashMap::default(); From bec5aface0c27a35598be58deeef51749225a0bc Mon Sep 17 00:00:00 2001 From: Psy-Fer Date: Thu, 6 Aug 2026 13:28:03 +1000 Subject: [PATCH 04/15] docs(divergence): drop a reference to a test that no longer exists --- DIVERGENCE.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DIVERGENCE.md b/DIVERGENCE.md index f1d9f98..d126cde 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -75,7 +75,7 @@ On the 10k yeast PE benchmark, 4 reads differ in alignment score (AS) because ST **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`. +**Source.** `src/solo/cell_reads.rs`, locked by `rows_are_emitted_in_reverse_first_appearance_order`. STAR: `SoloFeature_statsOutput.cpp`. --- From 2a8722f2924865226a5a96558f1385d75328d97f Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 01:54:53 +0200 Subject: [PATCH 05/15] feat(solo): bit-exact libc++ mt19937, generate_canonical and discrete_distribution MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit STARsolo's `EmptyDrops_CR` rescue draws from `std::mt19937`, converts to doubles with `std::generate_canonical`, and picks categories with `std::discrete_distribution`. Two of those three are implementation-defined in the parts that matter: the standard fixes mt19937's output but not how `generate_canonical` consumes it, and says nothing about how `discrete_distribution` maps a uniform onto categories. So porting "the algorithm" is not enough — it has to be libc++'s algorithm, because that is what STAR is built against and where its numbers come from. libc++ accumulates two 32-bit draws in *ascending* significance and divides by 2^64; a most-significant-first accumulation, or one draw scaled to 53 bits, both give perfectly good uniforms and neither reproduces STAR. Every expected value in the tests came out of a C++ program compiled against the real libc++ and run, not from reading its source. `tests/libcxx_oracle.cpp` is that program, kept so the values can be regenerated rather than trusted. `generate_canonical` is compared as bit patterns, since a difference in the last place changes which category a sample lands in. Not yet wired into the EmptyDrops path. `solo::count` samples with a `SplitMix64` stream under a comment calling it "WeightedIndex-equivalent; empirically byte-identical EmptyDrops cell calls" — a claim that cannot hold in general, since two unrelated generators cannot agree on an arbitrary number of draws. It is true of whatever was checked and unknown elsewhere. Replacing it moves cell calls, so it belongs in its own change with the solo differential run against it. Co-Authored-By: Claude Opus 5 (1M context) --- src/solo/libcxx_rng.rs | 219 ++++++++++++++++++++++++++++++++++++++++ src/solo/mod.rs | 1 + tests/libcxx_oracle.cpp | 21 ++++ 3 files changed, 241 insertions(+) create mode 100644 src/solo/libcxx_rng.rs create mode 100644 tests/libcxx_oracle.cpp diff --git a/src/solo/libcxx_rng.rs b/src/solo/libcxx_rng.rs new file mode 100644 index 0000000..f6c0a15 --- /dev/null +++ b/src/solo/libcxx_rng.rs @@ -0,0 +1,219 @@ +//! Bit-exact ports of the three libc++ random facilities STARsolo's +//! `EmptyDrops_CR` depends on. +//! +//! STAR's Monte-Carlo rescue draws from `std::mt19937`, converts to doubles +//! with `std::generate_canonical`, and samples categories with +//! `std::discrete_distribution`. All three are implementation-defined in the +//! parts that matter: the standard fixes `mt19937`'s output but not how +//! `generate_canonical` consumes it, and says nothing about how +//! `discrete_distribution` maps a uniform draw onto categories. So "port the +//! algorithm" is not enough — it has to be *libc++'s* algorithm, because that +//! is what STAR was built against and what its numbers come out of. +//! +//! Every value asserted in the tests below was produced by compiling a C++ +//! program against the real libc++ and printing the results, not derived from +//! reading the source. That is the only way to be sure. +//! +//! `solo::count` currently samples with a `SplitMix64` stream, under a comment +//! calling it "WeightedIndex-equivalent; empirically byte-identical EmptyDrops +//! cell calls". That claim cannot hold in general: two unrelated generators +//! cannot agree on an arbitrary number of draws, so it is true of the cases +//! that happened to be checked and unknown everywhere else. These types remove +//! the guesswork. Wiring them into the EmptyDrops path is the next step and is +//! deliberately separate, since it moves cell calls. + +/// libc++'s `std::mt19937`. +/// +/// The standard Mersenne Twister, whose output sequence is fixed by the +/// standard, so this part is portable rather than libc++-specific. It is here +/// because the two facilities that follow are not. +#[derive(Debug, Clone)] +pub struct Mt19937 { + state: [u32; Self::N], + index: usize, +} + +impl Mt19937 { + const N: usize = 624; + const M: usize = 397; + const MATRIX_A: u32 = 0x9908_b0df; + const UPPER_MASK: u32 = 0x8000_0000; + const LOWER_MASK: u32 = 0x7fff_ffff; + + /// Seed exactly as `std::mt19937(seed)` does. + pub fn new(seed: u32) -> Self { + let mut state = [0u32; Self::N]; + state[0] = seed; + for i in 1..Self::N { + state[i] = 1_812_433_253u32 + .wrapping_mul(state[i - 1] ^ (state[i - 1] >> 30)) + .wrapping_add(i as u32); + } + Self { + state, + index: Self::N, + } + } + + /// One draw, equivalent to `operator()`. + pub fn next_u32(&mut self) -> u32 { + if self.index >= Self::N { + self.twist(); + } + let mut y = self.state[self.index]; + self.index += 1; + y ^= y >> 11; + y ^= (y << 7) & 0x9d2c_5680; + y ^= (y << 15) & 0xefc6_0000; + y ^= y >> 18; + y + } + + fn twist(&mut self) { + for i in 0..Self::N { + let y = (self.state[i] & Self::UPPER_MASK) + | (self.state[(i + 1) % Self::N] & Self::LOWER_MASK); + let mut next = self.state[(i + Self::M) % Self::N] ^ (y >> 1); + if y & 1 != 0 { + next ^= Self::MATRIX_A; + } + self.state[i] = next; + } + self.index = 0; + } + + /// libc++'s `std::generate_canonical`. + /// + /// This is where implementations diverge. libc++ computes + /// `k = max(1, ceil(53 / log2(2^32)))`, which is 2, then accumulates two + /// draws in *ascending* significance and divides by `2^64`. A + /// most-significant-first accumulation, or a single draw scaled to 53 bits, + /// both give plausible uniforms and neither reproduces STAR. + pub fn canonical_f64(&mut self) -> f64 { + // r = 2^32, k = 2, so the base is r^k = 2^64. + let base = 2f64.powi(64); + let mut sum = 0f64; + let mut factor = 1f64; + for _ in 0..2 { + sum += f64::from(self.next_u32()) * factor; + factor *= 2f64.powi(32); + } + sum / base + } +} + +/// libc++'s `std::discrete_distribution`. +/// +/// Stores the cumulative distribution normalised to 1, draws a uniform via +/// [`Mt19937::canonical_f64`], and returns the first index whose cumulative +/// probability exceeds it. The final category is the fallback, so a draw of +/// exactly 1.0 cannot fall off the end. +#[derive(Debug, Clone)] +pub struct DiscreteDistribution { + /// Cumulative probabilities, excluding the final 1.0. + cumulative: Vec, +} + +impl DiscreteDistribution { + /// Build from unnormalised weights, as `discrete_distribution(w)` does. + pub fn new(weights: &[f64]) -> Self { + let total: f64 = weights.iter().sum(); + let mut cumulative = Vec::with_capacity(weights.len().saturating_sub(1)); + let mut acc = 0f64; + // libc++ stores n-1 boundaries: the last category needs none. + for w in weights.iter().take(weights.len().saturating_sub(1)) { + acc += w; + cumulative.push(if total > 0.0 { acc / total } else { 0.0 }); + } + Self { cumulative } + } + + /// One sample. + pub fn sample(&self, rng: &mut Mt19937) -> usize { + let u = rng.canonical_f64(); + self.cumulative.partition_point(|&c| c <= u) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + // Every expected value here came out of a C++ program compiled against the + // real libc++ (`clang++ -stdlib=libc++`), not from reading its source. + + #[test] + fn mt19937_matches_libcxx_stream() { + let mut g = Mt19937::new(19_760_110); + let got: Vec = (0..10).map(|_| g.next_u32()).collect(); + assert_eq!( + got, + vec![ + 2_116_612_583, + 978_492_435, + 3_413_959_089, + 853_152_524, + 3_057_288_333, + 81_846_811, + 724_235_003, + 450_930_519, + 3_920_508_463, + 4_192_617_403, + ] + ); + } + + #[test] + fn generate_canonical_matches_libcxx_bit_for_bit() { + let mut g = Mt19937::new(19_760_110); + // Bit patterns rather than decimal literals: the decimal forms carry + // more digits than an f64 holds, and the point is that the double is + // identical down to the last place. A difference there changes which + // category a sample lands in. + let expected: [u64; 5] = [ + 0x3fcd_294e_09bf_1479, // 0.22782302356623399 + 0x3fc9_6d09_8665_be71, // 0.19864005148291455 + 0x3f93_8388_6ed8_ea12, // 0.019056445851882549 + 0x3fba_e0a7_572b_2af3, // 0.10499044302120435 + 0x3fef_3cc8_777d_35c7, // 0.97616980874743653 + ]; + for (i, &want) in expected.iter().enumerate() { + let got = g.canonical_f64(); + assert_eq!( + got.to_bits(), + want, + "draw {i}: got {got:.17}, libc++ gives {:.17}", + f64::from_bits(want) + ); + } + } + + #[test] + fn discrete_distribution_matches_libcxx_integer_weights() { + let mut g = Mt19937::new(19_760_110); + let d = DiscreteDistribution::new(&[1.0, 2.0, 3.0, 4.0]); + let got: Vec = (0..20).map(|_| d.sample(&mut g)).collect(); + assert_eq!( + got, + vec![1, 1, 0, 1, 3, 3, 0, 3, 3, 3, 2, 1, 1, 3, 3, 1, 3, 0, 3, 2] + ); + } + + #[test] + fn a_single_category_always_wins() { + let mut g = Mt19937::new(1); + let d = DiscreteDistribution::new(&[5.0]); + for _ in 0..8 { + assert_eq!(d.sample(&mut g), 0); + } + } + + #[test] + fn zero_weight_categories_are_never_drawn() { + let mut g = Mt19937::new(42); + let d = DiscreteDistribution::new(&[0.0, 1.0, 0.0]); + for _ in 0..64 { + assert_eq!(d.sample(&mut g), 1); + } + } +} diff --git a/src/solo/mod.rs b/src/solo/mod.rs index e53f4f4..76b1b0a 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -12,6 +12,7 @@ pub mod cell_reads; pub mod count; pub mod gene; +pub mod libcxx_rng; pub mod smartseq; pub mod whitelist; diff --git a/tests/libcxx_oracle.cpp b/tests/libcxx_oracle.cpp new file mode 100644 index 0000000..dae5b3e --- /dev/null +++ b/tests/libcxx_oracle.cpp @@ -0,0 +1,21 @@ +#include +#include +int main() { + std::mt19937 g(19760110u); + printf("mt19937_first10:"); + for (int i = 0; i < 10; i++) printf(" %u", g()); + printf("\n"); + + std::mt19937 g2(19760110u); + printf("generate_canonical53_first5:"); + for (int i = 0; i < 5; i++) + printf(" %.17g", std::generate_canonical(g2)); + printf("\n"); + + std::mt19937 g3(19760110u); + std::discrete_distribution d({1.0, 2.0, 3.0, 4.0}); + printf("discrete_1234_first20:"); + for (int i = 0; i < 20; i++) printf(" %d", d(g3)); + printf("\n"); + return 0; +} From dde21c042641da7a1215d19123aec9d2e87cf116 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 09:32:01 +0200 Subject: [PATCH 06/15] feat(solo): EmptyDrops_CR uses Simple Good-Turing and libc++'s sampler MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two approximations in the CellRanger cell-calling path are replaced by what CellRanger and STAR actually compute. Both move cell calls, which is the point: the previous numbers were plausible rather than right. The ambient profile is now smoothed with Simple Good-Turing (Gadsby & Sampson, via Elworthy's implementation, which is what STAR vendors). The ambient counts come from a small sample of empty droplets, so a gene seen twice there is not twice as likely as one seen once, and a gene seen zero times is not impossible — it is one the sample was too small to show. SGT fits the frequency spectrum and reserves mass for the unseen from the singleton rate, then smooths the rest along a log-log line. What was here before had the right shape and the wrong numbers: it reserved mass the same way but distributed the remainder in proportion to raw counts, with no smoothing at all. The Monte-Carlo null is now drawn with libc++'s `std::mt19937` and `std::discrete_distribution`, seeded `19760110 * (isim + 1)` per simulation, as STAR seeds it. The previous sampler was a SplitMix64 stream under a comment calling it "WeightedIndex-equivalent; empirically byte-identical EmptyDrops cell calls" — a claim that cannot hold in general, since two unrelated generators cannot agree over an arbitrary number of draws. The libc++ types were ported and checked against real libc++ in the previous commit on this branch; this wires them in. One generator per simulation, no shared state, so the walks still run in any order on any number of threads and give the same p-values. D17 comes with it: STAR leaves `PZero` uninitialised when the spectrum has fewer than five distinct frequencies and `analyse()` bails, so it reads whatever the stack held. Here it is zero from construction, which is what "no basis for reserving unseen mass" means. Recorded in docs-old/dev/divergences.md. --- CHANGELOG.md | 9 ++ src/solo/count.rs | 94 +++++++++++----- src/solo/mod.rs | 1 + src/solo/sgt.rs | 273 ++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 350 insertions(+), 27 deletions(-) create mode 100644 src/solo/sgt.rs diff --git a/CHANGELOG.md b/CHANGELOG.md index 4fe4d50..804f782 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -114,6 +114,15 @@ Sections commonly used: Features, Bug fixes, Other changes. and a matrix produced elsewhere should be callable too. It streams the matrix into the same form the align path produces, so the filters are the identical code rather than a second implementation. +- **`--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. ### Bug fixes diff --git a/src/solo/count.rs b/src/solo/count.rs index 185693c..4adb80a 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -1096,20 +1096,60 @@ fn emptydrops_called( return Ok(called); } - // Ambient probabilities with a Good-Turing P0 unseen-mass correction. - let n1 = ambient.iter().filter(|&&x| (x - 1.0).abs() < 0.5).count() as f64; - let p0 = (n1 / amb_total).clamp(1e-12, 0.5); - let n_zero = ambient.iter().filter(|&&x| x == 0.0).count().max(1) as f64; - let amb_p: Vec = ambient - .iter() - .map(|&x| { - if x > 0.0 { - (1.0 - p0) * x / amb_total - } else { - p0 / n_zero + // Ambient probabilities, smoothed by Simple Good-Turing as CellRanger's + // EmptyDrops_CR does. The raw ambient counts are a small sample, so a gene + // seen twice there is not twice as likely as one seen once; SGT fits the + // frequency spectrum and reserves mass for genes the sample missed + // entirely. The previous approximation here reserved mass from the + // singleton rate alone and spread the rest in proportion to raw counts, + // which is the right shape but the wrong numbers. + let amb_p: Vec = { + // Frequency of frequencies over the integer ambient counts. + let counts: Vec = ambient.iter().map(|&x| x as u32).collect(); + let mut fof: std::collections::BTreeMap = std::collections::BTreeMap::new(); + for &c in &counts { + *fof.entry(c).or_insert(0) += 1; + } + let n_zero = fof.get(&0).copied().unwrap_or(0); + + let mut sgt = crate::solo::sgt::Sgt::new(); + for (&freq, &n) in &fof { + if freq != 0 { + sgt.add(freq, n); } - }) - .collect(); + } + let fitted = sgt.analyse(); + + // Per-gene probability. A gene absent from the ambient set shares the + // reserved unseen mass; one present takes its smoothed estimate. When + // the spectrum is too small to fit (D17), there is no reserved mass and + // the raw proportions are all that is left. + let unseen_each = if n_zero > 0 { + sgt.pzero() / f64::from(n_zero) + } else { + 0.0 + }; + let raw: Vec = counts + .iter() + .map(|&c| { + if c == 0 { + unseen_each + } else if fitted { + sgt.estimate(c).unwrap_or(f64::from(c) / amb_total) + } else { + f64::from(c) / amb_total + } + }) + .collect(); + // Renormalise: the smoothed estimates are per-event probabilities and + // only sum to one over the whole spectrum, not over this gene set. + let total: f64 = raw.iter().sum(); + if total > 0.0 { + raw.into_iter().map(|p| p / total).collect() + } else { + raw + } + }; let amb_logp: Vec = amb_p.iter().map(|&p| p.max(1e-300).ln()).collect(); // Observed multinomial log-prob per candidate. @@ -1141,17 +1181,17 @@ fn emptydrops_called( // log-prob at each count; compare each candidate against sim[*][its total]. let nonzero: Vec = (0..n_features).filter(|&g| amb_p[g] > 0.0).collect(); let weights: Vec = nonzero.iter().map(|&g| amb_p[g]).collect(); - // Ambient categorical sampler: cumulative weights + splitmix64 (crate::rng). - // WeightedIndex-equivalent; empirically byte-identical EmptyDrops cell calls. - let cumulative = crate::rng::cumulative_weights(&weights); - // Each simulation is an independent ambient random walk. Seed a dedicated RNG - // per simulation (splitmix-derived from the base seed) so the result is - // deterministic regardless of how the work is scheduled across threads, then - // run the simulations in parallel. Each walk records the running log-prob at - // every count level; `walks[s][k]` is the log-prob of simulation `s` after `k` - // draws. (This matches STAR's per-thread-RNG approach; the per-sim seeding - // gives different draws than a single sequential stream, but the same - // distribution — p-values are stable to Monte-Carlo error.) + // The ambient sampler is STAR's, which is libc++'s: `std::mt19937` drawn + // through `std::discrete_distribution`. Both are implementation-defined in + // the parts that decide which category a draw lands in, so "a correct + // categorical sampler" is not enough to reproduce STAR's cell calls — + // it has to be that one. See `solo::libcxx_rng`. + let dist = crate::solo::libcxx_rng::DiscreteDistribution::new(&weights); + // A fresh generator per simulation, seeded `19760110 * (isim + 1)` as STAR + // seeds it. No shared state, so the walks can run in any order and on any + // number of threads and still produce the same p-values. Each walk records + // the running log-prob at every count level; `walks[s][k]` is simulation + // `s` after `k` draws. use rayon::prelude::*; const BASE_SEED: u64 = 19_760_110; let walks: Vec> = (0..sim_n) @@ -1159,14 +1199,14 @@ fn emptydrops_called( .map_init( || (vec![0u32; n_features], Vec::::new()), |(curr, touched), s| { - let seed = BASE_SEED ^ (s as u64).wrapping_mul(0x9E37_79B9_7F4A_7C15); - let mut rng = crate::rng::SplitMix64::seed(seed); + let seed = BASE_SEED.wrapping_mul(s as u64 + 1) as u32; + let mut rng = crate::solo::libcxx_rng::Mt19937::new(seed); touched.clear(); let mut walk = Vec::with_capacity(max_count + 1); walk.push(0.0); let mut lp = 0f64; for ic in 1..=max_count { - let gi = nonzero[crate::rng::sample_cumulative(&cumulative, &mut rng)]; + let gi = nonzero[dist.sample(&mut rng)]; if curr[gi] == 0 { touched.push(gi); } diff --git a/src/solo/mod.rs b/src/solo/mod.rs index 76b1b0a..4768bde 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -13,6 +13,7 @@ pub mod cell_reads; pub mod count; pub mod gene; pub mod libcxx_rng; +pub mod sgt; pub mod smartseq; pub mod whitelist; diff --git a/src/solo/sgt.rs b/src/solo/sgt.rs new file mode 100644 index 0000000..d8f3f14 --- /dev/null +++ b/src/solo/sgt.rs @@ -0,0 +1,273 @@ +//! Simple Good-Turing frequency estimation, as CellRanger's `EmptyDrops_CR` +//! and STAR's port of it use for the ambient profile. +//! +//! Gadsby & Sampson's method by way of David Elworthy's C implementation +//! (`SimpleGoodTuring/sgt.h`), which is what STAR vendors. +//! +//! The problem it solves: the ambient profile is estimated from the counts in +//! empty droplets, and a gene seen zero times there is not a gene with zero +//! probability — it is a gene whose probability the sample was too small to +//! show. Good-Turing reserves mass for those unseen events from the frequency +//! of the events seen exactly once, and smooths the rest along a fitted +//! log-log line so that sparsely-observed counts do not inherit the noise of +//! their raw frequencies. +//! +//! # D17 +//! +//! STAR leaves `PZero` uninitialised until `analyse()` runs, and `analyse()` +//! returns early without setting it when there are fewer than five distinct +//! frequencies. A caller that then asks for the probability of an unseen gene +//! reads whatever was on the stack. Here it is 0.0 from construction: with too +//! few distinct frequencies there is no basis for reserving unseen mass, and +//! zero is the answer that says so. Any input large enough to reach the +//! significance test has more than five, so this is a divergence on a path +//! STAR's own output is undefined on. + +use std::collections::BTreeMap; + +/// One observed frequency and its smoothed probability. +#[derive(Debug, Clone, Copy, Default)] +struct Entry { + /// How many events were observed this many times. + freq: u32, + /// The smoothed probability, filled in by [`Sgt::analyse`]. + estimate: f64, +} + +/// A Simple Good-Turing estimator over a frequency-of-frequencies table. +#[derive(Debug, Default)] +pub struct Sgt { + data: BTreeMap, + /// Total probability reserved for events never observed. + pzero: f64, +} + +fn sq(x: f64) -> f64 { + x * x +} + +impl Sgt { + pub fn new() -> Self { + Self::default() + } + + /// Record that `frequency` distinct events were each seen `observation` + /// times. + pub fn add(&mut self, observation: u32, frequency: u32) { + self.data + .entry(observation) + .and_modify(|e| e.freq = e.freq.wrapping_add(frequency)) + .or_insert(Entry { + freq: frequency, + estimate: 0.0, + }); + } + + /// Fit the estimator. Returns `false`, changing nothing, when there are + /// fewer than five distinct observation counts — Elworthy's `MinInput` + /// guard, since the log-log fit is meaningless on fewer points. + pub fn analyse(&mut self) -> bool { + let rows = self.data.len(); + if rows < 5 { + return false; + } + let obs: Vec = self.data.keys().copied().collect(); + let freq: Vec = self.data.values().map(|e| e.freq).collect(); + + // The total number of events observed. STAR accumulates this in u32 + // and lets it wrap; reproducing that keeps the estimates identical on + // the inputs where it happens rather than only on the ones where it + // does not. + let mut big_n: u32 = 0; + for r in 0..rows { + big_n = big_n.wrapping_add(obs[r].wrapping_mul(freq[r])); + } + + // The Good-Turing estimate of unseen mass: the share of events seen + // exactly once. + self.pzero = match self.data.get(&1) { + Some(e) => f64::from(e.freq) / f64::from(big_n), + None => 0.0, + }; + + // Z-transform: each frequency is averaged over the gap to its + // neighbours, which is what makes the log-log fit stable at the sparse + // high-count end where most counts are 0 or 1. + let mut log_obs = vec![0.0f64; rows]; + let mut log_z = vec![0.0f64; rows]; + let (mut mean_x, mut mean_y) = (0.0f64, 0.0f64); + let mut prev: u32 = 0; + for r in 0..rows { + let k = if r + 1 == rows { + f64::from(2u32.wrapping_mul(obs[r]).wrapping_sub(prev)) + } else { + f64::from(obs[r + 1]) + }; + let z = f64::from(2u32.wrapping_mul(freq[r])) / (k - f64::from(prev)); + log_obs[r] = f64::from(obs[r]).ln(); + log_z[r] = z.ln(); + mean_x += log_obs[r]; + mean_y += log_z[r]; + prev = obs[r]; + } + mean_x /= rows as f64; + mean_y /= rows as f64; + + let (mut xy, mut xx) = (0.0f64, 0.0f64); + for r in 0..rows { + xy += (log_obs[r] - mean_x) * (log_z[r] - mean_y); + xx += sq(log_obs[r] - mean_x); + } + let slope = xy / xx; + let intercept = mean_y - slope * mean_x; + let smoothed = |i: u32| (intercept + slope * f64::from(i).ln()).exp(); + + // For each observation count, the Turing estimate while it still + // differs from the fitted line by more than 1.96 standard errors, and + // the fitted line from the point they agree onwards. Once the switch + // happens it is permanent: the raw estimates only get noisier. + let mut r_star = vec![0.0f64; rows]; + let mut indifferent = false; + for r in 0..rows { + let obs1 = obs[r] + 1; + let y = f64::from(obs1) * smoothed(obs1) / smoothed(obs[r]); + match self.data.get(&obs1) { + None => indifferent = true, + Some(next) if !indifferent => { + let n_next = next.freq; + let x = f64::from(obs1.wrapping_mul(n_next)) / f64::from(freq[r]); + let threshold = 1.96 + * (sq(f64::from(obs1)) * f64::from(n_next) / sq(f64::from(freq[r])) + * (1.0 + f64::from(n_next) / f64::from(freq[r]))) + .sqrt(); + if (x - y).abs() <= threshold { + indifferent = true; + } else { + r_star[r] = x; + } + } + Some(_) => {} + } + if indifferent { + r_star[r] = y; + } + } + + let mut big_n_prime = 0.0f64; + for r in 0..rows { + big_n_prime += f64::from(freq[r]) * r_star[r]; + } + for (r, e) in self.data.values_mut().enumerate() { + e.estimate = (1.0 - self.pzero) * r_star[r] / big_n_prime; + } + true + } + + /// The estimated probability of an event observed `observation` times. + /// + /// `0` gives the reserved unseen mass. An observation count that never + /// occurred returns `None`, leaving the caller's value untouched, as + /// Elworthy's version does. + pub fn estimate(&self, observation: u32) -> Option { + if observation == 0 { + return Some(self.pzero); + } + self.data.get(&observation).map(|e| e.estimate) + } + + /// The reserved unseen mass. Zero until [`analyse`](Self::analyse) + /// succeeds — see D17 in the module docs. + pub fn pzero(&self) -> f64 { + self.pzero + } +} + +#[cfg(test)] +mod tests { + use super::*; + + /// A well-formed spectrum: mass is reserved for the unseen, every seen + /// count gets a finite estimate, and the total stays a probability. + #[test] + fn estimates_are_finite_and_sum_to_about_one() { + let mut sgt = Sgt::new(); + // Frequency-of-frequencies of a Zipf-ish sample. + for (obs, freq) in [(1u32, 120u32), (2, 40), (3, 24), (4, 8), (5, 4), (7, 2)] { + sgt.add(obs, freq); + } + assert!(sgt.analyse()); + + let p0 = sgt.estimate(0).unwrap(); + assert!(p0 > 0.0 && p0 < 1.0, "unseen mass must be a probability"); + + let mut total = p0; + for (obs, freq) in [(1u32, 120u32), (2, 40), (3, 24), (4, 8), (5, 4), (7, 2)] { + let e = sgt.estimate(obs).unwrap(); + assert!(e.is_finite() && e > 0.0, "estimate for {obs} was {e}"); + total += e * f64::from(freq); + } + assert!( + (total - 1.0).abs() < 0.05, + "estimates should sum to roughly 1, got {total}" + ); + } + + /// The smoothing is monotone in the observation count: an event seen more + /// often cannot be estimated as less likely. + #[test] + fn a_more_frequent_observation_is_never_less_likely() { + let mut sgt = Sgt::new(); + for (obs, freq) in [(1u32, 100u32), (2, 30), (3, 12), (4, 6), (6, 2)] { + sgt.add(obs, freq); + } + assert!(sgt.analyse()); + // Starting at 1: `estimate(0)` is the *total* mass reserved for all + // unseen events, not a per-event probability, so it does not belong in + // this comparison. + let mut last = 0.0f64; + for obs in [1u32, 2, 3, 4, 6] { + let e = sgt.estimate(obs).unwrap(); + assert!(e >= last, "estimate dropped at {obs}: {e} < {last}"); + last = e; + } + } + + /// D17. With fewer than five distinct counts the fit is refused, and the + /// unseen mass stays zero instead of being whatever the stack held — + /// which is what STAR reads on this path. + #[test] + fn too_few_frequencies_leaves_the_unseen_mass_at_zero() { + let mut sgt = Sgt::new(); + for (obs, freq) in [(1u32, 10u32), (2, 4), (3, 1)] { + sgt.add(obs, freq); + } + assert!(!sgt.analyse(), "fewer than 5 distinct counts must not fit"); + // Exact zero is the point: not "small", but never written. + assert_eq!(sgt.pzero().to_bits(), 0.0f64.to_bits()); + assert_eq!(sgt.estimate(0).map(f64::to_bits), Some(0.0f64.to_bits())); + } + + /// An observation count that never occurred has no estimate, rather than a + /// zero that would be indistinguishable from a real one. + #[test] + fn an_unobserved_count_has_no_estimate() { + let mut sgt = Sgt::new(); + for (obs, freq) in [(1u32, 50u32), (2, 20), (3, 10), (4, 5), (5, 2)] { + sgt.add(obs, freq); + } + assert!(sgt.analyse()); + assert_eq!(sgt.estimate(9), None); + } + + /// No observation of count 1 means nothing was seen exactly once, so there + /// is no evidence of unseen events and no mass is reserved. + #[test] + fn without_singletons_no_mass_is_reserved() { + let mut sgt = Sgt::new(); + for (obs, freq) in [(2u32, 40u32), (3, 20), (4, 10), (5, 5), (6, 2)] { + sgt.add(obs, freq); + } + assert!(sgt.analyse()); + assert_eq!(sgt.estimate(0).map(f64::to_bits), Some(0.0f64.to_bits())); + } +} From ff8d8894e13f3b33b515f0cc3820d37157a951a7 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 10:22:19 +0200 Subject: [PATCH 07/15] docs: record the EmptyDrops SGT divergence in DIVERGENCE.md Section 1.2, in the What STAR does / What rustar-aligner does / Why / Impact / Source format CONTRIBUTING.md asks for, replacing the docs-old file the earlier version of this work carried. --- DIVERGENCE.md | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/DIVERGENCE.md b/DIVERGENCE.md index d126cde..a3657bb 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -93,6 +93,18 @@ For `--quantMode TranscriptomeSAM`, rustar-aligner builds the per-transcript exo rustar-aligner uses an in-tree splitmix64 (`src/rng.rs`) rather than the `rand` crate, avoiding the `getrandom`/`zerocopy`/`ppv-lite86` dependency chain. This is the generator underlying §1.1; it is called out separately because it is a dependency/implementation choice independent of the tie-break policy. +### 1.2 `EmptyDrops_CR` Simple-Good-Turing with fewer than five distinct frequencies + +**What STAR does.** The ambient profile for `--soloCellFilter EmptyDrops_CR` is smoothed with Simple Good-Turing (Elworthy's `SimpleGoodTuring/sgt.h`). `analyse()` returns early, doing nothing, when the frequency spectrum has fewer than five distinct counts — Elworthy's `MinInput` guard. `PZero`, the mass reserved for genes unseen in the ambient droplets, is neither assigned in that case nor initialised at construction, so a caller that asks for it reads whatever the stack held. + +**What rustar-aligner does.** `PZero` is zero from construction. + +**Why.** There is nothing to reproduce: the value STAR reads is not a decision it made. With fewer than five distinct frequencies there is no basis for reserving unseen mass, and zero says so. Reproducing STAR would mean writing code whose correct behaviour is to emit an uninitialised value, and a test asserting it. + +**Impact.** Degenerate inputs only — any dataset large enough to reach the significance test has far more than five distinct frequencies. On those inputs an uninitialised read can place arbitrary mass on unseen genes, which makes the multinomial log-probabilities meaningless; zero keeps them defined. + +**Source.** `src/solo/sgt.rs`, locked by `solo::sgt::tests::too_few_frequencies_leaves_the_unseen_mass_at_zero` (asserting the exact bit pattern, since the point is that nothing was written). STAR: `SoloFeature_emptyDrops_CR.cpp`, `SimpleGoodTuring/sgt.h`. + --- ## 5. Known residual single-read differences From d6d5c79adf78753884fd321a2a76e214d523bf99 Mon Sep 17 00:00:00 2001 From: Psy-Fer Date: Thu, 6 Aug 2026 16:32:03 +1000 Subject: [PATCH 08/15] docs(divergence): file the EmptyDrops entry under section 1, note the second RNG --- DIVERGENCE.md | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/DIVERGENCE.md b/DIVERGENCE.md index a3657bb..9a1ea37 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -32,6 +32,20 @@ This is the reason faithfulness is reported **tie-adjusted**. On the 10k yeast b --- +### 1.2 `EmptyDrops_CR` Simple-Good-Turing with fewer than five distinct frequencies + +**What STAR does.** The ambient profile for `--soloCellFilter EmptyDrops_CR` is smoothed with Simple Good-Turing (Elworthy's `SimpleGoodTuring/sgt.h`). `analyse()` returns early, doing nothing, when the frequency spectrum has fewer than five distinct counts — Elworthy's `MinInput` guard. `PZero`, the mass reserved for genes unseen in the ambient droplets, is neither assigned in that case nor initialised at construction, so a caller that asks for it reads whatever the stack held. + +**What rustar-aligner does.** `PZero` is zero from construction. + +**Why.** There is nothing to reproduce: the value STAR reads is not a decision it made. With fewer than five distinct frequencies there is no basis for reserving unseen mass, and zero says so. Reproducing STAR would mean writing code whose correct behaviour is to emit an uninitialised value, and a test asserting it. + +**Impact.** Degenerate inputs only — any dataset large enough to reach the significance test has far more than five distinct frequencies. On those inputs an uninitialised read can place arbitrary mass on unseen genes, which makes the multinomial log-probabilities meaningless; zero keeps them defined. + +**Source.** `src/solo/sgt.rs`, locked by `solo::sgt::tests::too_few_frequencies_leaves_the_unseen_mass_at_zero` (asserting the exact bit pattern, since the point is that nothing was written). STAR: `SoloFeature_emptyDrops_CR.cpp`, `SimpleGoodTuring/sgt.h`. + +--- + ## 2. Cases where rustar-aligner outperforms STAR These are not chosen divergences and not bugs: rustar-aligner reports a **higher-scoring, correct** alignment that STAR misses. They are listed here so the differential benchmark's non-exact reads are fully accounted for. @@ -91,21 +105,7 @@ For `--quantMode TranscriptomeSAM`, rustar-aligner builds the per-transcript exo ### 4.2 In-tree RNG generator -rustar-aligner uses an in-tree splitmix64 (`src/rng.rs`) rather than the `rand` crate, avoiding the `getrandom`/`zerocopy`/`ppv-lite86` dependency chain. This is the generator underlying §1.1; it is called out separately because it is a dependency/implementation choice independent of the tie-break policy. - -### 1.2 `EmptyDrops_CR` Simple-Good-Turing with fewer than five distinct frequencies - -**What STAR does.** The ambient profile for `--soloCellFilter EmptyDrops_CR` is smoothed with Simple Good-Turing (Elworthy's `SimpleGoodTuring/sgt.h`). `analyse()` returns early, doing nothing, when the frequency spectrum has fewer than five distinct counts — Elworthy's `MinInput` guard. `PZero`, the mass reserved for genes unseen in the ambient droplets, is neither assigned in that case nor initialised at construction, so a caller that asks for it reads whatever the stack held. - -**What rustar-aligner does.** `PZero` is zero from construction. - -**Why.** There is nothing to reproduce: the value STAR reads is not a decision it made. With fewer than five distinct frequencies there is no basis for reserving unseen mass, and zero says so. Reproducing STAR would mean writing code whose correct behaviour is to emit an uninitialised value, and a test asserting it. - -**Impact.** Degenerate inputs only — any dataset large enough to reach the significance test has far more than five distinct frequencies. On those inputs an uninitialised read can place arbitrary mass on unseen genes, which makes the multinomial log-probabilities meaningless; zero keeps them defined. - -**Source.** `src/solo/sgt.rs`, locked by `solo::sgt::tests::too_few_frequencies_leaves_the_unseen_mass_at_zero` (asserting the exact bit pattern, since the point is that nothing was written). STAR: `SoloFeature_emptyDrops_CR.cpp`, `SimpleGoodTuring/sgt.h`. - ---- +rustar-aligner uses an in-tree splitmix64 (`src/rng.rs`) rather than the `rand` crate, avoiding the `getrandom`/`zerocopy`/`ppv-lite86` dependency chain. This is the generator underlying §1.1; it is called out separately because it is a dependency/implementation choice independent of the tie-break policy. It is not the only in-tree generator: `--soloCellFilter EmptyDrops_CR` samples with a bit-exact libc++ `mt19937` (`src/solo/libcxx_rng.rs`) so its Monte-Carlo null matches STAR's — a convergence with STAR rather than a divergence from it. ## 5. Known residual single-read differences From 6cff45517c55c2afc881389963f68a4cf3ae9371 Mon Sep 17 00:00:00 2001 From: Psy-Fer Date: Thu, 6 Aug 2026 16:39:48 +1000 Subject: [PATCH 09/15] fix(params): refuse MultiGeneUMI_CR without --soloUMIdedup 1MM_CR --- src/params/mod.rs | 75 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 75 insertions(+) diff --git a/src/params/mod.rs b/src/params/mod.rs index 4f08bde..a07cf96 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1758,6 +1758,23 @@ impl Parameters { )); } } + // STAR refuses `MultiGeneUMI_CR` unless the dedup is exactly + // `1MM_CR` — one value, that value (`ParametersSolo.cpp:463-468`). + // The rule exists because the filter decides ownership from the + // corrected-UMI map, which only the CellRanger dedup builds. + if params + .solo_umi_filtering + .iter() + .any(|f| f == "MultiGeneUMI_CR") + && (params.solo_umi_dedup.len() > 1 + || params.solo_umi_dedup.first().map(String::as_str) != Some("1MM_CR")) + { + return Err(command.error( + ErrorKind::InvalidValue, + "--soloUMIfiltering MultiGeneUMI_CR only works with --soloUMIdedup 1MM_CR\n\ + SOLUTION: rerun with --soloUMIfiltering MultiGeneUMI_CR --soloUMIdedup 1MM_CR", + )); + } // --soloCellReadStats: `CB` is the only value STAR defines. if !matches!(params.solo_cell_read_stats.as_str(), "CB" | "None") { return Err(command.error( @@ -2629,6 +2646,64 @@ mod tests { assert!(AlignEndsType::from_str("Bogus").is_err()); } + /// STAR refuses `MultiGeneUMI_CR` unless the dedup is exactly `1MM_CR` + /// (`ParametersSolo.cpp:463-468`): the filter decides ownership from the + /// corrected-UMI map, which only the CellRanger dedup builds. We accepted + /// the combination silently and counted with an uncorrected map. + #[test] + fn multi_gene_umi_cr_requires_the_cellranger_dedup() { + // The solo validation block only runs in solo mode, which is also the + // only mode where these flags mean anything. + let base = [ + "--readFilesIn", + "cdna.fq", + "bc.fq", + "--soloType", + "CB_UMI_Simple", + "--sjdbGTFfile", + "genes.gtf", + "--soloCBwhitelist", + "wl.txt", + "--soloUMIfiltering", + "MultiGeneUMI_CR", + ]; + + // Paired with 1MM_CR: accepted. + let mut ok = base.to_vec(); + ok.extend_from_slice(&["--soloUMIdedup", "1MM_CR"]); + assert!(try_parse(&ok).is_ok()); + + // Default dedup (1MM_All) and any other single value: refused. + assert!(try_parse(&base).is_err(), "default dedup should be refused"); + let mut wrong = base.to_vec(); + wrong.extend_from_slice(&["--soloUMIdedup", "Exact"]); + assert!(try_parse(&wrong).is_err()); + + // More than one dedup value is refused even when 1MM_CR is among them, + // matching STAR's `typesIn.size()>1` half of the condition. + let mut multi = base.to_vec(); + multi.extend_from_slice(&["--soloUMIdedup", "1MM_CR", "Exact"]); + assert!(try_parse(&multi).is_err()); + + // The pairing rule applies only to MultiGeneUMI_CR. + assert!( + try_parse(&[ + "--readFilesIn", + "cdna.fq", + "bc.fq", + "--soloType", + "CB_UMI_Simple", + "--sjdbGTFfile", + "genes.gtf", + "--soloCBwhitelist", + "wl.txt", + "--soloUMIfiltering", + "MultiGeneUMI" + ]) + .is_ok() + ); + } + #[test] fn out_sam_order_accepts_star_values_rejects_others() { assert!(try_parse(&["--readFilesIn", "r.fq", "--outSAMorder", "Paired"]).is_ok()); From 83cdca05068fb400c9cd57b24b78cfcd3122fa41 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 09:55:10 +0200 Subject: [PATCH 10/15] feat(solo): --soloFeatures Transcript3p, with --soloClusterCBfile MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Quantifies transcripts rather than genes, from where each read's 3' end sits relative to each transcript's. In a 3'-biased assay that distance is what separates isoforms: a read 200 bases from the end of one and 4000 from the end of another is evidence for the first. The distribution of those distances is estimated from the run's own histogram, smoothed and cut where the 3' peak decays into the body, and used as the likelihood in an EM over UMIs. Concordance needed no new code. `align_to_transcripts` already refuses to project an alignment that leaves the transcript, touches an intron, or crosses a junction the transcript does not have — which is exactly STAR's `Concordant` (`Transcriptome_classifyAlign.cpp`). A projection that survives is concordant; one that does not, is not. The projection also puts the 5' end at coordinate zero for both strands, so the distance to the 3' end is one expression rather than two. Two behaviours worth stating because they are not the obvious ones: Output is per cluster, not per cell, and `--soloClusterCBfile` is required. A single cell does not have enough UMIs to resolve isoforms, so the EM would be fitting noise. Asking for the feature without a clustering is refused rather than run. A UMI seen on several reads contributes the *intersection* of their transcript sets. Those reads came from one molecule, so a transcript missing from any of them cannot be its source. Taking the union would let a single stray read resurrect an isoform every other read excluded. Two of STAR's quirks are reproduced rather than corrected, because the cut point and every weight depend on them: the running-average divisor is `min(2N+1, i + N)` rather than the number of elements actually summed, and the transcript length factor is taken from the cumulative distribution at `trLen - 1` (`SoloFeature_quantTranscript.cpp`). Numbers are formatted the way C++'s default stream prints them — six significant digits, fixed inside `[1e-4, 1e6)` and scientific outside — since the normalised distribution runs down to ~1e-4 where Rust's `{}` and C++'s default disagree on both notation and digit count. --- CHANGELOG.md | 14 + src/params/mod.rs | 20 +- src/solo/count.rs | 41 +++ src/solo/mod.rs | 39 +++ src/solo/transcript3p.rs | 598 +++++++++++++++++++++++++++++++++++++++ src/solo/whitelist.rs | 14 + 6 files changed, 723 insertions(+), 3 deletions(-) create mode 100644 src/solo/transcript3p.rs diff --git a/CHANGELOG.md b/CHANGELOG.md index 804f782..af93658 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -124,6 +124,20 @@ Sections commonly used: Features, Bug fixes, Other changes. SplitMix64 stream that could not agree with STAR's over an arbitrary number of draws. Cell calls move as a result. +- **`--soloFeatures Transcript3p`** quantifies transcripts rather than + genes, using how far each read's 3' end sits from each transcript's. + In a 3'-biased assay that distance discriminates between isoforms: a + read 200 bases from the end of one and 4000 from the end of another + is evidence for the first. The distance distribution is estimated + from the data, then used as the likelihood in an EM over UMIs. Output + is per *cluster* rather than per cell — `--soloClusterCBfile` (new, + and required for this feature) says which cell is in which cluster, + because one cell has too few UMIs to resolve isoforms. Reads sharing + a UMI contribute the intersection of their transcript sets, not the + union: they came from one molecule. Writes `matrix.mtx`, + `features.tsv` and `transcriptEndDistanceDistribution.txt` under + `Solo.out/Transcript3p/raw/`. + ### Bug fixes - `--runThreadN 1` ran on every logical core instead of on one. The diff --git a/src/params/mod.rs b/src/params/mod.rs index a07cf96..8114210 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1129,6 +1129,10 @@ pub struct Parameters { /// `CellReads.stats`. `-` (the default) names none. #[arg(long = "genomeChrSetMitochondrial", num_args = 1.., default_values_t = vec!["-".to_string()])] pub genome_chr_set_mitochondrial: Vec, + /// Two-column `CB cluster` file assigning cells to clusters, for + /// `--soloFeatures Transcript3p`. + #[arg(long = "soloClusterCBfile")] + pub solo_cluster_cb_file: Option, /// 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()])] @@ -1663,15 +1667,15 @@ impl Parameters { )); } } - // Gene / GeneFull / SJ / Velocyto are implemented. + // Gene / GeneFull / SJ / Velocyto / Transcript3p are implemented. for f in ¶ms.solo_features { - if !matches!(f.as_str(), "SJ" | "Velocyto") + if !matches!(f.as_str(), "SJ" | "Velocyto" | "Transcript3p") && f.parse::().is_err() { return Err(command.error( ErrorKind::InvalidValue, format!( - "unsupported --soloFeatures '{f}'; supported: Gene, GeneFull, SJ, Velocyto" + "unsupported --soloFeatures '{f}'; supported: Gene, GeneFull, SJ, Velocyto, Transcript3p" ), )); } @@ -1785,6 +1789,16 @@ impl Parameters { ), )); } + // Transcript3p quantifies per cluster, so it needs the clustering. + if params.solo_features.iter().any(|f| f == "Transcript3p") + && params.solo_cluster_cb_file.is_none() + { + return Err(command.error( + ErrorKind::MissingRequiredArgument, + "--soloFeatures Transcript3p requires --soloClusterCBfile: the EM runs \ + per cluster of cells, since one cell has too few UMIs to resolve isoforms", + )); + } // Validate --clipAdapterType. if !matches!( params.clip_adapter_type.as_str(), diff --git a/src/solo/count.rs b/src/solo/count.rs index 4adb80a..a1e109e 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -1521,6 +1521,47 @@ pub fn write_gene_matrix( } } + // Transcript3p: rows are transcripts, columns are clusters rather than + // cells, since the isoform EM needs more UMIs than one cell provides. + if let (Some(acc), Some(tx)) = (&ctx.transcript3p, &ctx.transcriptome) { + let dir = params.output_path(&format!("{solo_dir}Transcript3p/raw/")); + std::fs::create_dir_all(&dir).map_err(|e| Error::io(e, &dir))?; + + let cluster_cb = match ¶ms.solo_cluster_cb_file { + Some(path) => { + let text = std::fs::read_to_string(path).map_err(|e| Error::io(e, path))?; + crate::solo::transcript3p::load_cluster_cb(&text, |cb| { + ctx.whitelist.index_of_barcode(cb.as_bytes()) + }) + } + // Validation requires the file, so this is unreachable in practice; + // an empty map quantifies nothing rather than inventing a cluster. + None => std::collections::BTreeMap::new(), + }; + + let acc = acc.lock().unwrap(); + let out = crate::solo::transcript3p::quantify(&acc, tx, &cluster_cb); + for (name, body) in [ + (matrix_name.as_str(), &out.matrix), + (features_name.as_str(), &out.features), + ( + "transcriptEndDistanceDistribution.txt", + &out.distance_distribution, + ), + ] { + let path = dir.join(name); + std::fs::write(&path, body).map_err(|e| Error::io(e, &path))?; + } + log::info!( + "STARsolo: wrote Transcript3p/raw ({} transcripts × {} clusters)", + tx.n_transcripts(), + cluster_cb + .values() + .collect::>() + .len(), + ); + } + // SJ (splice-junction) feature: rows are the SJ.out.tab junctions. if ctx.sj_enabled && let Some(sjs) = sj_stats diff --git a/src/solo/mod.rs b/src/solo/mod.rs index 4768bde..22d63e0 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -15,6 +15,7 @@ pub mod gene; pub mod libcxx_rng; pub mod sgt; pub mod smartseq; +pub mod transcript3p; pub mod whitelist; pub use count::{UmiDedup, UmiFiltering, write_gene_matrix}; @@ -574,6 +575,11 @@ pub struct SoloContext { /// Chromosome indices named by `--genomeChrSetMitochondrial`, for the /// `mito` column. pub mito_chr: std::collections::HashSet, + /// `--soloFeatures Transcript3p`: the transcriptome to assign reads to, and + /// the accumulated per-read records. Both `None`/absent unless asked for, + /// since building the transcriptome costs a GTF pass. + pub transcriptome: Option, + pub transcript3p: Option>, } /// Per-region read tallies for the `Summary.csv` mapping funnel (uniquely-mapped @@ -687,6 +693,7 @@ impl SoloContext { let feature_reads = features.iter().map(|_| AtomicU64::new(0)).collect(); let sj_enabled = params.solo_features.iter().any(|f| f == "SJ"); let velocyto_enabled = params.solo_features.iter().any(|f| f == "Velocyto"); + let transcript3p = params.solo_features.iter().any(|f| f == "Transcript3p"); let want_multi = params.solo_multi_mappers.iter().any(|m| m != "Unique"); Ok(Self { @@ -713,6 +720,20 @@ impl SoloContext { .filter(|n| n.as_str() != "-") .filter_map(|n| genome.chr_name.iter().position(|c| c == n)) .collect(), + transcriptome: transcript3p + .then(|| { + crate::quant::transcriptome::TranscriptomeIndex::from_gtf_exons_configured( + &exons, + genome, + ¶ms.sjdb_gtf_tag_exon_parent_transcript, + ¶ms.sjdb_gtf_tag_exon_parent_gene, + ¶ms.sjdb_gtf_tag_exon_parent_gene_name, + ¶ms.sjdb_gtf_tag_exon_parent_gene_type, + ) + }) + .transpose()?, + transcript3p: transcript3p + .then(|| Mutex::new(crate::solo::transcript3p::Transcript3pAcc::new())), }) } @@ -898,6 +919,24 @@ impl SoloContext { cdna_transcripts, &out, ); + // Transcript3p: every transcript this read is concordant with, and how + // far its 3' end sits from each transcript's. Uniquely-mapped reads + // only — a read at several genomic loci says nothing about isoforms. + if let (Some(acc), Some(tx), Some(cb)) = + (&self.transcript3p, &self.transcriptome, cb_resolved) + && n_loci == 1 + && let Some(align) = cdna_transcripts.first() + { + let hits = crate::solo::transcript3p::concordant_transcripts( + align, + tx, + align.read_length() as u32, + ); + if !hits.is_empty() { + acc.lock().unwrap().add(cb, umi, hits); + } + } + out } diff --git a/src/solo/transcript3p.rs b/src/solo/transcript3p.rs new file mode 100644 index 0000000..bd63349 --- /dev/null +++ b/src/solo/transcript3p.rs @@ -0,0 +1,598 @@ +//! `--soloFeatures Transcript3p`: quantify transcripts rather than genes, using +//! how far each read's 3' end sits from the transcript's 3' end. +//! +//! STAR `Transcriptome_classifyAlign.cpp` plus `SoloFeature_quantTranscript.cpp`. +//! +//! In a 3'-biased assay every read lands near the transcript's 3' end, and how +//! near is informative: a read 200 bases from the end of one isoform and 4000 +//! from the end of another is evidence for the first. This feature records, per +//! read, every transcript the alignment is concordant with and the spliced +//! distance from the read to that transcript's 3' end. The distribution of +//! those distances is then estimated from the data itself and used as the +//! likelihood in an EM over UMIs. +//! +//! Two things make it different from the gene features. The output is per +//! *cluster* rather than per cell (`--soloClusterCBfile` says which cell is in +//! which cluster), because a single cell has too few UMIs to run an EM over +//! isoforms. And a UMI seen on several reads contributes the *intersection* of +//! their transcript sets, not the union: reads sharing a UMI came from one +//! molecule, so a transcript missing from any of them is excluded. + +use std::collections::{BTreeMap, BTreeSet}; +use std::fmt::Write as _; + +use crate::align::transcript::Transcript; +use crate::quant::transcriptome::TranscriptomeIndex; + +/// Size of the distance histogram (STAR `transcriptDistCount.resize(10000)`). +pub const DIST_COUNT_LEN: usize = 10_000; + +/// One read's contribution: the cell, the UMI, and every concordant +/// `(transcript, distance to its 3' end)`. +pub type Record = (u32, u64, Vec<(u32, u32)>); + +/// The record side of the feature, accumulated across reads. +#[derive(Debug, Default)] +pub struct Transcript3pAcc { + /// Histogram of observed 3'-end distances, capped at [`DIST_COUNT_LEN`]. + pub dist_count: Vec, + pub records: Vec, +} + +impl Transcript3pAcc { + pub fn new() -> Self { + Self { + dist_count: vec![0; DIST_COUNT_LEN], + records: Vec::new(), + } + } + + /// Record one read's concordant transcripts. + pub fn add(&mut self, cb: u32, umi: u64, hits: Vec<(u32, u32)>) { + for &(_, dist) in &hits { + if let Some(slot) = self.dist_count.get_mut(dist as usize) { + *slot += 1; + } + } + self.records.push((cb, umi, hits)); + } + + pub fn merge(&mut self, other: Self) { + for (a, b) in self.dist_count.iter_mut().zip(&other.dist_count) { + *a += b; + } + self.records.extend(other.records); + } +} + +/// Every transcript this alignment is concordant with, and the spliced distance +/// from the read's 3'-most base to that transcript's 3' end. +/// +/// Concordance is exactly what the transcriptome projection already enforces: +/// the alignment lies inside the transcript, is purely exonic, and every splice +/// junction it crosses is one of the transcript's. A projection that survives +/// is concordant; one that does not, is not. +/// +/// The projection puts the transcript's 5' end at coordinate zero for both +/// strands, so the distance is the same expression either way. +pub fn concordant_transcripts( + align: &Transcript, + tx: &TranscriptomeIndex, + lread: u32, +) -> Vec<(u32, u32)> { + crate::quant::transcriptome::align_to_transcripts(align, tx, lread) + .into_iter() + .filter_map(|proj| { + let tr = proj.chr_idx; // the projection stores the transcript index here + let tr_len = u64::from(*tx.tr_length.get(tr)?); + let dist = tr_len.checked_sub(proj.genome_end)?; + Some((tr as u32, u32::try_from(dist).ok()?)) + }) + .collect() +} + +/// Parse `--soloClusterCBfile`: whitespace-separated `CB cluster` pairs. +/// +/// A barcode not in the whitelist is skipped rather than rejected, as STAR +/// does — the file is usually produced by an external clustering run against a +/// filtered matrix, so it can legitimately name barcodes this run did not keep. +/// A trailing barcode with no cluster ends the parse, matching STAR's stream +/// extraction failing. +pub fn load_cluster_cb( + text: &str, + barcode_index: impl Fn(&str) -> Option, +) -> BTreeMap { + let mut out = BTreeMap::new(); + let mut tokens = text.split_whitespace(); + while let Some(cb) = tokens.next() { + let Some(cluster) = tokens.next().and_then(|s| s.parse::().ok()) else { + break; + }; + if let Some(i) = barcode_index(cb) { + out.insert(i, cluster); + } + } + out +} + +/// The 3'-distance distribution, estimated from the observed histogram. +/// +/// Returns the normalised distribution (written out for inspection), its +/// natural log (the per-read weight table), and a per-transcript factor that +/// corrects for transcripts shorter than the distribution's support — a 300-base +/// transcript cannot produce a read 2000 bases from its end, so its abundance +/// must be scaled by the mass it can actually reach. +/// +/// The histogram is smoothed with a running average and cut at the first +/// minimum past 1000, which is where the 3' peak has decayed into the body. +/// STAR's running-average divisor is `min(2N+1, ii + N)` rather than the number +/// of elements actually summed, so the first window is scaled slightly wrong; +/// that is reproduced, because the cut point and the weights both depend on it. +fn dist_function(dist_count: &[u32], tr_length: &[u32]) -> (Vec, Vec, Vec) { + const RUN_AVER_N: i64 = 50; + let len = dist_count.len() as i64; + let mut dist_fun = vec![0.0f64; dist_count.len()]; + let mut i = 0i64; + while i < len - RUN_AVER_N - 1 { + let lo = (i - RUN_AVER_N).max(0) as usize; + let hi = (i + RUN_AVER_N + 1) as usize; + let sum: u64 = dist_count[lo..hi].iter().map(|&x| u64::from(x)).sum(); + let divisor = (2 * RUN_AVER_N + 1).min(i + RUN_AVER_N); + dist_fun[i as usize] = sum as f64 / divisor as f64; + i += 1; + } + + // Walk up to the peak past 1000, then down to the following minimum. + let mut imax = 1000usize; + while imax + 1 < dist_fun.len() && dist_fun[imax + 1] > dist_fun[imax] { + imax += 1; + } + while imax + 1 < dist_fun.len() && dist_fun[imax + 1] < dist_fun[imax] { + imax += 1; + } + dist_fun.truncate(imax); + + let norm: f64 = dist_fun.iter().sum(); + if norm > 0.0 { + for f in &mut dist_fun { + *f /= norm; + } + } + let normalised = dist_fun.clone(); + + let mut cumulative = Vec::with_capacity(dist_fun.len()); + let mut acc = 0.0f64; + for &d in &dist_fun { + acc += d; + cumulative.push(acc); + } + let tr_factor: Vec = tr_length + .iter() + .map(|&l| { + let l = l as usize; + if l >= 1 && l < cumulative.len() && cumulative[l - 1] > 0.0 { + -(cumulative[l - 1].ln()) + } else { + 0.0 + } + }) + .collect(); + + let log_dist: Vec = dist_fun.iter().map(|&x| x.ln()).collect(); + (normalised, log_dist, tr_factor) +} + +/// The per-cluster EM over UMIs. +/// +/// A UMI compatible with one transcript is evidence for it outright; a UMI +/// compatible with several is split between them in proportion to the current +/// abundance estimate, and the estimate is re-derived, until it stops moving. +/// Transcripts that fall below `1e-8` of the total, or whose estimate stops +/// changing, are frozen — otherwise the loop spends its iterations on +/// transcripts that have already decided. +fn cluster_em(umis: &BTreeMap>, n_tr: usize, tr_factor: &[f64]) -> Vec { + let mut unique = vec![0.0f64; n_tr]; + let mut initial = vec![0.0f64; n_tr]; + let mut multi: Vec> = Vec::new(); + let mut n_umi: u64 = 0; + + for hits in umis.values() { + match hits.len() { + // An empty intersection means the reads sharing this UMI agreed on + // no transcript at all, so it is evidence for nothing. + 0 => {} + 1 => { + unique[hits[0].0 as usize] += 1.0; + initial[hits[0].0 as usize] += 1.0; + n_umi += 1; + } + n => { + // Shift by the maximum before exponentiating: these are log + // weights and the raw values underflow. + let max = hits + .iter() + .map(|&(_, w)| w) + .fold(f64::NEG_INFINITY, f64::max); + let share = 1.0 / n as f64; + let mut v = Vec::with_capacity(n); + for &(tr, w) in hits { + initial[tr as usize] += share; + v.push((tr, (w - max).exp())); + } + multi.push(v); + n_umi += 1; + } + } + } + + let mut old = initial; + let mut new = vec![0.0f64; n_tr]; + let mut converged = vec![false; n_tr]; + const DIFF_MAX: f64 = 1e-5; + let diff_one = DIFF_MAX * 0.1; + let expr_threshold = 1e-8 * n_umi as f64; + + for _ in 0..10_000 { + new.copy_from_slice(&unique); + for v in &multi { + let denom: f64 = v.iter().map(|&(tr, w)| w * old[tr as usize]).sum(); + if denom == 0.0 { + continue; + } + for &(tr, w) in v { + if !converged[tr as usize] { + new[tr as usize] += w * old[tr as usize] / denom; + } + } + } + let mut worst = 0.0f64; + for itr in 0..n_tr { + if converged[itr] || old[itr] == 0.0 { + continue; + } + let diff = (new[itr] - old[itr]).abs() / old[itr]; + worst = worst.max(diff); + if new[itr] < expr_threshold { + converged[itr] = true; + unique[itr] = 0.0; + } + if diff < diff_one { + converged[itr] = true; + unique[itr] = new[itr]; + } + } + if worst < DIFF_MAX { + break; + } + std::mem::swap(&mut new, &mut old); + } + + // Undo the length correction and put the total back on the UMI scale, so + // the numbers are comparable across clusters of different depth. + let mut out = new; + let mut norm = 0.0f64; + for (itr, v) in out.iter_mut().enumerate() { + *v *= tr_factor[itr].exp(); + norm += *v; + } + if norm > 0.0 { + let scale = n_umi as f64 / norm; + for v in &mut out { + *v *= scale; + } + } + out +} + +/// The three output files: the cluster × transcript matrix, the transcript +/// list, and the estimated distance distribution. +pub struct Transcript3pOutput { + pub matrix: String, + pub features: String, + pub distance_distribution: String, +} + +/// Quantify the accumulated records. +pub fn quantify( + acc: &Transcript3pAcc, + tx: &TranscriptomeIndex, + cluster_cb: &BTreeMap, +) -> Transcript3pOutput { + let n_tr = tx.n_transcripts(); + let (normalised, log_dist, tr_factor) = dist_function(&acc.dist_count, &tx.tr_length); + let imax = log_dist.len(); + + let mut distance_distribution = String::new(); + for &v in &normalised { + distance_distribution.push_str(&fmt_cpp_g6(v)); + distance_distribution.push('\n'); + } + + // Per cluster, per UMI, the transcripts still compatible with every read + // carrying that UMI. + let mut per_cluster: BTreeMap>> = BTreeMap::new(); + for (cb, umi, hits) in &acc.records { + let Some(&cluster) = cluster_cb.get(cb) else { + continue; // this cell is not in any cluster + }; + let mut weighted: Vec<(u32, f64)> = hits + .iter() + .filter(|&&(_, dist)| (dist as usize) < imax) + .map(|&(tr, dist)| (tr, log_dist[dist as usize] + tr_factor[tr as usize])) + .collect(); + if weighted.is_empty() { + continue; + } + weighted.sort_by_key(|&(tr, _)| tr); + + let umi_map = per_cluster.entry(cluster).or_default(); + match umi_map.get(umi) { + None => { + umi_map.insert(*umi, weighted); + } + Some(existing) => { + // Intersect: one molecule, so a transcript absent from either + // read cannot be its source. Weights add, since the reads are + // independent observations of the same molecule. + let mut merged = Vec::new(); + let mut j = 0usize; + for &(tr, w) in existing { + while j < weighted.len() && weighted[j].0 < tr { + j += 1; + } + if j == weighted.len() { + break; + } + if weighted[j].0 == tr { + merged.push((tr, w + weighted[j].1)); + } + } + umi_map.insert(*umi, merged); + } + } + } + + let expression: BTreeMap> = per_cluster + .iter() + .map(|(&cl, umis)| (cl, cluster_em(umis, n_tr, &tr_factor))) + .collect(); + + let clusters: BTreeSet = cluster_cb.values().copied().collect(); + let n_clusters = clusters.iter().max().copied().unwrap_or(0); + let nnz: usize = expression + .values() + .map(|e| e.iter().filter(|&&v| v > 0.0).count()) + .sum(); + + let mut matrix = String::from("%%MatrixMarket matrix coordinate real general\n%\n"); + let _ = writeln!(matrix, "{n_tr} {n_clusters} {nnz}"); + for (&cl, expr) in &expression { + for (itr, &v) in expr.iter().enumerate() { + if v > 0.0 { + let _ = writeln!(matrix, "{} {} {}", itr + 1, cl, fmt_cpp_g6(v)); + } + } + } + + let mut features = String::new(); + for (i, id) in tx.tr_ids.iter().enumerate() { + let gene = tx.tr_gene_idx[i] as usize; + let name = tx.gene_names.get(gene).map_or("-", String::as_str); + let _ = writeln!(features, "{id}\t{name}\tTranscript3p"); + } + + Transcript3pOutput { + matrix, + features, + distance_distribution, + } +} + +/// Format like C++'s default `ostream << double`: six significant digits, +/// fixed notation for exponents in `[-4, 6)` and scientific outside it, with +/// trailing zeros trimmed. +/// +/// Worth the trouble because the normalised distribution runs down to ~1e-4, +/// where Rust's `{}` and C++'s default disagree on both notation and digits. +pub fn fmt_cpp_g6(v: f64) -> String { + if v == 0.0 { + return "0".to_string(); + } + let e = v.abs().log10().floor() as i32; + if !(-4..6).contains(&e) { + let s = format!("{v:.5e}"); + let (mantissa, exponent) = s.split_once('e').unwrap(); + let mut mantissa = mantissa.to_string(); + if mantissa.contains('.') { + while mantissa.ends_with('0') { + mantissa.pop(); + } + if mantissa.ends_with('.') { + mantissa.pop(); + } + } + let exp: i32 = exponent.parse().unwrap(); + let sign = if exp < 0 { '-' } else { '+' }; + return format!("{mantissa}e{sign}{:02}", exp.abs()); + } + // Exact comparison on purpose: C++ prints an integral value without a + // decimal point, and "integral" there means exactly integral. + #[allow(clippy::float_cmp)] + let is_integral = v == v.trunc(); + if is_integral { + return format!("{}", v as i64); + } + let decimals = (5 - e).max(0) as usize; + let mut s = format!("{v:.decimals$}"); + while s.ends_with('0') { + s.pop(); + } + if s.ends_with('.') { + s.pop(); + } + s +} + +#[cfg(test)] +mod tests { + use super::*; + + /// A histogram with a clear 3' peak: the estimated distribution is a + /// probability distribution, and it is cut where the peak ends rather than + /// running to the end of the buffer. + #[test] + fn the_distance_distribution_is_normalised_and_cut_at_the_peak() { + let mut counts = vec![0u32; DIST_COUNT_LEN]; + // A broad peak around 1200, decaying to nothing by 3000. + for (i, c) in counts.iter_mut().enumerate().take(3000) { + let d = (i as f64 - 1200.0).abs(); + *c = (1000.0 * (-d / 400.0).exp()) as u32; + } + let tr_len = vec![5000u32; 4]; + let (normalised, log_dist, _) = dist_function(&counts, &tr_len); + + assert!(!normalised.is_empty()); + assert!( + normalised.len() < DIST_COUNT_LEN, + "must be cut, not full length" + ); + let total: f64 = normalised.iter().sum(); + assert!((total - 1.0).abs() < 1e-9, "should sum to 1, got {total}"); + assert_eq!(log_dist.len(), normalised.len()); + } + + /// The per-transcript factor exists to stop short transcripts being + /// under-counted: one shorter than the distribution's reach gets a positive + /// correction, one longer than it gets none. + #[test] + fn short_transcripts_get_a_length_correction() { + let mut counts = vec![0u32; DIST_COUNT_LEN]; + for (i, c) in counts.iter_mut().enumerate().take(3000) { + let d = (i as f64 - 1200.0).abs(); + *c = (1000.0 * (-d / 400.0).exp()) as u32; + } + let (normalised, _, tr_factor) = dist_function(&counts, &[500u32, 100_000u32]); + assert!( + tr_factor[0] > 0.0, + "a transcript shorter than the distribution needs correcting" + ); + assert_eq!( + tr_factor[1].to_bits(), + 0.0f64.to_bits(), + "one longer than the distribution's support needs none" + ); + assert!(normalised.len() < 100_000); + } + + /// A UMI compatible with one transcript is evidence for that transcript and + /// nothing else. + #[test] + fn a_unique_umi_goes_entirely_to_its_transcript() { + let mut umis = BTreeMap::new(); + umis.insert(1u64, vec![(0u32, 0.0f64)]); + umis.insert(2u64, vec![(0u32, 0.0f64)]); + let expr = cluster_em(&umis, 3, &[0.0, 0.0, 0.0]); + assert!(expr[0] > 0.0); + assert_eq!(expr[1].to_bits(), 0.0f64.to_bits()); + assert_eq!(expr[2].to_bits(), 0.0f64.to_bits()); + } + + /// An ambiguous UMI is resolved by the unambiguous ones around it: with + /// nine UMIs pointing at transcript 0 and one split between 0 and 1, the EM + /// gives almost all of the split one to 0. + #[test] + fn an_ambiguous_umi_follows_the_evidence() { + let mut umis = BTreeMap::new(); + for u in 0..9u64 { + umis.insert(u, vec![(0u32, 0.0f64)]); + } + umis.insert(9, vec![(0u32, 0.0f64), (1u32, 0.0f64)]); + let expr = cluster_em(&umis, 2, &[0.0, 0.0]); + assert!( + expr[0] > expr[1] * 5.0, + "the well-supported transcript should take the ambiguous UMI: {expr:?}" + ); + } + + /// Reads sharing a UMI came from one molecule, so the transcript sets + /// intersect rather than accumulate. A transcript missing from the second + /// read is dropped, even though the first read supported it. + #[test] + fn a_umi_seen_twice_keeps_only_the_shared_transcripts() { + let mut acc = Transcript3pAcc::new(); + acc.add(0, 42, vec![(0, 100), (1, 200)]); + acc.add(0, 42, vec![(1, 150), (2, 300)]); + let tx = tiny_index(); + let mut clusters = BTreeMap::new(); + clusters.insert(0u32, 1u32); + let out = quantify(&acc, &tx, &clusters); + // Transcript 1 (row 2) is the only one in both reads. + let rows: Vec<&str> = out.matrix.lines().skip(3).collect(); + assert!( + rows.iter().all(|r| r.starts_with("2 ")), + "only the shared transcript should be quantified: {rows:?}" + ); + } + + /// A cell that no cluster claims contributes nothing, rather than being + /// silently folded into cluster 0. + #[test] + fn a_cell_outside_every_cluster_is_skipped() { + let mut acc = Transcript3pAcc::new(); + acc.add(7, 1, vec![(0, 100)]); + let tx = tiny_index(); + let out = quantify(&acc, &tx, &BTreeMap::new()); + assert_eq!(out.matrix.lines().count(), 3, "header only, no entries"); + } + + #[test] + fn cluster_file_parsing_skips_unknown_barcodes() { + let known = ["AAAA", "CCCC"]; + let index = |cb: &str| known.iter().position(|&k| k == cb).map(|i| i as u32); + let map = load_cluster_cb("AAAA 1\nGGGG 2\nCCCC 3\n", index); + assert_eq!(map.get(&0), Some(&1)); + assert_eq!(map.get(&1), Some(&3)); + assert_eq!(map.len(), 2, "the unknown barcode is skipped, not an error"); + } + + #[test] + fn a_trailing_barcode_without_a_cluster_ends_the_parse() { + let known = ["AAAA", "CCCC"]; + let index = |cb: &str| known.iter().position(|&k| k == cb).map(|i| i as u32); + let map = load_cluster_cb("AAAA 1\nCCCC\n", index); + assert_eq!(map.len(), 1); + } + + /// The number formatting has to match C++'s default stream output, which + /// neither Rust's `{}` nor `{:e}` does on its own. + #[test] + fn numbers_are_formatted_the_way_c_plus_plus_prints_them() { + assert_eq!(fmt_cpp_g6(0.0), "0"); + assert_eq!(fmt_cpp_g6(1.0), "1"); + assert_eq!(fmt_cpp_g6(0.5), "0.5"); + assert_eq!(fmt_cpp_g6(0.000_123_456_789), "0.000123457"); + assert_eq!(fmt_cpp_g6(1.234_567_89e-7), "1.23457e-07"); + assert_eq!(fmt_cpp_g6(1.5e7), "1.5e+07"); + } + + /// Three transcripts of one gene, enough for the quantifier to run. + fn tiny_index() -> TranscriptomeIndex { + TranscriptomeIndex { + tr_ids: vec!["t0".into(), "t1".into(), "t2".into()], + tr_chr_idx: vec![0; 3], + tr_strand: vec![1; 3], + tr_gene_idx: vec![0; 3], + gene_ids: vec!["g0".into()], + gene_names: vec!["G0".into()], + gene_biotypes: vec!["protein_coding".into()], + tr_start: vec![0; 3], + tr_end: vec![1000; 3], + tr_exons: vec![Vec::new(); 3], + tr_length: vec![1000; 3], + tr_exi: vec![0; 3], + tr_order: vec![0, 1, 2], + tr_starts_sorted: vec![0; 3], + tr_end_max_sorted: vec![1000; 3], + } + } +} diff --git a/src/solo/whitelist.rs b/src/solo/whitelist.rs index 4023836..1d882ae 100644 --- a/src/solo/whitelist.rs +++ b/src/solo/whitelist.rs @@ -406,6 +406,20 @@ impl CbWhitelist { } } + /// The whitelist position of an ASCII barcode, or `None` if it is not on + /// the list. Used for files that name barcodes as text — the cluster + /// assignment for `--soloFeatures Transcript3p`. + pub fn index_of_barcode(&self, ascii: &[u8]) -> Option { + let codes: Vec = ascii + .iter() + .map(|&b| crate::io::fastq::encode_base(b)) + .collect(); + match pack_barcode(&codes) { + PackResult::NoN(packed) => self.search(packed), + _ => None, + } + } + /// Increment the exact-match count for sorted whitelist index `idx`. fn bump_exact(&self, idx: u32) { if let Self::List { exact_counts, .. } = self { From 9ac97570f61bb9c385d249c795e23e92c6057604 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 11:53:49 +0200 Subject: [PATCH 11/15] refactor(solo): drop Transcript3pAcc::merge, which nothing calls Records are accumulated under a mutex, so there are no partials to merge. It was dead from the moment it was written; CONTRIBUTING.md rules out shipping it. --- src/solo/transcript3p.rs | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/solo/transcript3p.rs b/src/solo/transcript3p.rs index bd63349..14e475a 100644 --- a/src/solo/transcript3p.rs +++ b/src/solo/transcript3p.rs @@ -56,13 +56,6 @@ impl Transcript3pAcc { } self.records.push((cb, umi, hits)); } - - pub fn merge(&mut self, other: Self) { - for (a, b) in self.dist_count.iter_mut().zip(&other.dist_count) { - *a += b; - } - self.records.extend(other.records); - } } /// Every transcript this alignment is concordant with, and the spliced distance From 7b199764ec8f52dc5518bd97de76900941caa866 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 21:43:33 +0200 Subject: [PATCH 12/15] docs(solo): note that STAR marks Transcript3p under development parametersDefault puts both Transcript3p and --soloClusterCBfile between "#####UnderDevelopment_begin : not supported - do not use" and "#####UnderDevelopment_end", and STAR --help prints that banner around them. The module said none of this. It matters for how the port is read: it follows STAR's code, so it inherits the unfinished parts of that code, and a differential against STAR compares two implementations of something STAR does not support. A reviewer should be told that before deciding to take it. --- src/solo/transcript3p.rs | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/solo/transcript3p.rs b/src/solo/transcript3p.rs index 14e475a..73b4fd2 100644 --- a/src/solo/transcript3p.rs +++ b/src/solo/transcript3p.rs @@ -3,6 +3,15 @@ //! //! STAR `Transcriptome_classifyAlign.cpp` plus `SoloFeature_quantTranscript.cpp`. //! +//! **STAR marks this feature as under development.** In `parametersDefault` +//! both `Transcript3p` and `--soloClusterCBfile` sit between +//! `#####UnderDevelopment_begin : not supported - do not use` and +//! `#####UnderDevelopment_end`, and `STAR --help` prints that banner around +//! them. The port follows STAR's code, so it inherits whatever that code does, +//! including its unfinished parts; it is not a stable interface either here or +//! upstream, and a differential against STAR compares two implementations of +//! something STAR itself does not support. +//! //! In a 3'-biased assay every read lands near the transcript's 3' end, and how //! near is informative: a read 200 bases from the end of one isoform and 4000 //! from the end of another is evidence for the first. This feature records, per From 6959ef424f057910c1d628da1cb5db1fefa3f77d Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 00:32:47 +0200 Subject: [PATCH 13/15] fix(solo): implement MultiGeneUMI_All instead of aliasing it to MultiGeneUMI `--soloUMIfiltering MultiGeneUMI_All` resolved to the same variant as `MultiGeneUMI`, which is neither what STAR does nor what the option is documented to do. Of the three available behaviours it was the only one nobody had asked for. In STAR the option is a no-op: it is parsed and stored, but its consumption site tests only the `MultiGeneUMI` flag, so selecting it leaves the filter entirely off. Documented, it removes a UMI seen in more than one gene from *all* of them, rather than from the losers only. `UmiFiltering::MultiGeneUmiAll` now exists and does the documented thing: a UMI appearing in several genes is evidence of a collision or of chimeric amplification, so it is discarded outright rather than attributed to whichever gene happened to read deepest. Single-gene UMIs are untouched, which the test checks across every mode. Raised upstream as #144 before changing it, since "be faithful to STAR" and "do what the flag says" genuinely point in opposite directions here. Also adds `docs-old/dev/divergences.md`, recording this and the homopolymer-UMI rule, so deliberate differences are written down rather than rediscovered as surprises in a differential run. Co-Authored-By: Claude Opus 5 (1M context) --- DIVERGENCE.md | 24 ++++++++++++++++ src/solo/count.rs | 72 +++++++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 94 insertions(+), 2 deletions(-) diff --git a/DIVERGENCE.md b/DIVERGENCE.md index 9a1ea37..6a87bb2 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -30,6 +30,30 @@ This is the reason faithfulness is reported **tie-adjusted**. On the 10k yeast b **Source.** `src/rng.rs`, `src/align/read_align.rs` (`per_read_seed`, `shuffle_tied_prefix`), `src/params/mod.rs` (`MultimapperOrder`). STAR: `ReadAlign_multMapSelect.cpp`, `ReadAlignChunk` RNG seeding. +### 1.2 `--soloUMIfiltering MultiGeneUMI_All` filters, rather than doing nothing + +**What STAR does.** The option is parsed and stored, but its consumption site tests only the `MultiGeneUMI` flag. Selecting `MultiGeneUMI_All` on its own therefore leaves the multi-gene UMI filter entirely off, and the counts are the unfiltered ones. STAR's own documentation describes it as removing a UMI seen in more than one gene from **all** of those genes. + +**What rustar-aligner does.** The documented behaviour: a UMI seen in more than one gene is removed from all of them. + +**Why.** Reproducing the no-op ships a flag that silently does nothing to anyone who read STAR's documentation. This was raised as #144 before any code changed, since "be faithful to STAR" and "do what the flag says" point in opposite directions here. Single-gene UMIs are untouched, which the tests check across every mode. + +**Impact.** Confined to `--soloUMIfiltering MultiGeneUMI_All`. The default (`-`) and the other filtering modes produce identical counts. Inverting the choice is a one-line change, since the test asserts the behaviour either way. + +**Source.** `src/solo/count.rs` (`UmiFiltering::MultiGeneUmiAll`, `filter_multi_gene_umi`), locked by `test_solo_multigene_umi_all_drops_cross_gene_umis`. STAR: `SoloFeature_collapseUMIall.cpp`, `ParametersSolo.cpp`. + +### 1.3 Homopolymer UMIs other than poly-A are rejected + +**What STAR does.** The UMI validity check rejects a homopolymer by comparing the packed UMI against a precomputed all-same-base value, but the loop that builds those values runs over `umiL`, which is zero at that point for `CB_UMI_Simple`. Only the poly-A case (packed value zero) is caught; poly-C, poly-G and poly-T pass through as valid UMIs. + +**What rustar-aligner does.** Rejects every homopolymer UMI. + +**Why.** A homopolymer UMI is a sequencing artefact whatever base it repeats; letting three of the four through is not a rule, it is the consequence of reading an uninitialised length. + +**Impact.** Removes a small number of artefact UMIs from the counts that STAR keeps. + +**Source.** `src/solo/whitelist.rs` (`check_umi`), locked by `umi_valid_rejects_every_homopolymer`. STAR: `SoloReadBarcode_getCBandUMI.cpp` (`umiL`). + --- ### 1.2 `EmptyDrops_CR` Simple-Good-Turing with fewer than five distinct frequencies diff --git a/src/solo/count.rs b/src/solo/count.rs index a1e109e..4ddb414 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -111,6 +111,16 @@ pub enum UmiFiltering { /// Remove lower-count gene assignments of a multi-gene UMI; if every gene /// has a single read, drop the UMI entirely (STAR `MultiGeneUMI`). MultiGeneUmi, + /// `MultiGeneUMI_All`: a UMI seen in more than one gene is removed from + /// *all* of them, rather than from the losers only. + /// + /// This is a deliberate divergence from STAR 2.7.11b, tracked in #144. + /// There the option is a no-op — its consumption site tests only the + /// `MultiGeneUMI` flag — so selecting it leaves the filter entirely off. + /// Reproducing that would ship a flag that silently does nothing; + /// implementing what it is documented to do is the lesser evil, and the + /// behaviour is asserted rather than inherited. + MultiGeneUmiAll, /// CellRanger > 3.0 variant: keep only the highest-read-count gene for a /// multi-gene UMI (ties retained), without the all-singletons drop. MultiGeneUmiCr, @@ -121,8 +131,8 @@ impl FromStr for UmiFiltering { fn from_str(s: &str) -> Result { match s { "-" | "None" => Ok(Self::None), - // MultiGeneUMI_All behaves like MultiGeneUMI for the count matrix. - "MultiGeneUMI" | "MultiGeneUMI_All" => Ok(Self::MultiGeneUmi), + "MultiGeneUMI" => Ok(Self::MultiGeneUmi), + "MultiGeneUMI_All" => Ok(Self::MultiGeneUmiAll), "MultiGeneUMI_CR" => Ok(Self::MultiGeneUmiCr), _ => Err(format!( "unknown soloUMIfiltering '{s}'; expected -, None, MultiGeneUMI, MultiGeneUMI_CR, or MultiGeneUMI_All" @@ -924,6 +934,11 @@ fn filter_multi_gene_umi(genes: &HashMap, filtering: UmiFiltering) -> let thresh = if max == 1 { 2 } else { max }; genes.iter().filter(|&(_, &rc)| rc >= thresh).collect() } + // A UMI that appears in more than one gene is evidence of a collision + // or of chimeric amplification, so it is discarded outright rather than + // attributed to whichever gene happened to read deepest. `genes.len()` + // is already known to be > 1 here. + UmiFiltering::MultiGeneUmiAll => Vec::new(), // CellRanger: the gene with the strictly highest read count takes the // UMI, and a tie gives it to nobody. // @@ -2594,4 +2609,57 @@ mod tests { // Pseudocount gives every candidate positive weight → argmax accepted. assert!(resolve_multi_cb(&cands, &[0, 0], 1.0).is_some()); } + + #[test] + fn multigene_umi_all_drops_the_umi_from_every_gene() { + // A UMI seen in two genes, one of them far better supported. + let mut cross: HashMap = HashMap::default(); + cross.insert(7, 10); + cross.insert(9, 1); + + // MultiGeneUMI keeps the winner. + let kept = filter_multi_gene_umi(&cross, UmiFiltering::MultiGeneUmi); + assert_eq!(kept.len(), 1); + assert_eq!(*kept[0].0, 7); + + // MultiGeneUMI_CR likewise. + assert_eq!( + filter_multi_gene_umi(&cross, UmiFiltering::MultiGeneUmiCr).len(), + 1 + ); + + // MultiGeneUMI_All discards it from both: a UMI in two genes is + // evidence of a collision, not of the deeper gene. + assert!(filter_multi_gene_umi(&cross, UmiFiltering::MultiGeneUmiAll).is_empty()); + + // A single-gene UMI is untouched by every mode, including _All. + let mut single: HashMap = HashMap::default(); + single.insert(7, 3); + for mode in [ + UmiFiltering::None, + UmiFiltering::MultiGeneUmi, + UmiFiltering::MultiGeneUmiCr, + UmiFiltering::MultiGeneUmiAll, + ] { + assert_eq!( + filter_multi_gene_umi(&single, mode).len(), + 1, + "{mode:?} must not touch a single-gene UMI" + ); + } + } + + #[test] + fn multigene_umi_all_parses_to_its_own_variant() { + // It used to alias to MultiGeneUMI, which was neither STAR's behaviour + // (a no-op) nor the documented one. + assert_eq!( + "MultiGeneUMI_All".parse::().unwrap(), + UmiFiltering::MultiGeneUmiAll + ); + assert_eq!( + "MultiGeneUMI".parse::().unwrap(), + UmiFiltering::MultiGeneUmi + ); + } } From 4cb4744bbed07fcb6742a2b7adb03bd1fc251a48 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 30 Jul 2026 20:12:39 +0200 Subject: [PATCH 14/15] docs(changelog): record the MultiGeneUMI_All fix Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index af93658..087aa49 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -149,6 +149,11 @@ Sections commonly used: Features, Bug fixes, Other changes. Since one read per gene is the ordinary shape of a multi-gene UMI, the flag removed nothing in practice. On a 20k-read 10x fixture the count matrix moves from 16 465 to 15 414 against STAR's 15 423. +- `--soloUMIfiltering MultiGeneUMI_All` was aliased to `MultiGeneUMI`, + which is neither STAR's behaviour nor the documented one: in STAR + 2.7.11b the variant is a no-op. It now removes a UMI seen in two or + more genes from **all** of them, the behaviour the option name + describes. Recorded in `DIVERGENCE.md` (closes #144). - **STARsolo `Gene` assignment now requires exon concordance**, matching STARsolo: a read counts toward a gene only when every aligned block From d9d33b9fca40db7e26816ee4dd2db6b54b58764b Mon Sep 17 00:00:00 2001 From: Psy-Fer Date: Thu, 6 Aug 2026 20:35:34 +1000 Subject: [PATCH 15/15] docs(divergence): correct the MultiGeneUMI_All entry, defer the homopolymer one --- DIVERGENCE.md | 32 ++++++++++++++------------------ 1 file changed, 14 insertions(+), 18 deletions(-) diff --git a/DIVERGENCE.md b/DIVERGENCE.md index 6a87bb2..66f4088 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -32,31 +32,27 @@ This is the reason faithfulness is reported **tie-adjusted**. On the 10k yeast b ### 1.2 `--soloUMIfiltering MultiGeneUMI_All` filters, rather than doing nothing -**What STAR does.** The option is parsed and stored, but its consumption site tests only the `MultiGeneUMI` flag. Selecting `MultiGeneUMI_All` on its own therefore leaves the multi-gene UMI filter entirely off, and the counts are the unfiltered ones. STAR's own documentation describes it as removing a UMI seen in more than one gene from **all** of those genes. +**What STAR does.** Nothing, in effect — but not because the rule is unimplemented. `SoloFeature_collapseUMIall.cpp:79-88` implements exactly the documented behaviour, zeroing every gene for any UMI seen in more than one: -**What rustar-aligner does.** The documented behaviour: a UMI seen in more than one gene is removed from all of them. - -**Why.** Reproducing the no-op ships a flag that silently does nothing to anyone who read STAR's documentation. This was raised as #144 before any code changed, since "be faithful to STAR" and "do what the flag says" point in opposite directions here. Single-gene UMIs are untouched, which the tests check across every mode. - -**Impact.** Confined to `--soloUMIfiltering MultiGeneUMI_All`. The default (`-`) and the other filtering modes produce identical counts. Inverting the choice is a one-line change, since the test asserts the behaviour either way. - -**Source.** `src/solo/count.rs` (`UmiFiltering::MultiGeneUmiAll`, `filter_multi_gene_umi`), locked by `test_solo_multigene_umi_all_drops_cross_gene_umis`. STAR: `SoloFeature_collapseUMIall.cpp`, `ParametersSolo.cpp`. +```cpp +if (pSolo.umiFiltering.MultiGeneUMI_All) { + for (auto &iu : umiGeneMapCount) + if (iu.second.size()>1) + for (auto &ig : iu.second) ig.second=0; //kill all genes for this UMI +}; +``` -### 1.3 Homopolymer UMIs other than poly-A are rejected +The site that acts on those zeroed counts, however, gates on a different flag (`:116`, `if (pSolo.umiFiltering.MultiGeneUMI && umiGeneMapCount[...]==0)`), and `MultiGeneUMI` and `MultiGeneUMI_All` are set in mutually exclusive branches (`ParametersSolo.cpp:457-462`). Selecting `MultiGeneUMI_All` therefore zeroes the counts and then never reads them, and the run reports unfiltered counts. -**What STAR does.** The UMI validity check rejects a homopolymer by comparing the packed UMI against a precomputed all-same-base value, but the loop that builds those values runs over `umiL`, which is zero at that point for `CB_UMI_Simple`. Only the poly-A case (packed value zero) is caught; poly-C, poly-G and poly-T pass through as valid UMIs. - -**What rustar-aligner does.** Rejects every homopolymer UMI. - -**Why.** A homopolymer UMI is a sequencing artefact whatever base it repeats; letting three of the four through is not a rule, it is the consequence of reading an uninitialised length. +**What rustar-aligner does.** The documented behaviour: a UMI seen in more than one gene is removed from all of them. -**Impact.** Removes a small number of artefact UMIs from the counts that STAR keeps. +**Why.** This is a one-line wiring bug in STAR, not a design decision: STAR's own code, immediately above, computes the documented result and then discards it. Matching the binary would mean shipping a flag that silently does nothing to anyone who read either the documentation or STAR's source, which is what #144 was raised about. So the divergence is from STAR's behaviour but *not* from its intent. Single-gene UMIs are untouched, which the tests check across every mode. -**Source.** `src/solo/whitelist.rs` (`check_umi`), locked by `umi_valid_rejects_every_homopolymer`. STAR: `SoloReadBarcode_getCBandUMI.cpp` (`umiL`). +**Impact.** Confined to `--soloUMIfiltering MultiGeneUMI_All`. The default (`-`) and the other filtering modes produce identical counts. Inverting the choice is a one-line change, since the test asserts the behaviour either way. ---- +**Source.** `src/solo/count.rs` (`UmiFiltering::MultiGeneUmiAll`, `filter_multi_gene_umi`), locked by `multigene_umi_all_drops_the_umi_from_every_gene` and `multigene_umi_all_parses_to_its_own_variant`. STAR: `SoloFeature_collapseUMIall.cpp`, `ParametersSolo.cpp`. -### 1.2 `EmptyDrops_CR` Simple-Good-Turing with fewer than five distinct frequencies +### 1.3 `EmptyDrops_CR` Simple-Good-Turing with fewer than five distinct frequencies **What STAR does.** The ambient profile for `--soloCellFilter EmptyDrops_CR` is smoothed with Simple Good-Turing (Elworthy's `SimpleGoodTuring/sgt.h`). `analyse()` returns early, doing nothing, when the frequency spectrum has fewer than five distinct counts — Elworthy's `MinInput` guard. `PZero`, the mass reserved for genes unseen in the ambient droplets, is neither assigned in that case nor initialised at construction, so a caller that asks for it reads whatever the stack held.