solo: --soloUMIfiltering MultiGeneUMI_All (closes #144)

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

DIVERGENCE.md

@@ -30,6 +30,18 @@ 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`.
+
 ---
 
 ## 2. Cases where rustar-aligner outperforms STAR

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<Self, Self::Err> {
         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<u32, u32>, 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<u32, u32> = 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<u32, u32> = 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::<UmiFiltering>().unwrap(),
+            UmiFiltering::MultiGeneUmiAll
+        );
+        assert_eq!(
+            "MultiGeneUMI".parse::<UmiFiltering>().unwrap(),
+            UmiFiltering::MultiGeneUmi
+        );
+    }
 }