perf(align): reuse cluster_seeds' window-bin map across reads (-7% median, output-neutral)

CHANGELOG.md

@@ -11,8 +11,57 @@ 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
+  `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
+  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
+  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
+  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
+  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
+  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).
@@ -99,6 +148,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

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]

DIVERGENCE.md

@@ -65,6 +65,195 @@ 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.96% above CellRanger to
+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
+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
+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.
+
+---
+
+### 3.3 `--soloOutLayout CellRanger` (non-STAR; on by default on 10x geometry)
+
+**What STAR does.** STARsolo writes
+`Solo.out/<feature>/{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.
+
+---
+
+### 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.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.
+
+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.
+
+---
+
+### 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.

src/align/stitch.rs

@@ -461,7 +461,7 @@ pub fn cluster_seeds(
 ) -> Vec<SeedCluster> {
     // 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<Window> = 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<FxHashMap<(bool, u64), usize>> =
+            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();
     }
 
@@ -1048,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<WindowAlignment> 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,
@@ -1064,6 +1085,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
 }