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 diff --git a/DIVERGENCE.md b/DIVERGENCE.md index bd957ed..a4c1a44 100644 --- a/DIVERGENCE.md +++ b/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 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 + ); + } }