From 1383b87c78376f78c14c8673a6a637cae981e216 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 00:32:47 +0200 Subject: [PATCH 1/3] 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 bd957ed..a553e2d 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`). + --- ## 2. Cases where rustar-aligner outperforms STAR diff --git a/src/solo/count.rs b/src/solo/count.rs index 9af3edf..42d9e8e 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" @@ -822,6 +832,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 > 3.0: keep the highest-read-count gene(s); no singleton drop. UmiFiltering::MultiGeneUmiCr => genes.iter().filter(|&(_, &rc)| rc >= max).collect(), UmiFiltering::None => unreachable!(), @@ -2099,4 +2114,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 5d0168c7e3e2040341495967c3f379f4b6747aa1 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 11:46:32 +0200 Subject: [PATCH 2/3] docs(divergence): keep only this PR's divergence The homopolymer-UMI rule describes behaviour already on main, not something this diff changes, and its locking test does not exist yet. It belongs in its own change rather than riding along here. --- DIVERGENCE.md | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/DIVERGENCE.md b/DIVERGENCE.md index a553e2d..a4c1a44 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -42,18 +42,6 @@ This is the reason faithfulness is reported **tie-adjusted**. On the 10k yeast b **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`). - --- ## 2. Cases where rustar-aligner outperforms STAR From b7800b2df58cd0c258d42fbafbc36bf7f7b8ac13 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 30 Jul 2026 20:12:39 +0200 Subject: [PATCH 3/3] docs(changelog): record the MultiGeneUMI_All fix Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index b6e4893..cbe09e8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -99,6 +99,12 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Bug fixes +- `--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 lies within the gene's exons, rather than merely overlapping one. This