From 51304b55078218f6efe258e6a1f649bb22c4b471 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 10:06:45 +0200 Subject: [PATCH 01/14] 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 02/14] 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 03/14] 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 04/14] 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 05/14] 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 06/14] 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 07/14] 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 08/14] 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, ]) From 0f5b44f72ef39b52b7059c05ee95ed99ac06913a Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 11:54:53 +0200 Subject: [PATCH 09/14] feat(solo): metrics_summary.csv under the CellRanger layout Writes CellRanger 10.0.0's 20 metrics, in its order and value formats (thousands separators, quoted when that adds a comma; percentages to one decimal). STARsolo's Summary.csv is unchanged and still written alongside. Twelve of the 20 match a real cellranger count 10.0.0 run exactly on the 20 000-read fixture. The other eight rest on an interpretation of a denominator CellRanger does not document, and each interpretation is stated in DIVERGENCE.md 3.4 rather than left implicit. Three metrics needed new counters: Q30Stats tallies Phred >= 30 bases over the barcode, the UMI and the cDNA read, populated only when the layout asks for the metrics. The exonic/intronic split needs the gene-body overlap query, so a Gene-only run under this layout now performs it. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 6 + DIVERGENCE.md | 50 ++++++++ src/lib.rs | 26 +++- src/solo/count.rs | 243 +++++++++++++++++++++++++++++++++++- src/solo/mod.rs | 66 +++++++++- tests/alignment_features.rs | 34 +++++ 6 files changed, 415 insertions(+), 10 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a8ba71b..c75867d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,12 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Features +- Under `--soloOutLayout CellRanger`, `metrics_summary.csv` is written with + CellRanger 10.0.0's 20 metrics, in its order and value formats. STARsolo's + `Summary.csv` is unchanged and still written alongside. Twelve of the 20 + match a real `cellranger count` run exactly on the test fixture; the + interpretation behind the other eight is in `DIVERGENCE.md` §3.4. + - `--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 diff --git a/DIVERGENCE.md b/DIVERGENCE.md index 79248cb..74a098d 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -158,6 +158,56 @@ taken from its source. --- +### 3.4 `metrics_summary.csv` under the CellRanger layout (non-STAR) + +**What STAR does.** STARsolo writes `Summary.csv`, its own metric set with its +own names. It has no `metrics_summary.csv`. + +**What rustar-aligner does.** The same, by default. Under +`--soloOutLayout CellRanger` it *additionally* writes +`metrics_summary.csv` with CellRanger 10.0.0's 20 metrics, in CellRanger's +order and value formats. `Summary.csv` is written unchanged alongside it, so +nothing STARsolo-faithful is altered. Under the CellRanger layout the older +`CellRanger.summary.csv` is not written: its four rows are a subset of the +metrics file. + +**Why.** The file is how a 10x pipeline reads run quality. Its 20 names are +also the clearest public statement of what CellRanger measures, which makes it +useful as a target even where our value differs. + +**Impact.** No count changes. It costs one extra gene-body overlap query per +read, because the exonic/intronic split needs it and a `Gene`-only run does not +otherwise do it. + +**Not all 20 are certainties.** Twelve match a real `cellranger count` 10.0.0 +run exactly on the 20 000-read fixture. The other eight follow from our own +tallies under a stated interpretation, because CellRanger does not document the +denominators: + +* `Reads Mapped Confidently to *` reads "confidently" as MAPQ 255, i.e. our + uniquely-mapped set. +* `Reads Mapped Confidently to Transcriptome` is the `Gene` feature's own + uniquely-assigned read tally. +* `Valid UMI Sequences` is measured over reads that reached the UMI check, i.e. + those with a valid barcode. +* `Sequencing Saturation` is `1 - molecules / reads` over the reads that + entered the matrix. This is the largest disagreement on the fixture: 12.8% + against CellRanger's 7.4%. +* `Mean Reads per Cell` is total reads over called cells, not reads-in-cells + over called cells, which is what reproduces CellRanger's value. + +The cell-count-dependent metrics (`Median Genes per Cell`, +`Median UMI Counts per Cell`, `Total Genes Detected`) differ by 2-4 on the +fixture, following the count differences recorded in §1.3 rather than a +different definition. + +**Source.** `src/solo/count.rs` (`write_metrics_summary`, `metric_int`, +`metric_pct`), `src/solo/mod.rs` (`Q30Stats`). CellRanger: the +`metrics_summary.csv` of 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/lib.rs b/src/lib.rs index 6fce173..29090c7 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -2172,8 +2172,13 @@ fn align_reads_solo( stats.record_alignment(0, max_multimaps); stats.record_unmapped_reason(crate::stats::UnmappedReason::Other); // No alignment → barcode still counts toward stats (unmapped → no gene). - let outcome = - solo.process_read(&[], 0, sread.barcode.as_ref(), &[]); + let outcome = solo.process_read( + &[], + 0, + sread.barcode.as_ref(), + &[], + &read.quality, + ); return Ok(SoloReadProduct { sam_records: buffer, per_feature: outcome.per_feature, @@ -2220,6 +2225,7 @@ fn align_reads_solo( transcripts.len(), sread.barcode.as_ref(), &junctions, + &read.quality, ); // Build SAM records for the cDNA alignment (same as SE path). @@ -2559,7 +2565,12 @@ fn align_reads_solo_pe( .iter() .map(|pa| (&pa.mate1_transcript, &pa.mate2_transcript)) .collect(); - solo.process_read_pe(&pairs, pread.barcode.as_ref(), &junctions) + solo.process_read_pe( + &pairs, + pread.barcode.as_ref(), + &junctions, + &pread.mate1.quality, + ) } else if let Some(PairedAlignmentResult::HalfMapped { mapped_transcript, .. @@ -2570,9 +2581,16 @@ fn align_reads_solo_pe( 1, pread.barcode.as_ref(), &junctions, + &pread.mate1.quality, ) } else { - solo.process_read(&[], 0, pread.barcode.as_ref(), &[]) + solo.process_read( + &[], + 0, + pread.barcode.as_ref(), + &[], + &pread.mate1.quality, + ) }; // SAM records (skipped under `--outSAMtype None`). diff --git a/src/solo/count.rs b/src/solo/count.rs index 613af35..784197c 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -1321,8 +1321,12 @@ pub fn write_gene_matrix( .position(|&x| x == f) .map_or(0, |i| ctx.feature_reads[i].load(Ordering::Relaxed)) }; - let have_funnel = ctx.features.contains(&crate::solo::SoloFeature::Gene) - && ctx.features.contains(&crate::solo::SoloFeature::GeneFull); + // The split needs the gene-body query as well as the exon one. Both `Gene` + // + `GeneFull` gives it, and so does `ctx.want_metrics`, which asks + // `process_read` for the body query on its own account. + let have_funnel = ctx.want_metrics + || (ctx.features.contains(&crate::solo::SoloFeature::Gene) + && ctx.features.contains(&crate::solo::SoloFeature::GeneFull)); let region = have_funnel.then(|| RegionFunnel { exonic: ctx.region_stats.exonic.load(Ordering::Relaxed), intronic: ctx.region_stats.intronic.load(Ordering::Relaxed), @@ -1520,10 +1524,32 @@ pub fn write_gene_matrix( reads_of(*feature), )?; log::info!("STARsolo: wrote {}/Summary.csv", feature.dir_name()); + // Under the CellRanger layout the funnel goes into `metrics_summary.csv` + // alongside the other 19 metrics, which is where CellRanger puts it. + if cr_layout && let Some(r) = region { + let invalid_umis = ctx.stats.n_in_umi.load(Ordering::Relaxed) + + ctx.stats.umi_homopolymer.load(Ordering::Relaxed); + write_metrics_summary( + &feature_dir.join("metrics_summary.csv"), + &mstats, + &ctx.q30, + total_reads, + valid_barcodes, + invalid_umis, + mapped_unique, + mapped_multi, + reads_of(crate::solo::SoloFeature::Gene), + r, + )?; + log::info!( + "STARsolo: wrote {}", + feature_dir.join("metrics_summary.csv").display() + ); + } // CellRanger-style mapping funnel goes in a SEPARATE additional file so the // faithful Summary.csv is never altered (PR #90 review: keep this release a // drop-in faithful port; output-changing features come later). - if let Some(r) = region { + if !cr_layout && let Some(r) = region { write_cellranger_summary( &feature_dir.join("CellRanger.summary.csv"), total_reads, @@ -1858,6 +1884,168 @@ struct RegionFunnel { antisense: u64, } +/// Format an integer the way CellRanger's `metrics_summary.csv` does: thousands +/// separated by commas, and quoted when that puts a comma in the field. +/// +/// `200` stays `200`; `20000` becomes `"20,000"`. +fn metric_int(n: u64) -> String { + let digits = n.to_string(); + if digits.len() <= 3 { + return digits; + } + let mut out = String::with_capacity(digits.len() + digits.len() / 3 + 2); + for (i, c) in digits.chars().enumerate() { + if i > 0 && (digits.len() - i).is_multiple_of(3) { + out.push(','); + } + out.push(c); + } + format!("\"{out}\"") +} + +/// Format a fraction as CellRanger does: one decimal place and a percent sign. +fn metric_pct(num: u64, den: u64) -> String { + let f = if den == 0 { + 0.0 + } else { + num as f64 / den as f64 + }; + format!("{:.1}%", f * 100.0) +} + +/// Write CellRanger's `metrics_summary.csv`: one header row of 20 metric names +/// and one row of values, in CellRanger 10.0.0's order. +/// +/// Every metric here is computed from tallies rustar already keeps, or from the +/// Q30 counters added alongside this function. Where CellRanger's definition is +/// not fully documented, the interpretation is stated in a comment on the row +/// rather than left implicit — `DIVERGENCE.md` §3.4 lists the ones that are an +/// interpretation rather than a certainty. +#[allow(clippy::too_many_arguments)] +fn write_metrics_summary( + path: &Path, + mstats: &MatrixStats, + q30: &crate::solo::Q30Stats, + total_reads: u64, + valid_barcodes: u64, + invalid_umis: u64, + mapped_unique: u64, + mapped_multi: u64, + transcriptome_reads: u64, + region: RegionFunnel, +) -> Result<(), Error> { + use std::sync::atomic::Ordering; + + // Cell calling: the same CR2.2 knee `Summary.csv` uses, so the two files + // never disagree about how many cells there are. + let mut umis_desc: Vec = mstats.cells.iter().map(|c| c.n_umis).collect(); + umis_desc.sort_unstable_by(|a, b| b.cmp(a)); + let thr = knee_cr22(&umis_desc, 3000, 0.99, 10.0); + let cells: Vec<&CellStat> = mstats.cells.iter().filter(|c| c.n_umis >= thr).collect(); + let n_cells = cells.len() as u64; + + let reads_counted: u64 = mstats.cells.iter().map(|c| c.n_reads).sum(); + let umis_all: u64 = mstats.cells.iter().map(|c| c.n_umis).sum(); + let reads_in_cells: u64 = cells.iter().map(|c| c.n_reads).sum(); + + let mut umis_sorted: Vec = cells.iter().map(|c| c.n_umis).collect(); + let mut genes_sorted: Vec = cells.iter().map(|c| c.n_genes as u64).collect(); + umis_sorted.sort_unstable(); + genes_sorted.sort_unstable(); + + // Saturation = the share of reads that added no new molecule. + let saturation_num = reads_counted.saturating_sub(umis_all); + let q = |n: &std::sync::atomic::AtomicU64| n.load(Ordering::Relaxed); + + let rows: [(&str, String); 20] = [ + ("Estimated Number of Cells", metric_int(n_cells)), + // CellRanger divides total reads by cells, not reads-in-cells. + ( + "Mean Reads per Cell", + metric_int(total_reads.checked_div(n_cells).unwrap_or(0)), + ), + ( + "Median Genes per Cell", + metric_int(median_sorted(&genes_sorted)), + ), + ("Number of Reads", metric_int(total_reads)), + ("Valid Barcodes", metric_pct(valid_barcodes, total_reads)), + // Of the reads that got as far as the UMI check, i.e. those with a + // valid barcode: the share whose UMI was neither N-containing nor a + // homopolymer. + ( + "Valid UMI Sequences", + metric_pct(valid_barcodes.saturating_sub(invalid_umis), valid_barcodes), + ), + ( + "Sequencing Saturation", + metric_pct(saturation_num, reads_counted), + ), + ( + "Q30 Bases in Barcode", + metric_pct(q(&q30.cb_q30), q(&q30.cb_bases)), + ), + ( + "Q30 Bases in RNA Read", + metric_pct(q(&q30.rna_q30), q(&q30.rna_bases)), + ), + ( + "Q30 Bases in UMI", + metric_pct(q(&q30.umi_q30), q(&q30.umi_bases)), + ), + ( + "Reads Mapped to Genome", + metric_pct(mapped_unique + mapped_multi, total_reads), + ), + // "Confidently" is CellRanger's word for MAPQ 255, which is our + // uniquely-mapped set. + ( + "Reads Mapped Confidently to Genome", + metric_pct(mapped_unique, total_reads), + ), + ( + "Reads Mapped Confidently to Intergenic Regions", + metric_pct(region.intergenic, total_reads), + ), + ( + "Reads Mapped Confidently to Intronic Regions", + metric_pct(region.intronic, total_reads), + ), + ( + "Reads Mapped Confidently to Exonic Regions", + metric_pct(region.exonic, total_reads), + ), + // Uniquely mapped *and* assigned to one gene, which is the `Gene` + // feature's own read tally. + ( + "Reads Mapped Confidently to Transcriptome", + metric_pct(transcriptome_reads, total_reads), + ), + ( + "Reads Mapped Antisense to Gene", + metric_pct(region.antisense, total_reads), + ), + ( + "Fraction Reads in Cells", + metric_pct(reads_in_cells, reads_counted), + ), + ( + "Total Genes Detected", + metric_int(mstats.genes_detected.into()), + ), + ( + "Median UMI Counts per Cell", + metric_int(median_sorted(&umis_sorted)), + ), + ]; + + let header: Vec<&str> = rows.iter().map(|(k, _)| *k).collect(); + let values: Vec<&str> = rows.iter().map(|(_, v)| v.as_str()).collect(); + let out = format!("{}\n{}\n", header.join(","), values.join(",")); + std::fs::write(path, out).map_err(|e| Error::io(e, path))?; + Ok(()) +} + /// Write the STARsolo-faithful `Summary.csv` for one feature: the sequencing / /// genome-mapping rows plus per-cell UMI/gene statistics over the CR2.2-knee-called /// cells. The CellRanger-style exonic/intronic/intergenic/antisense funnel is a @@ -2141,6 +2329,55 @@ mod tests { assert_eq!(std::fs::read_to_string(&sub).unwrap(), "TGCA-1\n"); } + /// The two `metrics_summary.csv` value formats, taken from a real + /// CellRanger 10.0.0 file: integers get thousands separators and are quoted + /// once that puts a comma in the field; fractions get one decimal and a `%`. + #[test] + fn metric_values_are_formatted_the_way_cellranger_formats_them() { + assert_eq!(metric_int(0), "0"); + assert_eq!(metric_int(200), "200"); // "Estimated Number of Cells,200" + assert_eq!(metric_int(999), "999"); + assert_eq!(metric_int(1000), "\"1,000\""); + assert_eq!(metric_int(20_000), "\"20,000\""); // "Number of Reads,\"20,000\"" + assert_eq!(metric_int(1_234_567), "\"1,234,567\""); + + assert_eq!(metric_pct(978, 1000), "97.8%"); // "Valid Barcodes,97.8%" + assert_eq!(metric_pct(1, 1), "100.0%"); + assert_eq!(metric_pct(0, 1000), "0.0%"); + assert_eq!(metric_pct(2, 1000), "0.2%"); + // No reads is 0%, not a division by zero. + assert_eq!(metric_pct(0, 0), "0.0%"); + } + + /// Q30 is Phred ≥ 30 on Phred+33 bytes, i.e. `'?'` and above, counted + /// separately for the barcode, the UMI and the cDNA read. + #[test] + fn q30_counts_bases_at_phred_30_and_above() { + use crate::solo::{CellBarcode, Q30Stats}; + use std::sync::atomic::Ordering; + + let bc = CellBarcode { + cb_seq: vec![0, 1, 2, 3], + // '>' is Phred 29, '?' is Phred 30, 'I' is Phred 40. + cb_qual: b">?II".to_vec(), + umi_seq: vec![0, 1], + umi_qual: b">>".to_vec(), + }; + let q = Q30Stats::default(); + q.record(Some(&bc), b"II>I"); + assert_eq!(q.cb_bases.load(Ordering::Relaxed), 4); + assert_eq!(q.cb_q30.load(Ordering::Relaxed), 3); + assert_eq!(q.umi_bases.load(Ordering::Relaxed), 2); + assert_eq!(q.umi_q30.load(Ordering::Relaxed), 0); + assert_eq!(q.rna_bases.load(Ordering::Relaxed), 4); + assert_eq!(q.rna_q30.load(Ordering::Relaxed), 3); + + // A read with no barcode still contributes its cDNA bases. + q.record(None, b"I"); + assert_eq!(q.cb_bases.load(Ordering::Relaxed), 4); + assert_eq!(q.rna_bases.load(Ordering::Relaxed), 5); + } + #[test] fn median_sorted_odd_even_empty() { assert_eq!(median_sorted(&[]), 0); diff --git a/src/solo/mod.rs b/src/solo/mod.rs index 2253dee..93ac590 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -534,6 +534,9 @@ impl SoloRecorder { /// Everything the alignment loop needs to quantify a solo run, shared as an /// `Arc` across rayon threads. The gene model is built from `--sjdbGTFfile`; /// the whitelist and stats are read concurrently (interior atomics). +// The bools are independent feature switches read on the hot path; grouping +// them into a struct would add an indirection for no clarity. +#[allow(clippy::struct_excessive_bools)] pub struct SoloContext { pub layout: SoloBarcodeLayout, pub whitelist: CbWhitelist, @@ -564,6 +567,48 @@ pub struct SoloContext { /// `--soloMultiMappers` includes a non-`Unique` method → capture gene- /// ambiguous reads for distribution into `UniqueAndMult-*.mtx`. pub want_multi: bool, + /// `--soloOutLayout CellRanger`: write `metrics_summary.csv`, which needs + /// the Q30 tallies and the full positional funnel. + pub want_metrics: bool, + /// Base-quality tallies for the three `Q30 Bases in ...` metrics. + pub q30: Q30Stats, +} + +/// Bases seen and bases at Phred ≥ 30, split the way CellRanger's +/// `metrics_summary.csv` splits them: the cell barcode, the UMI, and the cDNA +/// read. Only populated when `--soloOutLayout CellRanger` asks for the metrics. +#[derive(Default)] +pub struct Q30Stats { + pub cb_bases: AtomicU64, + pub cb_q30: AtomicU64, + pub umi_bases: AtomicU64, + pub umi_q30: AtomicU64, + pub rna_bases: AtomicU64, + pub rna_q30: AtomicU64, +} + +impl Q30Stats { + /// Tally one read's barcode-read and cDNA-read qualities. + /// + /// Qualities are raw FASTQ bytes, Phred+33, so a Phred of 30 is the byte + /// `'?'` (63). The cDNA quality is the **unclipped** read, which is what + /// CellRanger reports on. + pub fn record(&self, bc: Option<&CellBarcode>, rna_qual: &[u8]) { + const Q30: u8 = 30 + 33; + let tally = |bases: &AtomicU64, q30: &AtomicU64, q: &[u8]| { + if q.is_empty() { + return; + } + bases.fetch_add(q.len() as u64, Ordering::Relaxed); + let n = q.iter().filter(|&&b| b >= Q30).count() as u64; + q30.fetch_add(n, Ordering::Relaxed); + }; + if let Some(bc) = bc { + tally(&self.cb_bases, &self.cb_q30, &bc.cb_qual); + tally(&self.umi_bases, &self.umi_q30, &bc.umi_qual); + } + tally(&self.rna_bases, &self.rna_q30, rna_qual); + } } /// Per-region read tallies for the `Summary.csv` mapping funnel (uniquely-mapped @@ -678,6 +723,7 @@ impl SoloContext { let sj_enabled = params.solo_features.iter().any(|f| f == "SJ"); let velocyto_enabled = params.solo_features.iter().any(|f| f == "Velocyto"); let want_multi = params.solo_multi_mappers.iter().any(|m| m != "Unique"); + let want_metrics = params.solo_out_layout == "CellRanger"; Ok(Self { layout: SoloBarcodeLayout::from_params(params), @@ -695,6 +741,8 @@ impl SoloContext { velocyto_enabled, velocyto_records: Mutex::new(Vec::new()), want_multi, + want_metrics, + q30: Q30Stats::default(), }) } @@ -720,15 +768,26 @@ impl SoloContext { n_loci: usize, barcode: Option<&CellBarcode>, junctions: &[(u64, u64)], + rna_qual: &[u8], ) -> SoloReadOutcome { let mut out = SoloReadOutcome::default(); + // Base qualities for `metrics_summary.csv`. Tallied before any early + // return, so every read reaching this point is counted exactly once. + if self.want_metrics { + self.q30.record(barcode, rna_qual); + } + // One-pass classification: the two overlap queries are shared between the // per-feature gene assignment and the CellRanger-style mapping funnel, so // this is no more work than the old per-feature `assign_gene_se` calls. - let want_exon = self.features.contains(&SoloFeature::Gene); + let want_exon = self.features.contains(&SoloFeature::Gene) || self.want_metrics; // Velocyto assigns its gene by gene-body overlap, so it needs `want_body`. - let want_body = self.features.contains(&SoloFeature::GeneFull) || self.velocyto_enabled; + // `metrics_summary.csv` reports the exonic/intronic split, which is the + // same query, so it asks for it too. + let want_body = self.features.contains(&SoloFeature::GeneFull) + || self.velocyto_enabled + || self.want_metrics; let class = classify_read( cdna_transcripts, &self.gene_ann, @@ -887,6 +946,7 @@ impl SoloContext { pairs: &[(&Transcript, &Transcript)], barcode: Option<&CellBarcode>, junctions: &[(u64, u64)], + rna_qual: &[u8], ) -> SoloReadOutcome { // Effective transcripts: give both mates of a pair the pair's strand // (mate 1's), so `classify_read`'s per-transcript strand filter treats them @@ -898,7 +958,7 @@ impl SoloContext { eff.push((*m1).clone()); eff.push(m2c); } - self.process_read(&eff, pairs.len(), barcode, junctions) + self.process_read(&eff, pairs.len(), barcode, junctions, rna_qual) } } diff --git a/tests/alignment_features.rs b/tests/alignment_features.rs index a0f812d..7355401 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -1147,6 +1147,40 @@ fn test_solo_out_layout_cellranger() { 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"); + + // metrics_summary.csv: CellRanger 10.0.0's 20 metrics, in its order, as a + // header row and a value row. The header is compared against the literal + // string from a real CellRanger run. + let metrics = fs::read_to_string(output_dir.join("outs").join("metrics_summary.csv")).unwrap(); + let mut lines = metrics.lines(); + assert_eq!( + lines.next().unwrap(), + "Estimated Number of Cells,Mean Reads per Cell,Median Genes per Cell,\ + Number of Reads,Valid Barcodes,Valid UMI Sequences,Sequencing Saturation,\ + Q30 Bases in Barcode,Q30 Bases in RNA Read,Q30 Bases in UMI,\ + Reads Mapped to Genome,Reads Mapped Confidently to Genome,\ + Reads Mapped Confidently to Intergenic Regions,\ + Reads Mapped Confidently to Intronic Regions,\ + Reads Mapped Confidently to Exonic Regions,\ + Reads Mapped Confidently to Transcriptome,Reads Mapped Antisense to Gene,\ + Fraction Reads in Cells,Total Genes Detected,Median UMI Counts per Cell" + ); + let values: Vec<&str> = lines.next().unwrap().split(',').collect(); + // 20 metrics; 8 reads, so no field is large enough to be comma-quoted. + assert_eq!(values.len(), 20); + assert_eq!(values[0], "1", "one called cell"); + assert_eq!(values[3], "8", "8 reads"); + assert_eq!(values[4], "100.0%", "all barcodes valid"); + assert_eq!(values[5], "100.0%", "all UMIs valid"); + // 8 reads collapsing to 2 molecules: 6 of the 8 added nothing. + assert_eq!(values[6], "75.0%", "sequencing saturation"); + // The fixture writes 'I' (Phred 40) for every base. + assert_eq!(values[7], "100.0%", "Q30 in barcode"); + assert_eq!(values[8], "100.0%", "Q30 in RNA read"); + assert_eq!(values[9], "100.0%", "Q30 in UMI"); + assert_eq!(values[18], "1", "one gene detected"); + assert_eq!(values[19], "2", "median 2 UMIs per cell"); + assert!(lines.next().is_none(), "exactly two rows"); } // --------------------------------------------------------------------------- From 7ac900093a8e4af958e9f48bf48e4db2edaaf10f Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 15:04:52 +0200 Subject: [PATCH 10/14] docs(divergence): pin down what differs in Sequencing Saturation 10x define it as 1 - n_deduped_reads / n_reads over unique (barcode, UMI, gene) combinations among confidently mapped reads. Implementing that literally, counting distinct triples before UMI correction, gives 0.0% on the fixture: no two reads there share an exact triple. So the literal reading is wrong and their numerator is the corrected molecule count, like ours. What differs is the read denominator, which 3.4 now says instead of leaving the disagreement unexplained. Co-Authored-By: Claude Opus 5 (1M context) --- DIVERGENCE.md | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/DIVERGENCE.md b/DIVERGENCE.md index 74a098d..3dc5d62 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -191,8 +191,16 @@ denominators: * `Valid UMI Sequences` is measured over reads that reached the UMI check, i.e. those with a valid barcode. * `Sequencing Saturation` is `1 - molecules / reads` over the reads that - entered the matrix. This is the largest disagreement on the fixture: 12.8% - against CellRanger's 7.4%. + entered the matrix. This is the largest disagreement on the fixture: 12.7% + against CellRanger's 7.4%. 10x define it as + `1 - n_deduped_reads / n_reads`, with `n_deduped_reads` the number of unique + `(barcode, UMI, gene)` combinations among confidently mapped reads. Taking + that literally — counting distinct triples *before* UMI correction — gives + **0.0%** on this fixture, because no two reads here share an exact triple, so + the literal reading is ruled out and their numerator is the corrected + molecule count, as ours is. The residual is therefore the denominator: at + 7.4% theirs implies about 16 300 reads where ours counts about 17 300. The + difference is which reads are "confidently mapped", not the formula. * `Mean Reads per Cell` is total reads over called cells, not reads-in-cells over called cells, which is what reproduces CellRanger's value. From 9a3092e48f55e587e26a73bf0f88e8cf73da4398 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 14:15:08 +0200 Subject: [PATCH 11/14] feat(solo): --soloCellFilter OrdMag, CellRanger's cell call The same quantile-over-ratio rule STARsolo calls CellRanger2.2, but with the expected cell count searched for instead of fixed at 3 000: minimise (OrdMag(x) - x)^2 / x over x from 2 to maxExpectedCells, where OrdMag(x) is the number of cells the rule calls when told to expect x of them. The loss is small where the rule predicts itself. EmptyDrops_CR now uses it for its initial cell set, which is the order CellRanger runs the two steps in. CellRanger2.2 is untouched and stays the default. Two things 10x's page leaves unspecified are decided here and documented: every integer in range is evaluated rather than a geometric grid, so the result is the true minimum of the stated loss and not a nearby grid point; and ties go to the smaller x, so nothing depends on iteration order. On the 20 000-read fixture this changes nothing: CellRanger2.2, OrdMag and EmptyDrops_CR all call 200 cells, as does CellRanger 10.0.0. That fixture has a clean plateau where any of these rules works. The two part company on a graded distribution, which is what the unit tests cover. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 7 +++ DIVERGENCE.md | 38 ++++++++++++ src/params/mod.rs | 10 ++- src/solo/count.rs | 151 +++++++++++++++++++++++++++++++++++++++++++++- 4 files changed, 203 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c75867d..8f4c96d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,13 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Features +- `--soloCellFilter OrdMag` implements CellRanger's cell call: the same + quantile-over-ratio rule as `CellRanger2.2`, but searching for the expected + cell count by minimising `(OrdMag(x) - x)^2 / x` instead of fixing it at + 3 000. `EmptyDrops_CR` now uses it for its initial cell set, which is the + order CellRanger runs the two steps in; `CellRanger2.2` is unchanged and + remains the default. See `DIVERGENCE.md` §3.5. + - Under `--soloOutLayout CellRanger`, `metrics_summary.csv` is written with CellRanger 10.0.0's 20 metrics, in its order and value formats. STARsolo's `Summary.csv` is unchanged and still written alongside. Twelve of the 20 diff --git a/DIVERGENCE.md b/DIVERGENCE.md index 3dc5d62..7d6f44d 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -216,6 +216,44 @@ rather than taken from its source. --- +### 3.5 `--soloCellFilter OrdMag`, and EmptyDrops_CR's initial cell set (non-STAR) + +**What STAR does.** STARsolo's cell calling is `CellRanger2.2`: take the +`quantile` of the top `nExpectedCells` barcodes by UMI count and call every +barcode holding at least that over `maxMinRatio`, with `nExpectedCells` fixed +at 3 000. `EmptyDrops_CR` uses the same knee for the set of cells it guarantees +before the Monte-Carlo rescue. + +**What CellRanger does.** The same rule, but it *searches* for the expected +cell count instead of assuming one: it minimises `(OrdMag(x) - x)^2 / x` over +`x` from 2 to about 45 000, where `OrdMag(x)` is the number of cells the rule +calls when told to expect `x`. The loss is small where the rule predicts +itself. EmptyDrops then rescues barcodes below that cutoff. + +**What rustar-aligner does.** `--soloCellFilter OrdMag maxExpectedCells +quantile ratio` (default `45000 0.99 10`) implements the search, and +`EmptyDrops_CR` now uses it for its initial set, which is the order CellRanger +runs the two steps in. `CellRanger2.2` is untouched and remains the default. + +**Two things the 10x page does not specify, decided here.** Every integer in +range is evaluated rather than a geometric grid, so the result is the true +minimum of the stated loss rather than a nearby grid point; each evaluation is +a binary search over the sorted totals, so an exhaustive sweep costs nothing +measurable. And ties go to the smaller `x`, so the result does not depend on +iteration order — determinism, as everywhere else in this codebase. + +**Impact.** `EmptyDrops_CR`'s initial cell set can differ from what it was. +On the 20 000-read fixture it does not: `CellRanger2.2`, `OrdMag` and +`EmptyDrops_CR` all call 200 cells, as does CellRanger 10.0.0. That fixture has +a clean plateau, where any of these rules works; the two rules part company on +a graded distribution with no plateau, which is what the unit tests cover. + +**Source.** `src/solo/count.rs` (`ordmag_at`, `ordmag_threshold`). CellRanger: +the Gene Expression algorithm page, "Cell Calling", read rather than taken from +its source. Coverage of the rest of that page is tracked in #181. + +--- + ## 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 edfa3c9..b85e6ed 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1113,7 +1113,15 @@ pub struct Parameters { #[arg(long = "soloCBmatchWLtype", default_value = "1MM_multi")] pub solo_cb_match_wl_type: String, - /// Cell-calling / matrix filtering: None, CellRanger2.2, EmptyDrops_CR, TopCells. + /// Cell-calling / matrix filtering: None, CellRanger2.2, EmptyDrops_CR, + /// TopCells, OrdMag. + /// + /// `OrdMag` is CellRanger's own initial cell call and is **not a STAR + /// method**: the same quantile-over-ratio rule as `CellRanger2.2`, but with + /// the expected cell count searched for rather than fixed at 3 000. Its + /// arguments are `maxExpectedCells quantile ratio`, default + /// `45000 0.99 10`. `EmptyDrops_CR` now uses it for its initial set, which + /// is the order CellRanger runs the two steps in. #[arg(long = "soloCellFilter", num_args = 1.., default_values_t = vec!["CellRanger2.2".to_string(), "3000".to_string(), "0.99".to_string(), "10".to_string()])] pub solo_cell_filter: Vec, diff --git a/src/solo/count.rs b/src/solo/count.rs index 784197c..a4b9314 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -998,6 +998,61 @@ fn knee_cr22(umis_desc: &[u64], n_expected: usize, max_pct: f64, max_min_ratio: (robust_max / max_min_ratio).ceil() as u64 } +/// How many cells the order-of-magnitude rule calls if `n` cells are expected, +/// and the UMI cutoff it used. +/// +/// The rule, from 10x's Gene Expression algorithm page: take `m`, the +/// `quantile` of the top `n` barcodes by UMI count, and call every barcode +/// holding at least `m / ratio`. With the defaults that is "the 99th percentile +/// of the top n, divided by ten", the same shape as `knee_cr22` but with `n` a +/// variable rather than a fixed 3 000. +fn ordmag_at(umis_desc: &[u64], n: usize, quantile: f64, ratio: f64) -> (usize, u64) { + if umis_desc.is_empty() || n == 0 { + return (0, 0); + } + let idx = ((n as f64 * (1.0 - quantile)).round() as usize).min(umis_desc.len() - 1); + let cutoff = ((umis_desc[idx] as f64 / ratio).round() as u64).max(1); + // `umis_desc` is descending, so the called set is its prefix. + let called = umis_desc.partition_point(|&u| u >= cutoff); + (called, cutoff) +} + +/// CellRanger's OrdMag cell-calling threshold: the UMI cutoff at the expected +/// cell count that best predicts itself. +/// +/// 10x describes this as minimising `(OrdMag(x) - x)^2 / x` over a search from +/// 2 to about 45 000 cells. `OrdMag(x)` is the number of cells the rule above +/// calls when told to expect `x` of them, so the loss is small where the rule +/// is self-consistent: feed it the right number of cells and it gives that +/// number back. +/// +/// Two choices this makes explicit, because the page does not state them: +/// +/// * **Every integer in range is evaluated**, rather than a geometric grid +/// refined by search. Each evaluation is a binary search over the sorted +/// totals, so an exhaustive sweep of 45 000 candidates costs nothing +/// measurable and finds the true minimum of the stated loss rather than a +/// grid point near it. +/// * **Ties go to the smaller `x`**, so the result does not depend on iteration +/// order. +fn ordmag_threshold(umis_desc: &[u64], max_expected: usize, quantile: f64, ratio: f64) -> u64 { + if umis_desc.is_empty() { + return 0; + } + let hi = max_expected.min(umis_desc.len()).max(2); + let mut best: Option<(f64, u64)> = None; + for x in 2..=hi { + let (called, cutoff) = ordmag_at(umis_desc, x, quantile, ratio); + let d = called as f64 - x as f64; + let loss = d * d / x as f64; + // Strictly-less keeps the first (smallest) x on a tie. + if best.is_none_or(|(b, _)| loss < b) { + best = Some((loss, cutoff)); + } + } + best.map_or(0, |(_, cutoff)| cutoff) +} + /// Whitelist indices of called cells (sorted ascending) per `--soloCellFilter`. /// `None` → no filtered/ output. `EmptyDrops_CR` writes only the knee-guaranteed /// cells here (the Monte-Carlo rescue is the standalone `emptydrops` binary). @@ -1012,6 +1067,17 @@ fn called_cells(cells: &[CellStat], filter: &[String]) -> Option> { idx.sort_by(|a, b| b.n_umis.cmp(&a.n_umis).then(a.cb.cmp(&b.cb))); idx.into_iter().take(n).map(|c| c.cb).collect() } + // CellRanger's own initial cell set, searched rather than assumed. + "OrdMag" => { + let mut umis: Vec = cells.iter().map(|c| c.n_umis).collect(); + umis.sort_unstable_by(|a, b| b.cmp(a)); + let thr = ordmag_threshold(&umis, arg(1, 45000.0) as usize, arg(2, 0.99), arg(3, 10.0)); + cells + .iter() + .filter(|c| c.n_umis >= thr) + .map(|c| c.cb) + .collect() + } // EmptyDrops_CR is handled by `emptydrops_called`; the knee here is the // fallback / guaranteed-cell base. "CellRanger2.2" | "EmptyDrops_CR" => { @@ -1051,7 +1117,9 @@ fn emptydrops_called( .and_then(|s| s.parse::().ok()) .unwrap_or(d) }; - let (n_expected, max_pct, ratio) = (arg(1, 3000.0) as usize, arg(2, 0.99), arg(3, 10.0)); + // `nExpectedCells` (arg 1) is the 2.2 knee's fixed cell count; OrdMag + // searches for it instead, so this path no longer reads it. + let (max_pct, ratio) = (arg(2, 0.99), arg(3, 10.0)); let (ind_min, ind_max) = (arg(4, 45000.0) as usize, arg(5, 90000.0) as usize); let umi_min = arg(6, 500.0) as u64; let umi_min_frac = arg(7, 0.01); @@ -1063,7 +1131,12 @@ fn emptydrops_called( let mut order: Vec<&CellStat> = cells.iter().collect(); order.sort_by(|a, b| b.n_umis.cmp(&a.n_umis).then(a.cb.cmp(&b.cb))); let totals_desc: Vec = order.iter().map(|c| c.n_umis).collect(); - let thr = knee_cr22(&totals_desc, n_expected, max_pct, ratio); + // CellRanger's initial cell set is OrdMag, not the 2.2 knee: the same + // quantile-over-ratio rule, but with the expected cell count searched for + // rather than fixed. EmptyDrops then rescues barcodes below it. `ind_min` + // is already the ~45 000 upper bound 10x states for that search, so the + // parameter list is unchanged. + let thr = ordmag_threshold(&totals_desc, ind_min, max_pct, ratio); let n_simple = totals_desc.iter().take_while(|&&u| u >= thr).count(); let mut called: Vec = order.iter().take(n_simple).map(|c| c.cb).collect(); @@ -2442,6 +2515,80 @@ mod tests { ); } + /// OrdMag finds the expected cell count that predicts itself. On a clean + /// population of 100 cells over an ambient tail, the loss is zero at 100 + /// and the cutoff is a tenth of the plateau. + #[test] + fn ordmag_finds_the_self_consistent_cell_count() { + let mut umis: Vec = vec![1000; 100]; + umis.extend(std::iter::repeat_n(10u64, 5000)); + umis.sort_unstable_by(|a, b| b.cmp(a)); + + // At x = 100 the 99th percentile is umis[1] = 1000, cutoff 100, and + // exactly 100 barcodes clear it: the loss is 0. + let (called, cutoff) = ordmag_at(&umis, 100, 0.99, 10.0); + assert_eq!((called, cutoff), (100, 100)); + + assert_eq!(ordmag_threshold(&umis, 45000, 0.99, 10.0), 100); + assert_eq!(umis.iter().filter(|&&u| u >= 100).count(), 100); + } + + /// The 2.2 knee is OrdMag with the search removed, so they agree when the + /// fixed guess happens to be right and part company when it is not. Here + /// 20 000 cells sit far from the knee's hardcoded 3 000. + #[test] + fn ordmag_beats_a_fixed_expected_cell_count_when_the_guess_is_wrong() { + let mut umis: Vec = vec![5000; 20_000]; + umis.extend(std::iter::repeat_n(5u64, 100_000)); + umis.sort_unstable_by(|a, b| b.cmp(a)); + + // The knee looks at the top 3 000 only, so its "99th percentile" is + // umis[30], deep inside the plateau: same cutoff here, by luck of a + // flat population. + assert_eq!(knee_cr22(&umis, 3000, 0.99, 10.0), 500); + assert_eq!(ordmag_threshold(&umis, 45000, 0.99, 10.0), 500); + // Both call all 20 000, which is the point: on an easy distribution the + // search changes nothing. It earns its keep on the graded one below. + assert_eq!(umis.iter().filter(|&&u| u >= 500).count(), 20_000); + } + + /// A graded distribution with no plateau, where the fixed guess and the + /// search disagree. This is the case the search exists for. + #[test] + fn ordmag_and_the_fixed_knee_disagree_on_a_graded_distribution() { + // 10 000 barcodes falling geometrically, then ambient. + let mut umis: Vec = (0..10_000) + .map(|i| (100_000.0 * 0.9995_f64.powi(i)) as u64) + .collect(); + umis.extend(std::iter::repeat_n(2u64, 50_000)); + umis.sort_unstable_by(|a, b| b.cmp(a)); + + let knee = knee_cr22(&umis, 3000, 0.99, 10.0); + let ord = ordmag_threshold(&umis, 45000, 0.99, 10.0); + assert_ne!(knee, ord, "the search should not reproduce the fixed guess"); + // The self-consistency check: the count OrdMag calls is what OrdMag was + // solving for, within the granularity of the distribution. + let called = umis.iter().filter(|&&u| u >= ord).count(); + let (recall, _) = ordmag_at(&umis, called, 0.99, 10.0); + assert!( + (recall as i64 - called as i64).abs() * 20 <= called as i64, + "OrdMag called {called} cells but predicts {recall} from that count" + ); + } + + /// Degenerate inputs return a threshold rather than panicking on an empty + /// slice or a zero expected count. + #[test] + fn ordmag_handles_empty_and_tiny_inputs() { + assert_eq!(ordmag_threshold(&[], 45000, 0.99, 10.0), 0); + assert_eq!(ordmag_at(&[], 10, 0.99, 10.0), (0, 0)); + assert_eq!(ordmag_at(&[100], 0, 0.99, 10.0), (0, 0)); + // One barcode: the cutoff is a tenth of it, and it clears its own bar. + assert_eq!(ordmag_at(&[100], 1, 0.99, 10.0), (1, 10)); + // The cutoff never drops below 1, so a zero-UMI barcode is never a cell. + assert_eq!(ordmag_at(&[1, 0], 2, 0.99, 10.0), (1, 1)); + } + #[test] fn knee_cr22_threshold() { // 100 cells at 1000 UMI, then a long ambient tail at 10. From 0e71544603b174ac46987e1ace77ed521dd78217 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 20:31:16 +0200 Subject: [PATCH 12/14] feat(solo): count intronic reads by default on 10x geometry CellRanger has included intronic reads in its gene counts since v7.0. STARsolo's Gene feature is exonic-only, so a 10x run compared against CellRanger came out 30.5% low. --soloFeatures now defaults to GeneFull on 10x geometry, alongside the other CellRanger defaults; naming the flag explicitly still wins. Measured on 10x's own pbmc_1k_v3, 20 M read pairs, refdata-gex-GRCh38-2024-A on both sides, against cellranger count 10.0.0: Gene 6 382 961 counts, -30.5%, 52.9% of entries identical GeneFull 9 039 161 counts, -1.6%, 91.8% of entries identical The yeast fixture could not show this: yeast has almost no introns, and there GeneFull counted fewer than Gene because its dense overlapping genes turn unambiguous exonic reads into ambiguous gene-body ones. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 7 +++++++ src/params/mod.rs | 40 +++++++++++++++++++++++++++++++++++++++- 2 files changed, 46 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8f4c96d..7b46fd3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,13 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Features +- On 10x geometry, `--soloFeatures` now defaults to `GeneFull` rather than + `Gene`: CellRanger has counted intronic reads toward the gene since v7.0. + Measured on 10x's `pbmc_1k_v3`, exonic-only counting is 30.5% below + CellRanger and `GeneFull` brings it to 1.6%. **This changes default counts + on 10x runs**; naming `--soloFeatures` explicitly still wins. See + `DIVERGENCE.md` §1.3. + - `--soloCellFilter OrdMag` implements CellRanger's cell call: the same quantile-over-ratio rule as `CellRanger2.2`, but searching for the expected cell count by minimising `(OrdMag(x) - x)^2 / x` instead of fixing it at diff --git a/src/params/mod.rs b/src/params/mod.rs index b85e6ed..b87b2c4 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1921,7 +1921,7 @@ 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); 6] = [ +const CELLRANGER_DEFAULTS: [(&str, &str); 7] = [ ("clip_adapter_type", "CellRanger4"), ("out_filter_score_min", "30"), ("solo_cb_match_wl_type", "1MM_multi_Nbase_pseudocounts"), @@ -1930,6 +1930,10 @@ const CELLRANGER_DEFAULTS: [(&str, &str); 6] = [ // 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"), + // CellRanger has counted intronic reads toward the gene since v7.0. + // STARsolo's `Gene` is exonic-only, which on human data is 30% below + // CellRanger; `GeneFull` is the equivalent of its default. + ("solo_features", "GeneFull"), ]; /// Does this look like a 10x Chromium run? @@ -1982,6 +1986,7 @@ fn apply_cellranger_defaults_on_10x(params: &mut Parameters, matches: &clap::Arg "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(), + "solo_features" => params.solo_features = vec![value.to_string()], _ => continue, } applied.push(value); @@ -2062,6 +2067,39 @@ mod tests { 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"); + assert_eq!(p.solo_features, vec!["GeneFull".to_string()]); + } + + /// `--soloFeatures` given explicitly wins, so a user who wants STARsolo's + /// exonic-only counting on 10x data still gets it. + #[test] + fn an_explicit_solo_features_beats_the_10x_intron_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", + "--soloFeatures", + "Gene", + ]) + .unwrap(); + assert_eq!(p.solo_features, vec!["Gene".to_string()]); + // The rest of the CellRanger defaults still apply. + assert_eq!(p.solo_umi_dedup, vec!["1MM_CR".to_string()]); } /// `--soloOutLayout CellRanger` pulls three existing flags and the output From 8dec2fdabbc464fa39ae65e93338e2c6970f1b27 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 21:29:40 +0200 Subject: [PATCH 13/14] perf(align): reuse cluster_seeds' window-bin map across reads cluster_seeds built a fresh FxHashMap per read. The pre-sizing to anchor_indices.len() * 2 was only a floor: merging two windows re-keys every bin in the merged span, so a read with wide windows inserts far more entries than it has anchors and the map rehashes. Sampling a human 10x run put hashbrown::reserve_rehash at 2.1% of on-CPU time, all of it here. The map is now kept per thread and cleared per read, so its capacity settles at what the workload needs and the growth is paid once per thread. It is taken out of the thread-local rather than borrowed across the body, so a nested call would get a fresh map instead of panicking on an active borrow. Output-neutral: matrix.mtx, barcodes.tsv, features.tsv and SJ.out.tab are byte-identical before and after on 20 M read pairs of 10x pbmc_1k_v3 against GRCh38, compared decompressed. Interleaved A/B, 8 pairs: new faster in 8 of 8, median -7.0%. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 8 ++++++++ src/align/stitch.rs | 27 +++++++++++++++++++++++---- 2 files changed, 31 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7b46fd3..02a15ae 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,14 @@ Sections commonly used: Features, Bug fixes, Other changes. ## [Unreleased] +### Other changes + +- `cluster_seeds` reuses its window-bin map across reads on a thread instead + of rebuilding it per read. Merging two windows re-keys every bin in the + merged span, so the per-read pre-sizing was only a floor and the map + rehashed; profiling a human 10x run put that rehash at 2.1% of on-CPU time. + Output-neutral, verified by an empty diff on 20 M read pairs. + ### Features - On 10x geometry, `--soloFeatures` now defaults to `GeneFull` rather than diff --git a/src/align/stitch.rs b/src/align/stitch.rs index 1c25ecc..f2361bf 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -461,7 +461,7 @@ pub fn cluster_seeds( ) -> Vec { // Integer-keyed maps in this hot path (the #1 align hotspot, ~19% self-time): // FxHash, not the default SipHash, and pre-sized to avoid rehashing. - use rustc_hash::{FxBuildHasher, FxHashMap}; + use rustc_hash::FxHashMap; let win_bin_nbits = params.win_bin_nbits; let win_anchor_dist_nbins = params.win_anchor_dist_nbins; @@ -574,9 +574,25 @@ pub fn cluster_seeds( // At most one window per anchor position — pre-size to avoid growth reallocs. let mut windows: Vec = Vec::with_capacity(anchor_indices.len()); // winBin: (strand, bin) → window_index - // Chromosome is implicit since bins are from absolute forward positions - let mut win_bin: FxHashMap<(bool, u64), usize> = - FxHashMap::with_capacity_and_hasher(anchor_indices.len() * 2, FxBuildHasher); + // Chromosome is implicit since bins are from absolute forward positions. + // + // Reused across reads on this thread. The pre-sizing below is only a floor: + // merging two windows re-keys every bin in the merged span, so a read with + // wide windows inserts far more entries than it has anchors and the map + // rehashes. Profiling a human 10x run put `hashbrown::reserve_rehash` at + // 2.1% of on-CPU time, all of it here. Keeping the allocation between reads + // lets the capacity settle at whatever the workload needs and pays that + // cost once per thread instead of once per read. + // + // Taken out of the cell rather than borrowed across the body, so a nested + // call (there is none today) would get a fresh map instead of panicking. + thread_local! { + static WIN_BIN: std::cell::RefCell> = + std::cell::RefCell::new(FxHashMap::default()); + } + let mut win_bin = WIN_BIN.with(|c| std::mem::take(&mut *c.borrow_mut())); + win_bin.clear(); + win_bin.reserve(anchor_indices.len() * 2); for &anchor_idx in &anchor_indices { let anchor = &seeds[anchor_idx]; @@ -711,6 +727,7 @@ pub fn cluster_seeds( } if windows.iter().all(|w| !w.alive) { + WIN_BIN.with(|c| *c.borrow_mut() = win_bin); return Vec::new(); } @@ -1064,6 +1081,8 @@ pub fn cluster_seeds( }); } + // Hand the map back with its capacity, for the next read on this thread. + WIN_BIN.with(|c| *c.borrow_mut() = win_bin); clusters } From babe77c1d465d0f900d4d281b42260a4cac1c06a Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 31 Jul 2026 22:31:19 +0200 Subject: [PATCH 14/14] perf(align): move cluster_seeds' alignments into their cluster The per-window Vec was cloned into the SeedCluster while `windows` was dropped one statement later, so every read paid a full copy per window to throw the original away. std::mem::take moves it instead. Output-neutral, verified by an empty diff on 20 M read pairs of 10x pbmc_1k_v3 against GRCh38. Co-Authored-By: Claude Opus 5 (1M context) --- src/align/stitch.rs | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/align/stitch.rs b/src/align/stitch.rs index f2361bf..63aa90c 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -1065,13 +1065,17 @@ pub fn cluster_seeds( // Phase 5: Build SeedCluster output let mut clusters = Vec::with_capacity(windows.len()); - for window in &windows { + // `windows` is dropped at the end of this function, so each window's + // alignments move into its cluster rather than being cloned. The clone was + // a full Vec copy per window per read, thrown away one + // statement later. + for window in &mut windows { if !window.alive || window.alignments.is_empty() { continue; } clusters.push(SeedCluster { - alignments: window.alignments.clone(), + alignments: std::mem::take(&mut window.alignments), chr_idx: window.chr_idx, genome_start: window.actual_start, genome_end: window.actual_end,