STARsolo Velocyto and stats files as STARsolo

CHANGELOG.md

@@ -213,6 +213,86 @@ Sections commonly used: Features, Bug fixes, Other changes.
 
 ### Bug fixes
 
+- **STARsolo's stats files as STAR writes them.** `Solo.out/Barcodes.stats`,
+  `<feature>/Features.stats` and `<feature>/UMIperCellSorted.txt` were not
+  written, SJ and Velocyto had no `Summary.csv`, and the Gene/GeneFull
+  `Summary.csv` lacked the Q30 and "Unique+Multiple" lines, printed fixed
+  decimals, and computed several values its own way. They now come from
+  STAR's counters (`SoloReadBarcodeStats`, `SoloReadFeatureStats`) and its
+  `statsOutput`:
+  - each read counts once in `Barcodes.stats`, under its final `cbMatch`
+    code (a UMI rejection replaces the barcode's); a missing adapter or UMI
+    is counted as STAR counts it;
+  - the per-feature counts come from STAR's three stages: per read
+    (`record`), at counting (oneExact and 1MM_multi rejections) and at UMI
+    collapse; Velocyto keeps only the UMIs and cells, as `countVelocyto`
+    never sums the rest;
+  - "Reads With Valid Barcodes" subtracts both barcode and feature
+    rejections; "Reads Mapped to <feature>: Unique" is the counted reads;
+    means are integer, medians the element at `nCells/2`, the UMI median
+    from every barcode's UMIs;
+  - doubles, here and in `UniqueAndMult-*.mtx` (which printed `0.50000`
+    for STAR's `0.5`), are written as C++'s `ostream` does;
+  - `CB_samTagOut` writes `Barcodes.stats`, and so does a run with
+    `--soloCBwhitelist None` (all zeros, as STAR's: it bins reads only
+    against a whitelist). That mode still writes no matrices.
+
+  Also: a barcode's exact-match count (used by oneExact and the 1MM_multi
+  posterior) now counts only reads whose UMI passed too, as STAR's
+  `addStats` does; it counted exact barcodes with a rejected UMI. Across 30
+  STARsolo configurations of the 10x v2 mouse chr19 fixture, 453 of 457
+  `Solo.out` files are byte-identical to STARsolo's; the other four are the
+  `CB_UMI_Complex` homopolymer UMIs of DIVERGENCE.md §1.6.
+
+  STAR's per-line SJ counts are ported too: a junction that is not a matrix
+  row counts nowhere, a read with no junction is `noNoFeature` (a
+  1MM_multi barcode with junctions was too), and 1MM_multi barcodes on SJ
+  records, which were dropped, are resolved by the posterior as for Gene. `CB_UMI_Complex` reads rejected
+  twice can land in a different `Barcodes.stats` bin from STAR's
+  (DIVERGENCE.md §3.4).
+
+  `CellRanger.summary.csv`, a rustar-aligner extension written whenever
+  Gene and GeneFull both ran, is now opt-in: `--soloCellRangerSummary Yes`
+  (DIVERGENCE.md §1.8). By default `Solo.out` holds only STARsolo's files.
+
+- **`--soloCellFilter TopCells N` calls cells as STARsolo.** It took the N
+  barcodes with the most UMIs. STAR keeps every barcode with at least the
+  UMIs of the barcode at index N of the descending list, the (N+1)th, so
+  ties are kept (`SoloFeature_cellFiltering.cpp:46-59`). With `TopCells 100`
+  on the 10x v2 mouse chr19 fixture STAR calls 106 cells and rustar-aligner
+  called 100; the filtered matrices are now identical.
+
+- **`--soloFeatures Velocyto` counts as STARsolo.** It used its own rule (a
+  junction means spliced, an intronic block unspliced, anything else
+  ambiguous), read by read. It now ports STAR's:
+  - each uniquely mapped read is compared with every transcript that
+    contains it on the `--soloStrand` strand (`alignToTranscriptMinOverlap`,
+    velocyto's 5-base flank, already ported for bulk use in
+    `quant/splice_status.rs`);
+  - a molecule (cell, UMI, taken from the Gene feature's readInfo) keeps the
+    transcripts all its reads share; it is spliced, unspliced or ambiguous by
+    STAR's model flags, and dropped if its transcripts span two genes
+    (`SoloFeature_countVelocyto.cpp`);
+  - all three matrices hold the same (cell, gene) rows, zeros included;
+  - `Velocyto/filtered/` is written, for the Gene feature's called cells that
+    Velocyto detected (`SoloFeature_cellFiltering.cpp:16-30`), unless there
+    are none (`:106-108`).
+
+  On the 10x v2 mouse chr19 fixture the raw and filtered matrices are now
+  byte-identical to STARsolo's (spliced/unspliced/ambiguous UMIs 369/57/32,
+  against 52/100/362 before), in 10 configurations: Unstranded, Reverse,
+  `1MM_multi`, GeneFull listed first, with `CB`/`UB` tags, `1MM_Directional`
+  and `Exact` dedup, EmptyDrops_CR, and paired-end cDNA. STAR keeps a readInfo
+  per feature (`CB`/`UB` read the first feature's, Velocyto the Gene
+  feature's), and so does rustar-aligner now. `--soloFeatures Velocyto`
+  without `Gene`, on which STAR crashes, is refused; with `--soloType
+  SmartSeq` it is refused with STAR's message.
+
+  `--soloVelocytoAmbiguous` is removed. It was a rustar-aligner extension
+  that folded ambiguous molecules into spliced, on the premise that exon-only
+  reads land in ambiguous; under STAR's model they are spliced, and ambiguous
+  means exonic in one isoform and intronic in another.
+
 - **Paired-end reads with mates of 50 bases or fewer align as STAR.** STAR
   caps the spacing of its seed search starts by the whole read, both mates
   and the spacer (`ReadAlign_mapOneRead.cpp:41-42`); we capped each mate by

DIVERGENCE.md

@@ -130,7 +130,7 @@ So when `umiL` is still 0 at construction, the loop body never runs and all four
 
 Poly-A being special in RNA does not apply here: the UMI comes from the barcode read, not the cDNA. In oligo-dT chemistries the UMI sits next to the dT, so a UMI that runs into it reads as poly-T, the one homopolymer STAR's complex path lets through. Matching the accident would count poly-T UMIs as molecules.
 
-**Impact.** Under `--soloType CB_UMI_Complex`, reads whose UMI is a non-A homopolymer are dropped here and counted by STAR. On real data these are rare and are sequencing artefacts.
+**Impact.** Under `--soloType CB_UMI_Complex`, reads whose UMI is a non-A homopolymer are dropped here and counted by STAR. On real data these are rare and are sequencing artefacts. The stats files show it too: such reads count as `noUMIhomopolymer` in `Barcodes.stats` here, and under their barcode's match (or `noNoWLmatch`) in STAR's, which also moves `Features.stats` and the "Reads With Valid Barcodes" line of `Summary.csv`. On the 10x v2 mouse chr19 fixture run as `CB_UMI_Complex`, 43 reads; the count matrices are unchanged.
 
 **Source.** `src/solo/whitelist.rs` (`check_umi`, `umi_valid_rejects_every_homopolymer`). STAR: `SoloReadBarcode.cpp:16-21`, `SoloReadBarcode_getCBandUMI.cpp:139` and `:353-354`.
 
@@ -159,6 +159,18 @@ Every STAR record matches the alignment STAR computed only once SEQ is moved by
 
 **Source.** `src/clip/mod.rs` (`clip_mate_cr4`). STAR: `ClipMate_clip.cpp:14-24, 50-57, 64-70`, `readLoad.cpp:60-61`, `ReadAlign_alignBAM.cpp:224-230`, `ReadAlign_calcCIGAR.cpp:16-22`.
 
+### 1.8 `--soloCellRangerSummary`: a CellRanger-style mapping summary
+
+A rustar-aligner extension with no STAR equivalent, off by default (`No`). With the default no output changes.
+
+**What STAR does.** STARsolo writes one `Summary.csv` per feature (`SoloFeature_statsOutput.cpp`), with reads mapped to the genome and to the feature, but no split of the reads by region.
+
+**What rustar-aligner does.** With `--soloCellRangerSummary Yes` and both `Gene` and `GeneFull` in `--soloFeatures`, it also writes `Solo.out/<feature>/CellRanger.summary.csv`: of all input reads, the fractions that map uniquely and fall in exonic, intronic or intergenic regions, or antisense to a gene, under the names CellRanger's `metrics_summary.csv` uses ("Reads Mapped Confidently to Exonic Regions" and so on). The definitions are rustar-aligner's own (exonic and intronic come from the two features' overlaps) and are not claimed to equal CellRanger's. `Summary.csv` itself is never changed.
+
+**Why.** CellRanger users look at these fractions, and the two overlaps needed are computed anyway when both features run.
+
+**Source.** `src/solo/count.rs` (`write_cellranger_summary`), `src/solo/gene.rs` (`Region`).
+
 ---
 
 ## 2. Cases where rustar-aligner outperforms STAR
@@ -248,6 +260,22 @@ On the 10k yeast PE benchmark, 4 reads differ in alignment score (AS) because ST
 
 ---
 
+### 3.4 `CB_UMI_Complex`: the `Barcodes.stats` bin of two kinds of rejected read
+
+**What STAR does.** Under `CB_UMI_Complex`, STAR checks the UMI first and then matches each cell-barcode segment to its own whitelist, in order (`SoloReadBarcode_getCBandUMI.cpp:356-420`). The first failure sets `cbMatch`, which picks the read's `Barcodes.stats` bin. Two consequences:
+- a read with one mismatch in each of two segments is `cbMatch` -12, counted as `noTooManyMM`;
+- a read with a rejected UMI keeps the UMI's code (`noNinUMI` or `noUMIhomopolymer`) even when a segment cannot be extracted.
+
+**What rustar-aligner does.** It matches the concatenated segments against the product of the segment whitelists. A read with mismatches in two segments is two mismatches from every entry, so it is counted as `noNoWLmatch`. A barcode read too short for a segment is counted as `noNoCB` before its UMI is checked.
+
+**Why.** Not chosen: per-segment matching is not ported yet. The product whitelist gives the same match as STAR for a read with at most one mismatched segment, which is every read STAR counts.
+
+**Impact.** Only the bins of reads that both tools reject: the count matrices, `Features.stats` and `Summary.csv` are unchanged. On the 10x v2 mouse chr19 fixture run as `CB_UMI_Complex`, `Barcodes.stats` differs from STAR's only by the reads of §1.6.
+
+**Source.** `src/solo/mod.rs` (`SoloBarcodeLayout::extract`), `src/solo/whitelist.rs` (`CbWhitelist::load_complex`). STAR: `SoloReadBarcode_getCBandUMI.cpp:356-420`.
+
+---
+
 ## 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.

ROADMAP.md

@@ -312,7 +312,7 @@ See [docs-old/phase17_features.md](docs-old/phase17_features.md) for sub-phase t
 
 **Prerequisite met**: position agreement >99% (SE 99.815% tie-adj, PE 99.883%). Phase unblocked 2026-06-10.
 
-**Status**: A STARsolo-faithful single-cell pipeline, Gene count matrix **byte-identical to STARsolo's**, with GeneFull/SJ/Velocyto features, CB_UMI_Simple/Complex/SmartSeq chemistries, EmptyDrops_CR + CellRanger2.2 cell calling, multi-mapper resolution, `Summary.csv`, and the per-read STARsolo SAM tags. Validated in a native three-way benchmark against STARsolo and CellRanger (see below). Remaining: the Velocyto feature's counts differ from STARsolo's, and `Barcodes.stats`/`Features.stats` are not written.
+**Status**: A STARsolo-faithful single-cell pipeline, Gene count matrix **byte-identical to STARsolo's**, with GeneFull/SJ/Velocyto features, CB_UMI_Simple/Complex/SmartSeq chemistries, EmptyDrops_CR + CellRanger2.2 cell calling, multi-mapper resolution, `Summary.csv`, and the per-read STARsolo SAM tags. Validated in a native three-way benchmark against STARsolo and CellRanger (see below).
 
 Single-cell quantification layered around the existing aligner: the cDNA read aligns through the normal SE path; a paired **barcode read** (R1 = cell barcode + UMI) is parsed, corrected against a whitelist, assigned to a gene, UMI-deduplicated, and emitted as a sparse per-cell count matrix. Target: faithful port of STARsolo (all features). See [docs-old/phase14_starsolo.md](docs-old/phase14_starsolo.md) for the full design and sub-phase tracking.
 

docs/src/content/docs/reference/cli-parameters.md

@@ -162,6 +162,7 @@ Run `rustar-aligner --help` for the full machine-generated listing.
 | `--quantGeneCoverage` | `No` | rustar-aligner extension. `Yes` writes per-gene exonic coverage to `GeneCoverage.out.bin` for DegNorm. See [DegNorm](/rustar-aligner/guides/degnorm/). |
 | `--degNormSampleId` | basename of `--outFileNamePrefix` | rustar-aligner extension. Sample name stored in `GeneCoverage.out.bin`. |
 | `--degNormCoverageFiles` | — | rustar-aligner extension. Two or more `GeneCoverage.out.bin` files: runs the DegNorm merge (no reads, no `--runMode`) and writes `DegNorm.out/`. Tuning: `--degNormIter`, `--degNormNmfIter`, `--degNormDownsampleRate`, `--degNormMinimaxCoverage`, `--degNormSkipBaselineSelection`, `--degNormBins`, `--degNormMinHighCoverage` (see the guide). |
+| `--soloCellRangerSummary` | `No` | rustar-aligner extension. `Yes` writes `Solo.out/<feature>/CellRanger.summary.csv` for `Gene` and `GeneFull` (both needed): the fractions of all reads that map uniquely in exonic, intronic or intergenic regions, or antisense to a gene, under CellRanger's metric names. STARsolo's `Summary.csv` is unchanged. |
 
 ## Two-pass mode