From 51304b55078218f6efe258e6a1f649bb22c4b471 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 10:06:45 +0200 Subject: [PATCH 1/8] 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 9af3edf..7d906c8 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!(), } } @@ -2057,6 +2092,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 df80683ad2629da4f985d1c74d15ed232fc45a8e Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 10:06:45 +0200 Subject: [PATCH 2/8] docs(changelog): record the MultiGeneUMI_CR tie 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..91d4370 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -99,6 +99,12 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Bug fixes +- `--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 lies within the gene's exons, rather than merely overlapping one. This From a8f774e55c72c4dd8d30e155324575227744eb4a Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 10:42:46 +0200 Subject: [PATCH 3/8] feat(solo): --soloOutRawBarcodes Observed, for a CellRanger-shaped raw matrix MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit STARsolo's raw matrix has a column per whitelist barcode. For 10x v3 that is 3 686 400 columns and a 62 MB `barcodes.tsv`, nearly all zeros. CellRanger's `raw_feature_bc_matrix` has a column per *observed* barcode. The two files therefore share no keys, which is not a rounding difference in a comparison, it is zero overlap: comparing our raw output against a real `cellranger count` run gave 0 identical entries out of 27 396 until the columns were reconciled. `--soloOutRawBarcodes Observed` narrows the raw matrix to the barcodes that carry a count. Default `Whitelist` keeps what STARsolo writes, so nothing changes for anyone not asking. Measured on the 20 000-read fixture: Whitelist 3 686 400 barcodes barcodes.tsv 62 668 800 bytes Observed 200 barcodes barcodes.tsv 3 400 bytes with identical counts on both sides: 13 937 entries, 15 414 counts. `finalize_matrix` already took a column remap for the filtered matrix, so this reuses it rather than adding a second path. The observed set is read back from the streamed body, which costs one pass and only when the flag is on. This is a **non-STAR flag** and needs sign-off; recorded in `DIVERGENCE.md` §3.2 rather than presented as parity. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 6 +++++ DIVERGENCE.md | 27 +++++++++++++++++++++ src/params/mod.rs | 15 ++++++++++++ src/solo/count.rs | 62 +++++++++++++++++++++++++++++++++++++++++------ 4 files changed, 102 insertions(+), 8 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 91d4370..073d614 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,12 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Features +- `--soloOutRawBarcodes Observed` writes the raw matrix with one column + per *observed* barcode instead of one per whitelist barcode, matching + what CellRanger's `raw_feature_bc_matrix` contains. Counts are + unchanged; on a 200-cell run `barcodes.tsv` goes from 62 MB to 3.4 kB. + **Not a STAR parameter**; default `Whitelist` keeps STARsolo behaviour. + - **STARsolo single-cell quantification (`--soloType`)** — the 10x Chromium / plate-based count-matrix pipeline, ported from STAR and verified against real STARsolo (#90). diff --git a/DIVERGENCE.md b/DIVERGENCE.md index bd957ed..5f8e6e6 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -65,6 +65,33 @@ On the 10k yeast PE benchmark, 4 reads differ in alignment score (AS) because ST --- +### 3.2 `--soloOutRawBarcodes Observed` (opt-in, non-STAR) + +**What STAR does.** STARsolo's raw matrix has one column per whitelist +barcode, whether or not any read carried it. For 10x v3 that is 3 686 400 +columns and a 62 MB `barcodes.tsv`, nearly all of it zeros. + +**What rustar-aligner does.** The same, by default. `--soloOutRawBarcodes +Observed` narrows the raw matrix to the barcodes that actually hold a count, +which is what CellRanger's `raw_feature_bc_matrix` contains. On a 200-cell +fixture that is 200 columns and a 3.4 kB `barcodes.tsv`. + +**Why.** Someone comparing our raw matrix against CellRanger's finds no +overlapping keys at all, because the two files mean different things by "raw". +The flag makes the comparison possible without changing what STARsolo users +get. + +**Impact.** The counts are identical either way — same entries, same values, +verified on the fixture — only the columns present differ. This is a non-STAR +flag and needs maintainer sign-off; it is off by default so STARsolo parity is +untouched. + +**Source.** `src/solo/count.rs` (`observed_barcodes`), `src/params/mod.rs` +(`solo_out_raw_barcodes`). CellRanger: `outs/raw_feature_bc_matrix/` from a +`cellranger count` run, observed directly rather than taken from its source. + +--- + ## 4. Implementation divergences (no intended output difference) These differ in *how* a result is produced, not *what* is produced. They are documented so a reviewer chasing a discrepancy knows the mechanism differs by design. diff --git a/src/params/mod.rs b/src/params/mod.rs index 0536a85..9fc858f 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1133,6 +1133,21 @@ pub struct Parameters { #[arg(long = "soloOutGzip", default_value = "no")] pub solo_out_gzip: String, + /// Which barcodes the **raw** matrix has columns for. **Not a STAR + /// parameter**; a rustar-aligner addition, default `Whitelist`, which is + /// what STARsolo writes. + /// + /// `Whitelist` gives one column per whitelist barcode — 3.7 million of them + /// for 10x v3, whether or not a read ever carried them. `Observed` gives + /// one column per barcode that actually holds a count, which is what + /// CellRanger's `raw_feature_bc_matrix` contains, and turns a + /// hundreds-of-megabytes `barcodes.tsv` into a few kilobytes. + /// + /// The counts are identical either way; only the columns present differ. + #[arg(long = "soloOutRawBarcodes", default_value = "Whitelist", + value_parser = ["Whitelist", "Observed"])] + pub solo_out_raw_barcodes: String, + /// Velocyto ambiguous-molecule handling (rustar extension beyond STARsolo). /// `yes` (default) writes the three `spliced`/`unspliced`/`ambiguous` matrices /// like STARsolo — exon-only molecules with no junction/intron evidence stay in diff --git a/src/solo/count.rs b/src/solo/count.rs index 7d906c8..d766641 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -476,6 +476,27 @@ fn build_matrix_body( )) } +/// The whitelist indices that actually appear as a column in the streamed +/// matrix body, ascending. +/// +/// Reads the body once rather than tracking the set during counting, so the +/// default path pays nothing for a feature it does not use. +fn observed_barcodes(body: &tempfile::NamedTempFile) -> Result, Error> { + let reader = + BufReader::new(std::fs::File::open(body.path()).map_err(|e| Error::io(e, body.path()))?); + let mut seen: std::collections::BTreeSet = std::collections::BTreeSet::new(); + for line in reader.lines() { + let line = line.map_err(|e| Error::io(e, body.path()))?; + // " ", the layout `finalize_matrix` also parses. + if let Some(cb1) = line.split(' ').nth(1) + && let Ok(cb) = cb1.parse::() + { + seen.insert(cb.saturating_sub(1)); + } + } + Ok(seen.into_iter().collect()) +} + /// Write a final `matrix.mtx[.gz]` = MatrixMarket header + (optionally /// cb-remapped/filtered) body. With `remap = None` the body is copied verbatim /// (raw); with `Some(map)` only columns in the map survive, renumbered to the @@ -1235,20 +1256,45 @@ pub fn write_gene_matrix( &ctx.gene_ann.gene_names, gzip, )?; - write_barcodes( - &raw_dir.join(&barcodes_name), - &ctx.whitelist, - sorted.len(), - gzip, - )?; + // `--soloOutRawBarcodes Observed` narrows the raw matrix to the + // barcodes that actually carry a count, which is what CellRanger's + // `raw_feature_bc_matrix` holds. STARsolo's raw matrix has a column per + // whitelist barcode, so the default keeps that. + let observed: Option> = if params.solo_out_raw_barcodes == "Observed" { + Some(observed_barcodes(&body)?) + } else { + None + }; + let (raw_cols, raw_remap) = match &observed { + Some(cbs) => { + let map: HashMap = cbs + .iter() + .enumerate() + .map(|(col, &cb)| (cb, col as u32 + 1)) + .collect(); + (cbs.len(), Some(map)) + } + None => (sorted.len(), None), + }; + match &observed { + Some(cbs) => { + write_barcodes_subset(&raw_dir.join(&barcodes_name), &ctx.whitelist, cbs, gzip)?; + } + None => write_barcodes( + &raw_dir.join(&barcodes_name), + &ctx.whitelist, + sorted.len(), + gzip, + )?, + } finalize_matrix( &body, &raw_dir.join(&matrix_name), gzip, n_genes, - sorted.len(), + raw_cols, mstats.nnz, - None, + raw_remap.as_ref(), )?; log::info!( "STARsolo: wrote {}/raw matrix ({} genes × {} barcodes, {} entries){}", From abfcf3e2e39d346efde2868a6b87a5be110d9b50 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 10:53:46 +0200 Subject: [PATCH 4/8] 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 d766641..903514f 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; @@ -829,6 +849,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)> { @@ -2148,6 +2250,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 9f823e2bfd6a5152ae9949315c29212be9f47580 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 11:06:14 +0200 Subject: [PATCH 5/8] feat(solo): CellRanger behaviour by default on 10x geometry MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Aligning 10x data and comparing the count matrix against CellRanger gave a successful run and different numbers, with nothing in the output pointing at the five flags that explain the difference. Measured against CellRanger 10.0.0 on a 20 000-read fixture, those flags are the whole gap: 8.9% away without them, 0.03% with them. When the geometry is unambiguously 10x — `CB_UMI_Simple`, a whitelist, a 16-base cell barcode, a 10- or 12-base UMI — the five now default to their CellRanger values: --clipAdapterType CellRanger4 --outFilterScoreMin 30 --soloCBmatchWLtype 1MM_multi_Nbase_pseudocounts --soloUMIfiltering MultiGeneUMI_CR --soloUMIdedup 1MM_CR A flag given on the command line always wins, including when the value asked for is STARsolo's own default: `value_source` distinguishes an explicit flag from a default, so the divergence is escapable by naming what you want. Every substitution is logged at INFO with the geometry that triggered it. **This changes default output behaviour on 10x runs and diverges from STARsolo**, which is why it is confined to a geometry nothing else in common use shares, why it is announced on every run it touches, and why it is in `DIVERGENCE.md` §1.3 as the largest entry in that file. It needs sign-off. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 6 ++ DIVERGENCE.md | 28 +++++++ src/params/mod.rs | 208 ++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 242 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 073d614..45de81a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,12 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Features +- On 10x geometry (`CB_UMI_Simple`, a whitelist, 16 bp CB, 10 or 12 bp + UMI), the five CellRanger-matching flags now **default** to their + CellRanger values. Any flag named on the command line wins, and the + substitution is logged. **This changes default output on 10x runs** and + diverges from STARsolo; see `DIVERGENCE.md` §1.3. + - `--soloOutRawBarcodes Observed` writes the raw matrix with one column per *observed* barcode instead of one per whitelist barcode, matching what CellRanger's `raw_feature_bc_matrix` contains. Counts are diff --git a/DIVERGENCE.md b/DIVERGENCE.md index 5f8e6e6..8a0d92c 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -65,6 +65,34 @@ On the 10k yeast PE benchmark, 4 reads differ in alignment score (AS) because ST --- +### 1.3 CellRanger behaviour is the default on 10x geometry + +**What STAR does.** STARsolo's defaults are its own (`1MM_multi`, +`1MM_All`, no UMI filtering, `Hamming` clipping, `outFilterScoreMin 0`) +whatever the barcode geometry. Matching CellRanger requires passing five flags, +listed in STAR's `docs/STARsolo.md`. + +**What rustar-aligner does.** When the run is unambiguously 10x — +`CB_UMI_Simple`, a whitelist, a 16-base CB and a 10- or 12-base UMI — those +five flags default to their CellRanger values. Any flag given on the command +line wins, and the substitution is logged in full. + +**Why.** A user aligning 10x data and comparing against CellRanger otherwise +gets a successful run and different numbers, with nothing pointing at the five +flags that explain it. Measured against CellRanger 10.0.0 on a 20 000-read +fixture, those flags move the count matrix from 8.9% away to 0.03%. + +**Impact.** This is a **change of default output behaviour** and therefore the +largest divergence in this file. It is confined to a geometry nothing else in +common use shares, it is escapable by naming any flag explicitly, and it is +announced at `INFO` on every run it touches. It needs maintainer sign-off. + +**Source.** `src/params/mod.rs` (`looks_like_10x`, +`apply_cellranger_defaults_on_10x`). STAR: `docs/STARsolo.md`, "Matching +CellRanger 4.x and 5.x results". + +--- + ### 3.2 `--soloOutRawBarcodes Observed` (opt-in, non-STAR) **What STAR does.** STARsolo's raw matrix has one column per whitelist diff --git a/src/params/mod.rs b/src/params/mod.rs index 9fc858f..4ace622 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1363,6 +1363,8 @@ impl Parameters { let matches = command.clone().get_matches_from(args.iter()); let mut params = ::from_arg_matches(&matches)?; + apply_cellranger_defaults_on_10x(&mut params, &matches); + params.command_line = { let args: Vec<_> = args.iter().map(|s| s.to_string_lossy()).collect(); shlex::try_join(args.iter().map(AsRef::as_ref)).ok() @@ -1885,8 +1887,214 @@ impl Parameters { // Tests // --------------------------------------------------------------------------- +/// The flags STAR documents for matching CellRanger 4.x/5.x +/// (`docs/STARsolo.md`), applied by default when the run is a 10x one. +const CELLRANGER_DEFAULTS: [(&str, &str); 5] = [ + ("clip_adapter_type", "CellRanger4"), + ("out_filter_score_min", "30"), + ("solo_cb_match_wl_type", "1MM_multi_Nbase_pseudocounts"), + ("solo_umi_filtering", "MultiGeneUMI_CR"), + ("solo_umi_dedup", "1MM_CR"), +]; + +/// Does this look like a 10x Chromium run? +/// +/// `CB_UMI_Simple` with a whitelist, a 16-base cell barcode, and a UMI of 10 +/// (v2) or 12 (v3) bases. That is the geometry of every 10x 3'/5' gene +/// expression chemistry, and nothing else in common use shares it. +fn looks_like_10x(params: &Parameters) -> bool { + params.solo_type == SoloType::CbUmiSimple + && params.solo_cb_len == 16 + && (params.solo_umi_len == 10 || params.solo_umi_len == 12) + && params + .solo_cb_whitelist + .first() + .is_some_and(|w| w != "None" && w != "-") +} + +/// On a 10x run, default to CellRanger's behaviour rather than STARsolo's. +/// +/// **This diverges from STAR by default**, which is why it is confined to a +/// geometry that is unambiguously 10x, and why every flag it changes is +/// logged. A flag given on the command line always wins, so the change is +/// invisible to anyone who states what they want. +/// +/// The rationale is that a user aligning 10x data and comparing against +/// CellRanger currently gets a successful run and different numbers, with +/// nothing pointing at the five flags that explain the difference. Measured on +/// a 20 000-read fixture, those flags move the count matrix from 8.9% away +/// from CellRanger to 0.03%. +/// +/// Recorded in `DIVERGENCE.md`; it needs maintainer sign-off. +fn apply_cellranger_defaults_on_10x(params: &mut Parameters, matches: &clap::ArgMatches) { + use clap::parser::ValueSource; + + if !looks_like_10x(params) { + return; + } + + let given = |id: &str| matches.value_source(id) == Some(ValueSource::CommandLine); + + let mut applied: Vec<&str> = Vec::new(); + for (id, value) in CELLRANGER_DEFAULTS { + if given(id) { + continue; + } + match id { + "clip_adapter_type" => params.clip_adapter_type = value.to_string(), + "out_filter_score_min" => params.out_filter_score_min = 30, + "solo_cb_match_wl_type" => params.solo_cb_match_wl_type = value.to_string(), + "solo_umi_filtering" => params.solo_umi_filtering = vec![value.to_string()], + "solo_umi_dedup" => params.solo_umi_dedup = vec![value.to_string()], + _ => continue, + } + applied.push(value); + } + + if !applied.is_empty() { + log::info!( + "10x geometry detected (CB {} + UMI {} with a whitelist): defaulting to \ + CellRanger behaviour [{}]. Pass the flags explicitly to override; this \ + differs from STARsolo's defaults.", + params.solo_cb_len, + params.solo_umi_len, + applied.join(", ") + ); + } +} + #[cfg(test)] mod tests { + + /// 10x geometry with a whitelist gets CellRanger's five flags without the + /// user naming any of them. This is a deliberate divergence from STARsolo's + /// defaults, so the test states the whole set rather than spot-checking one. + #[test] + fn ten_x_geometry_defaults_to_cellranger_behaviour() { + let p = Parameters::try_parse_from([ + "rustar-aligner", + "--readFilesIn", + "cdna.fq", + "cb.fq", + "--sjdbGTFfile", + "genes.gtf", + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + "wl.txt", + "--soloCBstart", + "1", + "--soloCBlen", + "16", + "--soloUMIstart", + "17", + "--soloUMIlen", + "12", + ]) + .unwrap(); + assert_eq!(p.clip_adapter_type, "CellRanger4"); + assert_eq!(p.out_filter_score_min, 30); + assert_eq!(p.solo_cb_match_wl_type, "1MM_multi_Nbase_pseudocounts"); + assert_eq!(p.solo_umi_filtering, vec!["MultiGeneUMI_CR".to_string()]); + assert_eq!(p.solo_umi_dedup, vec!["1MM_CR".to_string()]); + } + + /// A flag given on the command line always wins, including when the value + /// asked for is STARsolo's own default. Without this the divergence would + /// be inescapable, which is a different and much worse thing than a + /// divergent default. + #[test] + fn an_explicit_flag_beats_the_10x_default() { + let p = Parameters::try_parse_from([ + "rustar-aligner", + "--readFilesIn", + "cdna.fq", + "cb.fq", + "--sjdbGTFfile", + "genes.gtf", + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + "wl.txt", + "--soloCBstart", + "1", + "--soloCBlen", + "16", + "--soloUMIstart", + "17", + "--soloUMIlen", + "12", + "--soloUMIdedup", + "1MM_All", + "--clipAdapterType", + "Hamming", + ]) + .unwrap(); + assert_eq!(p.solo_umi_dedup, vec!["1MM_All".to_string()]); + assert_eq!(p.clip_adapter_type, "Hamming"); + // The ones not named still take the CellRanger value. + assert_eq!(p.out_filter_score_min, 30); + } + + /// Geometry that is not 10x is left alone: a 12-base barcode is not any + /// Chromium chemistry, so nothing is overridden. + #[test] + fn non_10x_geometry_keeps_starsolo_defaults() { + let p = Parameters::try_parse_from([ + "rustar-aligner", + "--readFilesIn", + "cdna.fq", + "cb.fq", + "--sjdbGTFfile", + "genes.gtf", + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + "wl.txt", + "--soloCBstart", + "1", + "--soloCBlen", + "12", + "--soloUMIstart", + "13", + "--soloUMIlen", + "8", + ]) + .unwrap(); + assert_eq!(p.clip_adapter_type, "Hamming"); + assert_eq!(p.out_filter_score_min, 0); + assert_eq!(p.solo_cb_match_wl_type, "1MM_multi"); + } + + /// No whitelist means no 10x run, whatever the lengths say. (Without a + /// whitelist the CB-match type must be Exact anyway, which is unrelated + /// validation that predates this and is stated here so the test reads.) + #[test] + fn ten_x_lengths_without_a_whitelist_keep_starsolo_defaults() { + let p = Parameters::try_parse_from([ + "rustar-aligner", + "--readFilesIn", + "cdna.fq", + "cb.fq", + "--sjdbGTFfile", + "genes.gtf", + "--soloType", + "CB_UMI_Simple", + "--soloCBstart", + "1", + "--soloCBlen", + "16", + "--soloUMIstart", + "17", + "--soloUMIlen", + "12", + "--soloCBmatchWLtype", + "Exact", + ]) + .unwrap(); + assert_eq!(p.clip_adapter_type, "Hamming"); + assert_eq!(p.out_filter_score_min, 0); + } use super::*; /// Helper: parse a STAR-style command line (without program name). From 014d3a7929133e5c6d7cba8ae9c266defb5008c8 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 11:48:48 +0200 Subject: [PATCH 6/8] docs(divergence): correct the CellRanger gap figure in 1.3 Re-derived both sides from one clean state: the five flags move the matrix from 8.96% above CellRanger to 2.17% above it, not to 0.03%. The earlier figure compared a rustar run against a CellRanger run built from a different state of the fixture. STAR 2.7.11b with the same flags is at +0.09%, so the remaining gap is open rather than closed. Co-Authored-By: Claude Opus 5 (1M context) --- DIVERGENCE.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/DIVERGENCE.md b/DIVERGENCE.md index 8a0d92c..592aa1d 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -80,7 +80,9 @@ line wins, and the substitution is logged in full. **Why.** A user aligning 10x data and comparing against CellRanger otherwise gets a successful run and different numbers, with nothing pointing at the five flags that explain it. Measured against CellRanger 10.0.0 on a 20 000-read -fixture, those flags move the count matrix from 8.9% away to 0.03%. +fixture, those flags move the count matrix from 8.96% above CellRanger to +2.17% above it. The remaining 2.17% is an open divergence: STAR 2.7.11b with +the same flags is at +0.09%, so this closes most of the gap and not all of it. **Impact.** This is a **change of default output behaviour** and therefore the largest divergence in this file. It is confined to a geometry nothing else in From 925da7e8e5d35f05270d80614131f711bb337b15 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 14:21:10 +0200 Subject: [PATCH 7/8] docs(divergence): restore the measured CellRanger gap in 1.3 The figure I replaced this with was measured without #165, whose cbMinP posterior threshold is a precondition for it. With #165 the five flags move the matrix from 8.96% above CellRanger to 0.03% above it, which is what the original text said. Co-Authored-By: Claude Opus 5 (1M context) --- DIVERGENCE.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/DIVERGENCE.md b/DIVERGENCE.md index 592aa1d..60326d9 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -81,8 +81,9 @@ line wins, and the substitution is logged in full. gets a successful run and different numbers, with nothing pointing at the five flags that explain it. Measured against CellRanger 10.0.0 on a 20 000-read fixture, those flags move the count matrix from 8.96% above CellRanger to -2.17% above it. The remaining 2.17% is an open divergence: STAR 2.7.11b with -the same flags is at +0.09%, so this closes most of the gap and not all of it. +0.03% above it, once #165's `cbMinP` posterior threshold is also applied. +STAR 2.7.11b with the same flags is at +0.09%, so all three agree to within a +fraction of a percent. **Impact.** This is a **change of default output behaviour** and therefore the largest divergence in this file. It is confined to a geometry nothing else in From df2253a45f43dee6395f587d0bf762b70041c1bb Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 11:27:58 +0200 Subject: [PATCH 8/8] feat(solo): --soloOutLayout CellRanger, CellRanger's output shape Writes the same numbers where `cellranger count` writes them: outs/{raw,filtered}_feature_bc_matrix/, gzipped, a -1 GEM-well suffix on every barcode, one raw column per observed barcode, and no per-feature subdirectory when a single feature is requested. The layout implies --soloOutGzip yes, --soloOutRawBarcodes Observed and an outs/ output directory; each is still overridable on the command line. On 10x geometry it is the default, alongside the five flags from the previous commit, so it changes where output files are written on those runs. Counts are untouched: 13 959 entries and 15 439 counts in both layouts on the 20 000-read fixture, compared entry by entry. Against a real cellranger count 10.0.0 run on the same fixture the raw barcode sets match exactly, 200 of 200. Recorded in DIVERGENCE.md 3.3. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 9 ++ Cargo.toml | 3 + DIVERGENCE.md | 35 ++++++++ src/params/mod.rs | 163 +++++++++++++++++++++++++++++++++++- src/solo/count.rs | 103 +++++++++++++++++++---- tests/alignment_features.rs | 151 +++++++++++++++++++++++++++++++++ 6 files changed, 447 insertions(+), 17 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 45de81a..a8ba71b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,15 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Features +- `--soloOutLayout CellRanger` writes the solo matrices in the shape + `cellranger count` produces: `outs/raw_feature_bc_matrix/` and + `outs/filtered_feature_bc_matrix/`, gzipped, with a `-1` GEM-well suffix + on every barcode and one raw column per observed barcode. It implies + `--soloOutGzip yes`, `--soloOutRawBarcodes Observed` and an `outs/` + output directory, each still overridable on the command line. Counts are + unchanged. It is the default on 10x geometry, which **changes where + output files are written** on those runs; see `DIVERGENCE.md` §3.3. + - On 10x geometry (`CB_UMI_Simple`, a whitelist, 16 bp CB, 10 or 12 bp UMI), the five CellRanger-matching flags now **default** to their CellRanger values. Any flag named on the command line wins, and the diff --git a/Cargo.toml b/Cargo.toml index 8a4638f..444553d 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -68,6 +68,9 @@ noodles-bgzf = { version = "0.49", features = ["libdeflate"] } [dev-dependencies] assert_cmd = "2" +# Already a runtime dependency at the same version; listed here so the +# integration tests can read the gzipped CellRanger-layout matrix files. +flate2 = { version = "1", default-features = false, features = ["zlib-rs"] } predicates = "3" [build-dependencies] diff --git a/DIVERGENCE.md b/DIVERGENCE.md index 60326d9..79248cb 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -123,6 +123,41 @@ untouched. --- +### 3.3 `--soloOutLayout CellRanger` (non-STAR; on by default on 10x geometry) + +**What STAR does.** STARsolo writes +`Solo.out//{raw,filtered}/{matrix.mtx, barcodes.tsv, features.tsv}`, +uncompressed, with bare barcodes. + +**What rustar-aligner does.** The same, by default, on non-10x geometry. +`--soloOutLayout CellRanger` writes the same numbers in the shape +`cellranger count` produces: `outs/{raw,filtered}_feature_bc_matrix/`, all three +files gzipped, a `-1` GEM-well suffix on every barcode, one raw column per +observed barcode, and no per-feature subdirectory when a single feature is +requested. It implies `--soloOutGzip yes`, `--soloOutRawBarcodes Observed` and +`--soloOutFileNames outs/ ...`, each still overridable on the command line. + +Under §1.3 the same 10x detection turns this on by default, so a bare 10x run +lands in CellRanger's layout. + +**Why.** Tools written against CellRanger's `outs/` (scanpy's `read_10x_mtx`, +Seurat's `Read10X`, any in-house loader) key on those directory names and on +the `-1` suffix. Without them the numbers are right and nothing downstream can +read them without a rename step. + +**Impact.** No count changes: verified entry by entry on the 20 000-read +fixture, 13 959 entries and 15 439 counts in both layouts. On 10x geometry it +**changes where output files are written**, which needs maintainer sign-off +alongside §1.3. Against a real `cellranger count` run on the same fixture the +raw barcode sets match exactly, 200 of 200. + +**Source.** `src/solo/count.rs` (`write_gene_matrix`, `write_one_barcode`), +`src/params/mod.rs` (`solo_out_layout`, `apply_cellranger_layout`). CellRanger: +`outs/` from a `cellranger count` 10.0.0 run, observed directly rather than +taken from its source. + +--- + ## 4. Implementation divergences (no intended output difference) These differ in *how* a result is produced, not *what* is produced. They are documented so a reviewer chasing a discrepancy knows the mechanism differs by design. diff --git a/src/params/mod.rs b/src/params/mod.rs index 4ace622..edfa3c9 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1148,6 +1148,29 @@ pub struct Parameters { value_parser = ["Whitelist", "Observed"])] pub solo_out_raw_barcodes: String, + /// Shape of the solo matrix output on disk. **Not a STAR parameter**; a + /// rustar-aligner addition, default `STARsolo`, which changes nothing. + /// + /// `CellRanger` lays the same numbers out the way `cellranger count` does, + /// so a tool written against CellRanger's `outs/` reads a rustar run + /// unmodified. It implies, unless the corresponding flag is given + /// explicitly: + /// + /// * directories `raw_feature_bc_matrix/` and `filtered_feature_bc_matrix/` + /// instead of `raw/` and `filtered/`; + /// * a `-1` GEM-well suffix on every barcode; + /// * gzip on all three files (`--soloOutGzip yes`); + /// * `--soloOutRawBarcodes Observed`, since CellRanger's raw matrix has one + /// column per observed barcode, not one per whitelist entry; + /// * the output directory `outs/` instead of `Solo.out/`, with no + /// per-feature subdirectory when exactly one feature is requested. + /// + /// Counts are untouched. Only where the bytes land, and how the barcodes + /// are spelled, changes. + #[arg(long = "soloOutLayout", default_value = "STARsolo", + value_parser = ["STARsolo", "CellRanger"])] + pub solo_out_layout: String, + /// Velocyto ambiguous-molecule handling (rustar extension beyond STARsolo). /// `yes` (default) writes the three `spliced`/`unspliced`/`ambiguous` matrices /// like STARsolo — exon-only molecules with no junction/intron evidence stay in @@ -1364,6 +1387,7 @@ impl Parameters { let mut params = ::from_arg_matches(&matches)?; apply_cellranger_defaults_on_10x(&mut params, &matches); + apply_cellranger_layout(&mut params, &matches); params.command_line = { let args: Vec<_> = args.iter().map(|s| s.to_string_lossy()).collect(); @@ -1889,12 +1913,15 @@ impl Parameters { /// The flags STAR documents for matching CellRanger 4.x/5.x /// (`docs/STARsolo.md`), applied by default when the run is a 10x one. -const CELLRANGER_DEFAULTS: [(&str, &str); 5] = [ +const CELLRANGER_DEFAULTS: [(&str, &str); 6] = [ ("clip_adapter_type", "CellRanger4"), ("out_filter_score_min", "30"), ("solo_cb_match_wl_type", "1MM_multi_Nbase_pseudocounts"), ("solo_umi_filtering", "MultiGeneUMI_CR"), ("solo_umi_dedup", "1MM_CR"), + // Not a STAR flag: the output layout, so a 10x run lands where a tool + // written against `cellranger count` expects to find it. + ("solo_out_layout", "CellRanger"), ]; /// Does this look like a 10x Chromium run? @@ -1946,6 +1973,7 @@ fn apply_cellranger_defaults_on_10x(params: &mut Parameters, matches: &clap::Arg "solo_cb_match_wl_type" => params.solo_cb_match_wl_type = value.to_string(), "solo_umi_filtering" => params.solo_umi_filtering = vec![value.to_string()], "solo_umi_dedup" => params.solo_umi_dedup = vec![value.to_string()], + "solo_out_layout" => params.solo_out_layout = value.to_string(), _ => continue, } applied.push(value); @@ -1963,6 +1991,34 @@ fn apply_cellranger_defaults_on_10x(params: &mut Parameters, matches: &clap::Arg } } +/// `--soloOutLayout CellRanger` implies three existing flags and the output +/// directory name. Each is still overridable: a flag named on the command line +/// keeps its value, so the layout can be adopted piecemeal. +/// +/// Runs after `apply_cellranger_defaults_on_10x`, so it sees the layout whether +/// it was asked for or inferred from 10x geometry. +fn apply_cellranger_layout(params: &mut Parameters, matches: &clap::ArgMatches) { + use clap::parser::ValueSource; + + if params.solo_out_layout != "CellRanger" { + return; + } + let given = |id: &str| matches.value_source(id) == Some(ValueSource::CommandLine); + + if !given("solo_out_gzip") { + params.solo_out_gzip = "yes".to_string(); + } + if !given("solo_out_raw_barcodes") { + params.solo_out_raw_barcodes = "Observed".to_string(); + } + if !given("solo_out_file_names") + && let Some(dir) = params.solo_out_file_names.first_mut() + { + // CellRanger writes its matrices under `outs/`, not `Solo.out/`. + *dir = "outs/".to_string(); + } +} + #[cfg(test)] mod tests { @@ -1997,6 +2053,111 @@ mod tests { assert_eq!(p.solo_cb_match_wl_type, "1MM_multi_Nbase_pseudocounts"); assert_eq!(p.solo_umi_filtering, vec!["MultiGeneUMI_CR".to_string()]); assert_eq!(p.solo_umi_dedup, vec!["1MM_CR".to_string()]); + assert_eq!(p.solo_out_layout, "CellRanger"); + } + + /// `--soloOutLayout CellRanger` pulls three existing flags and the output + /// directory with it, so the layout is one decision rather than four. + #[test] + fn cellranger_layout_implies_gzip_observed_barcodes_and_outs_dir() { + let p = Parameters::try_parse_from([ + "rustar-aligner", + "--readFilesIn", + "cdna.fq", + "cb.fq", + "--sjdbGTFfile", + "genes.gtf", + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + "wl.txt", + "--soloCBstart", + "1", + "--soloCBlen", + "12", + "--soloUMIstart", + "13", + "--soloUMIlen", + "10", + "--soloOutLayout", + "CellRanger", + ]) + .unwrap(); + // 12-base CB, so the 10x autodetection is not what set this. + assert_eq!(p.solo_out_layout, "CellRanger"); + assert_eq!(p.solo_out_gzip, "yes"); + assert_eq!(p.solo_out_raw_barcodes, "Observed"); + assert_eq!(p.solo_out_file_names.first().unwrap(), "outs/"); + } + + /// The layout is adoptable piecemeal: each flag it implies is still + /// overridable on the command line. + #[test] + fn an_explicit_flag_beats_what_the_cellranger_layout_implies() { + let p = Parameters::try_parse_from([ + "rustar-aligner", + "--readFilesIn", + "cdna.fq", + "cb.fq", + "--sjdbGTFfile", + "genes.gtf", + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + "wl.txt", + "--soloCBstart", + "1", + "--soloCBlen", + "16", + "--soloUMIstart", + "17", + "--soloUMIlen", + "12", + "--soloOutGzip", + "no", + "--soloOutRawBarcodes", + "Whitelist", + "--soloOutFileNames", + "Solo.out/", + "features.tsv", + "barcodes.tsv", + "matrix.mtx", + ]) + .unwrap(); + assert_eq!(p.solo_out_layout, "CellRanger"); + assert_eq!(p.solo_out_gzip, "no"); + assert_eq!(p.solo_out_raw_barcodes, "Whitelist"); + assert_eq!(p.solo_out_file_names.first().unwrap(), "Solo.out/"); + } + + /// The default is STARsolo's layout: no `-1`, no `outs/`, no gzip. + #[test] + fn non_10x_geometry_keeps_the_starsolo_layout() { + let p = Parameters::try_parse_from([ + "rustar-aligner", + "--readFilesIn", + "cdna.fq", + "cb.fq", + "--sjdbGTFfile", + "genes.gtf", + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + "wl.txt", + "--soloCBstart", + "1", + "--soloCBlen", + "12", + "--soloUMIstart", + "13", + "--soloUMIlen", + "10", + ]) + .unwrap(); + assert_eq!(p.solo_out_layout, "STARsolo"); + assert_eq!(p.solo_out_gzip, "no"); + assert_eq!(p.solo_out_raw_barcodes, "Whitelist"); + assert_eq!(p.solo_out_file_names.first().unwrap(), "Solo.out/"); } /// A flag given on the command line always wins, including when the value diff --git a/src/solo/count.rs b/src/solo/count.rs index 903514f..613af35 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -1334,10 +1334,28 @@ pub fn write_gene_matrix( let n_genes = ctx.gene_ann.gene_ids.len(); let multi_methods = MultiMethod::parse_list(¶ms.solo_multi_mappers); + // `--soloOutLayout CellRanger`: CellRanger's directory names, and a `-1` + // GEM-well suffix on every barcode. The counts are the same either way. + let cr_layout = params.solo_out_layout == "CellRanger"; + let (raw_name, filt_name) = if cr_layout { + ("raw_feature_bc_matrix", "filtered_feature_bc_matrix") + } else { + ("raw", "filtered") + }; + let cb_suffix = if cr_layout { "-1" } else { "" }; + // CellRanger has no per-feature directory. Drop ours when there is exactly + // one feature; keep it when there are several, because collapsing them + // would have each feature silently overwrite the last. + let per_feature_dir = !cr_layout || ctx.features.len() > 1; + // One {prefix}{soloOutFileNames[0]}/{raw,filtered}/ per feature. for (feature, recorder) in ctx.features.iter().zip(&ctx.recorders) { - let feature_dir = params.output_path(&format!("{solo_dir}{}/", feature.dir_name())); - let raw_dir = feature_dir.join("raw"); + let feature_dir = if per_feature_dir { + params.output_path(&format!("{solo_dir}{}/", feature.dir_name())) + } else { + params.output_path(&solo_dir) + }; + let raw_dir = feature_dir.join(raw_name); std::fs::create_dir_all(&raw_dir).map_err(|e| Error::io(e, &raw_dir))?; // Stream the deduplicated counts into a shared temp body, then finalize @@ -1380,13 +1398,20 @@ pub fn write_gene_matrix( }; match &observed { Some(cbs) => { - write_barcodes_subset(&raw_dir.join(&barcodes_name), &ctx.whitelist, cbs, gzip)?; + write_barcodes_subset( + &raw_dir.join(&barcodes_name), + &ctx.whitelist, + cbs, + gzip, + cb_suffix, + )?; } None => write_barcodes( &raw_dir.join(&barcodes_name), &ctx.whitelist, sorted.len(), gzip, + cb_suffix, )?, } finalize_matrix( @@ -1399,10 +1424,10 @@ pub fn write_gene_matrix( raw_remap.as_ref(), )?; log::info!( - "STARsolo: wrote {}/raw matrix ({} genes × {} barcodes, {} entries){}", - feature.dir_name(), + "STARsolo: wrote {} matrix ({} genes × {} barcodes, {} entries){}", + raw_dir.display(), n_genes, - sorted.len(), + raw_cols, mstats.nnz, if gzip { " [gzip]" } else { "" }, ); @@ -1426,7 +1451,7 @@ pub fn write_gene_matrix( if let Some(cbs) = called && !cbs.is_empty() { - let filt_dir = feature_dir.join("filtered"); + let filt_dir = feature_dir.join(filt_name); std::fs::create_dir_all(&filt_dir).map_err(|e| Error::io(e, &filt_dir))?; let remap: HashMap = cbs .iter() @@ -1439,7 +1464,13 @@ pub fn write_gene_matrix( &ctx.gene_ann.gene_names, gzip, )?; - write_barcodes_subset(&filt_dir.join(&barcodes_name), &ctx.whitelist, &cbs, gzip)?; + write_barcodes_subset( + &filt_dir.join(&barcodes_name), + &ctx.whitelist, + &cbs, + gzip, + cb_suffix, + )?; let fnnz = finalize_matrix( &body, &filt_dir.join(&matrix_name), @@ -1450,8 +1481,8 @@ pub fn write_gene_matrix( Some(&remap), )?; log::info!( - "STARsolo: wrote {}/filtered matrix ({} cells, {} entries)", - feature.dir_name(), + "STARsolo: wrote {} matrix ({} cells, {} entries)", + filt_dir.display(), cbs.len(), fnnz, ); @@ -1522,11 +1553,14 @@ pub fn write_gene_matrix( write_file(&sj_dir.join(&features_name), gzip, |w| { sjs.write_sj_lines(w, genome, params).map(|_| ()) })?; + // SJ has no CellRanger counterpart, but every barcode written by one run + // is spelled the same way, so the suffix applies here too. write_barcodes( &sj_dir.join(&barcodes_name), &ctx.whitelist, sorted.len(), gzip, + cb_suffix, )?; let umi_len = params.solo_umi_len as usize; let nnz = build_sj_matrix( @@ -1562,6 +1596,7 @@ pub fn write_gene_matrix( &ctx.whitelist, sorted.len(), gzip, + cb_suffix, )?; let umi_len = params.solo_umi_len as usize; // `--soloVelocytoAmbiguous no` folds exon-only molecules into spliced and @@ -2016,16 +2051,21 @@ fn write_features( Ok(()) } -/// Unpack `cb` into `line` (with trailing newline) and write it. +/// Unpack `cb` into `line` (with `suffix` and a trailing newline) and write it. +/// +/// `suffix` is `"-1"` under `--soloOutLayout CellRanger` (the GEM-well tag +/// CellRanger appends to every barcode) and `""` otherwise. fn write_one_barcode( w: &mut dyn std::io::Write, whitelist: &CbWhitelist, cb: u32, line: &mut Vec, path: &Path, + suffix: &str, ) -> Result<(), Error> { line.clear(); whitelist.unpack_barcode_into(cb, line); + line.extend_from_slice(suffix.as_bytes()); line.push(b'\n'); w.write_all(line).map_err(|e| Error::io(e, path)) } @@ -2033,12 +2073,18 @@ fn write_one_barcode( /// `barcodes.tsv`: full whitelist in sorted order (matches the raw matrix /// columns). Lists millions of lines, so the writer is buffered and the barcode /// is unpacked into a reused scratch buffer (no per-line allocation). -fn write_barcodes(path: &Path, whitelist: &CbWhitelist, n: usize, gzip: bool) -> Result<(), Error> { +fn write_barcodes( + path: &Path, + whitelist: &CbWhitelist, + n: usize, + gzip: bool, + suffix: &str, +) -> Result<(), Error> { let len = whitelist.barcode_len(); write_file(path, gzip, |w| { - let mut line: Vec = Vec::with_capacity(len + 1); + let mut line: Vec = Vec::with_capacity(len + suffix.len() + 1); for i in 0..n { - write_one_barcode(w, whitelist, i as u32, &mut line, path)?; + write_one_barcode(w, whitelist, i as u32, &mut line, path, suffix)?; } Ok(()) })?; @@ -2052,12 +2098,13 @@ fn write_barcodes_subset( whitelist: &CbWhitelist, cbs: &[u32], gzip: bool, + suffix: &str, ) -> Result<(), Error> { let len = whitelist.barcode_len(); write_file(path, gzip, |w| { - let mut line: Vec = Vec::with_capacity(len + 1); + let mut line: Vec = Vec::with_capacity(len + suffix.len() + 1); for &cb in cbs { - write_one_barcode(w, whitelist, cb, &mut line, path)?; + write_one_barcode(w, whitelist, cb, &mut line, path, suffix)?; } Ok(()) })?; @@ -2070,6 +2117,30 @@ mod tests { use crate::io::fastq::encode_base; use crate::solo::whitelist::pack_barcode; + /// `--soloOutLayout CellRanger` appends the `-1` GEM-well suffix that + /// CellRanger puts on every barcode; the default writes the bare barcode. + #[test] + fn barcodes_carry_the_gem_well_suffix_only_under_the_cellranger_layout() { + let dir = tempfile::tempdir().unwrap(); + let wl_path = dir.path().join("wl.txt"); + std::fs::write(&wl_path, "ACGT\nTGCA\n").unwrap(); + let wl = CbWhitelist::load(&wl_path).unwrap(); + + let plain = dir.path().join("plain.tsv"); + write_barcodes(&plain, &wl, wl.len(), false, "").unwrap(); + assert_eq!(std::fs::read_to_string(&plain).unwrap(), "ACGT\nTGCA\n"); + + let cr = dir.path().join("cr.tsv"); + write_barcodes(&cr, &wl, wl.len(), false, "-1").unwrap(); + assert_eq!(std::fs::read_to_string(&cr).unwrap(), "ACGT-1\nTGCA-1\n"); + + // The subset writer (filtered matrix, and the raw one under + // `--soloOutRawBarcodes Observed`) takes the same suffix. + let sub = dir.path().join("sub.tsv"); + write_barcodes_subset(&sub, &wl, &[1], false, "-1").unwrap(); + assert_eq!(std::fs::read_to_string(&sub).unwrap(), "TGCA-1\n"); + } + #[test] fn median_sorted_odd_even_empty() { assert_eq!(median_sorted(&[]), 0); diff --git a/tests/alignment_features.rs b/tests/alignment_features.rs index 2db7d39..a0f812d 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -951,6 +951,11 @@ fn test_starsolo_gene_matrix() { "Gene", "--sjdbGTFfile", gtf.to_str().unwrap(), + // This fixture's 16 bp CB + 12 bp UMI is 10x geometry, which now + // defaults to CellRanger's output layout; the assertions below are + // about STARsolo's, so state it. + "--soloOutLayout", + "STARsolo", "--outFileNamePrefix", &prefix, ]) @@ -1028,6 +1033,122 @@ fn test_starsolo_gene_matrix() { ); } +// --------------------------------------------------------------------------- +// Test 9a'' — --soloOutLayout CellRanger writes the same numbers where +// `cellranger count` writes them: outs/{raw,filtered}_feature_bc_matrix/, +// gzipped, with a -1 GEM-well suffix on every barcode. +// --------------------------------------------------------------------------- + +#[test] +fn test_solo_out_layout_cellranger() { + use std::io::Read; + + let tmpdir = TempDir::new().unwrap(); + let genome = build_genome(); + let fasta = write_fasta(&tmpdir, &genome); + let gtf = write_gtf(&tmpdir); + + let genome_dir = tmpdir.path().join("genome"); + build_index(&fasta, &genome_dir, "7", Some(>f)); + + // Same fixture as test_starsolo_gene_matrix: 8 reads, one cell, two UMI + // clouds, so the expected counts are already known. + let cdna_path = tmpdir.path().join("cdna.fq"); + let barcode_path = tmpdir.path().join("barcode.fq"); + let wl_path = tmpdir.path().join("whitelist.txt"); + let cb = "AAAACCCCGGGGTTTT"; + let umi_a = "ACGTACGTAC"; + let umi_b = "TGCATGCATG"; + { + let mut cf = fs::File::create(&cdna_path).unwrap(); + let mut bf = fs::File::create(&barcode_path).unwrap(); + let exon1 = &genome[10000..10050]; + for i in 0..8usize { + writeln!(cf, "@read{i}").unwrap(); + cf.write_all(exon1).unwrap(); + writeln!(cf, "\n+\n{}", "I".repeat(50)).unwrap(); + let umi = if i < 4 { umi_a } else { umi_b }; + writeln!(bf, "@read{i}").unwrap(); + writeln!(bf, "{cb}{umi}").unwrap(); + writeln!(bf, "+\n{}", "I".repeat(26)).unwrap(); + } + } + { + let mut wf = fs::File::create(&wl_path).unwrap(); + writeln!(wf, "{cb}").unwrap(); + writeln!(wf, "CCCCGGGGTTTTAAAA").unwrap(); + writeln!(wf, "GGGGTTTTAAAACCCC").unwrap(); + } + + let output_dir = tmpdir.path().join("out_crlayout"); + fs::create_dir_all(&output_dir).unwrap(); + let prefix = format!("{}/", output_dir.display()); + + // No --soloOutLayout on the command line: 16 bp CB + 12 bp UMI with a + // whitelist is 10x geometry, so the CellRanger layout is the default. + cargo_bin_cmd!("rustar-aligner") + .args([ + "--runMode", + "alignReads", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--readFilesIn", + cdna_path.to_str().unwrap(), + barcode_path.to_str().unwrap(), + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + wl_path.to_str().unwrap(), + "--soloFeatures", + "Gene", + "--sjdbGTFfile", + gtf.to_str().unwrap(), + "--outFileNamePrefix", + &prefix, + ]) + .assert() + .success(); + + let gunzip = |p: &std::path::Path| -> String { + let f = fs::File::open(p).unwrap_or_else(|e| panic!("{}: {e}", p.display())); + let mut s = String::new(); + flate2::read::GzDecoder::new(f) + .read_to_string(&mut s) + .unwrap(); + s + }; + + // CellRanger's directory names, under outs/, with no per-feature level. + let raw = output_dir.join("outs").join("raw_feature_bc_matrix"); + assert!(raw.is_dir(), "expected {}", raw.display()); + assert!( + !output_dir.join("Solo.out").exists(), + "Solo.out/ should not be written under the CellRanger layout" + ); + + // Every barcode carries the -1 GEM-well suffix, and the raw matrix has one + // column per observed barcode (one), not one per whitelist entry (three). + let barcodes = gunzip(&raw.join("barcodes.tsv.gz")); + assert_eq!(barcodes.lines().count(), 1); + assert_eq!(barcodes.lines().next().unwrap(), format!("{cb}-1")); + + let features = gunzip(&raw.join("features.tsv.gz")); + assert!(features.starts_with("G1\tG1\tGene Expression")); + + // The 2 deduped molecules test_starsolo_gene_matrix asserts, in a matrix + // that is now 1 gene × 1 observed barcode. + let matrix = gunzip(&raw.join("matrix.mtx.gz")); + let dims = matrix.lines().find(|l| !l.starts_with('%')).unwrap(); + assert_eq!(dims, "1 1 1", "unexpected matrix dimensions"); + assert_eq!(matrix.lines().last().unwrap(), "1 1 2"); + + let filt = output_dir.join("outs").join("filtered_feature_bc_matrix"); + let f_barcodes = gunzip(&filt.join("barcodes.tsv.gz")); + assert_eq!(f_barcodes.lines().next().unwrap(), format!("{cb}-1")); + let f_matrix = gunzip(&filt.join("matrix.mtx.gz")); + assert_eq!(f_matrix.lines().last().unwrap(), "1 1 2"); +} + // --------------------------------------------------------------------------- // Test 9a' — Summary.csv stays STARsolo-faithful; the CellRanger mapping funnel // (exonic/intronic/intergenic/antisense) is split out into a separate @@ -1091,6 +1212,11 @@ fn test_starsolo_summary_split() { "Forward", "--sjdbGTFfile", gtf.to_str().unwrap(), + // This fixture's 16 bp CB + 12 bp UMI is 10x geometry, which now + // defaults to CellRanger's output layout; the assertions below are + // about STARsolo's, so state it. + "--soloOutLayout", + "STARsolo", "--outFileNamePrefix", &prefix, ]) @@ -1179,6 +1305,11 @@ fn test_starsolo_sj_feature() { "Forward", "--sjdbGTFfile", gtf.to_str().unwrap(), + // This fixture's 16 bp CB + 12 bp UMI is 10x geometry, which now + // defaults to CellRanger's output layout; the assertions below are + // about STARsolo's, so state it. + "--soloOutLayout", + "STARsolo", "--outFileNamePrefix", &prefix, ]) @@ -1284,6 +1415,11 @@ fn test_starsolo_multimappers() { "Uniform", "--sjdbGTFfile", gtf.to_str().unwrap(), + // This fixture's 16 bp CB + 12 bp UMI is 10x geometry, which now + // defaults to CellRanger's output layout; the assertions below are + // about STARsolo's, so state it. + "--soloOutLayout", + "STARsolo", "--outFileNamePrefix", &prefix, ]) @@ -1523,6 +1659,11 @@ fn test_starsolo_velocyto() { "Forward", "--sjdbGTFfile", gtf.to_str().unwrap(), + // This fixture's 16 bp CB + 12 bp UMI is 10x geometry, which now + // defaults to CellRanger's output layout; the assertions below are + // about STARsolo's, so state it. + "--soloOutLayout", + "STARsolo", "--outFileNamePrefix", &prefix, ]) @@ -1611,6 +1752,11 @@ fn test_starsolo_velocyto_fold_ambiguous() { "Forward", "--sjdbGTFfile", gtf.to_str().unwrap(), + // This fixture's 16 bp CB + 12 bp UMI is 10x geometry, which now + // defaults to CellRanger's output layout; the assertions below are + // about STARsolo's, so state it. + "--soloOutLayout", + "STARsolo", "--outFileNamePrefix", &prefix, ]) @@ -1900,6 +2046,11 @@ fn test_starsolo_cellranger_style_matrix() { "1MM_CR", "--outSAMtype", "SAM", + // This fixture's 16 bp CB + 12 bp UMI is 10x geometry, which now + // defaults to CellRanger's output layout; the assertions below are + // about STARsolo's, so state it. + "--soloOutLayout", + "STARsolo", "--outFileNamePrefix", &prefix, ])