feat(solo): CellRanger behaviour by default on 10x geometry (changes default output, stacked on #175) - #176
Conversation
`--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) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…w matrix
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) <noreply@anthropic.com>
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 scverse#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 scverse#172, with scverse#165 and scverse#173 also applied: identical entries CellRanger rustar scverse#165 + scverse#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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
Correction: the CellRanger numbers in this PR do not reproduceI re-derived the whole comparison from a clean state today, both sides in one pass, and the absolute figures I published are wrong. The relative effects hold. The headline claim does not, and it was the headline, so this needs saying plainly rather than in a footnote. Everything below is from today, one fixture, one
What survives
What does not survive"within 0.03%", "15 116", and "rustar is closer to CellRanger than STAR is". All three are wrong. rustar is at +2.17%; STAR is at +0.09%. Compared directly against STAR under identical flags, rustar is +315 counts with 273 entries STAR does not have. That is a real open divergence, and I reported it as closed. How it happenedThe fixture reference was regenerated between that measurement and this one, and I compared a rustar run against a CellRanger run derived from a different state of it. The relative deltas were unaffected, which is why the per-PR effects still reproduce and I did not notice. I should have re-derived both sides in one pass before publishing an absolute figure, and from here I will. What this changes about the PRNothing about the code, and nothing about the argument for it: the flags still move the matrix from 8.96% to 2.17% off CellRanger, which is most of the gap and is what the change is for. What changes is that the remaining 2.17% is open, not closed, and the 273 entries rustar has that STAR does not are the next thing to chase. I would rather that be visible in the review than discovered later. I have edited the numbers in the PR description to match this table. |
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) <noreply@anthropic.com>
Retracting the correction aboveThe correction I posted earlier was itself wrong, and I would rather say so immediately than leave it standing. The original numbers in this PR were right. What I failed to hold constant when I re-measured was #165, which is an open PR and not in this stack's base. Re-measured today, both sides in one pass, one fixture, one
The last row reproduces the original claim exactly, including the 13 676 / 13 709 identical entries. The PR body already named the precondition — it said "with #165 and #173" — and I dropped that condition when re-measuring, then published the resulting worse number as a correction of a claim that was never wrong. So, to be unambiguous about what holds:
Practical consequence for review: #165 should merge before this stack, or the numbers in these PRs will not reproduce. I have restored the original figures in the description and made that dependency explicit rather than parenthetical. The lesson I am taking from it, since it cost you reading time twice: when a measurement disagrees with a published one, the first thing to check is what moved between them, not the published number. |
The figure I replaced this with was measured without scverse#165, whose cbMinP posterior threshold is a precondition for it. With scverse#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) <noreply@anthropic.com>
Human data: the +0.03% is a yeast artefact, and the real gap is intronsThe 20 000-read yeast fixture is the only dataset these PRs were measured on, and I flagged more than once that a second, differently-shaped one was needed. I ran it: 10x's own Raw matrix, against CellRanger
Our default is 30.5% below CellRanger on real human data. The +0.03% I measured on yeast is correct for that fixture and does not generalise at all — yeast has almost no introns, so exonic-only counting happened to agree. On human, CellRanger has included intronic reads by default since v7.0 and we do not. That is one flag:
|
| metric | CellRanger | rustar |
|---|---|---|
| Q30 Bases in Barcode | 94.2% | 94.2% |
| Q30 Bases in RNA Read | 90.0% | 90.0% |
| Q30 Bases in UMI | 92.8% | 92.8% |
| Reads Mapped Confidently to Intergenic Regions | 3.7% | 3.7% |
| Sequencing Saturation | 42.0% | 42.1% |
| Valid Barcodes | 97.3% | 97.5% |
| Reads Mapped Confidently to Intronic Regions | 31.0% | 30.6% |
| Reads Mapped Antisense to Gene | 7.9% | 7.7% |
Sequencing Saturation was the metric I said I was least sure of. On yeast it read 12.7% against CellRanger's 7.4% and I could not explain it. On human it is 42.1% against 42.0%. The yeast disagreement was a property of a fixture where no two reads share a UMI, not a defect in the formula. DIVERGENCE.md §3.4 needs updating to say that, and I will.
The metrics that remain far apart — Reads Mapped Confidently to Transcriptome 81.3% vs 55.2%, Total Genes Detected 23 893 vs 19 531, Median UMI Counts per Cell 6 353 vs 4 555 — are all downstream of the same intron difference, not separate problems.
What I think this means for the stack
Nothing here invalidates #173–#179: the rules they implement are still the rules, and the yeast numbers in their descriptions are still what that fixture produces. What changes is the headline. "Within 0.03% of CellRanger" is true of one small yeast fixture and false of real human data, and I would rather say so here than have a reviewer find it.
The concrete consequence is that --soloFeatures GeneFull belongs in CELLRANGER_DEFAULTS in #176, alongside the five flags already there. It is a one-line change and I have not made it, because it changes default output again and #176 is already asking for sign-off on that. Your call; the measurement is above.
1 participant
Stacked on #175 → #174 → #173. Merge in that order.
This changes default output behaviour on 10x runs. I am putting that in the first line rather than the last, because it is the thing to argue about and everything else in the PR is arithmetic.
The problem
A user aligns 10x data, compares the count matrix against CellRanger, and gets a successful run with different numbers. Nothing in the output points at the five flags that explain the difference; they are in STAR's
docs/STARsolo.mdand nowhere the run can see.Measured against CellRanger 10.0.0 on the 20 000-read fixture from #172, those five flags are the whole gap:
What this does
On geometry that is unambiguously 10x —
CB_UMI_Simple, a whitelist, a 16-base cell barcode, and a 10- or 12-base UMI — these default to their CellRanger values:and the run says so:
Verified end to end: with no CellRanger flags on the command line, the output is byte-for-byte what passing all five produces.
The three things that make it defensible, and the one that does not
It is escapable. A flag named on the command line always wins, including when the value asked for is STARsolo's own default.
value_sourcedistinguishes an explicit flag from a default, so--soloUMIdedup 1MM_Allgets1MM_All. A divergence you cannot turn off would be a different and much worse thing than a divergent default.It is narrow. 16 bp CB with a 10 or 12 bp UMI and a whitelist is the geometry of every 10x 3'/5' gene expression chemistry and nothing else in common use. Non-10x geometry is untouched, and there is a test for that.
It is loud. Every run it touches announces the substitution and the geometry that triggered it.
What it is not: STAR-faithful. It is a deliberate change of default output, recorded in
DIVERGENCE.md§1.3 as the largest entry in that file.CONTRIBUTING.mdsays this needs sign-off, and I am asking for it rather than assuming it.If you would rather have the detection without the action — same geometry check, same message, no flags changed — that is a two-line change to this PR and I will make it on request. It was my own recommendation; the project owner asked for the default instead, and I would rather implement what was asked and let you rule on it with the code in front of you than argue it in the abstract.
Verification
ten_x_geometry_defaults_to_cellranger_behaviour— all five, from a bare 10x command linean_explicit_flag_beats_the_10x_default— including asking for STARsolo's own defaultnon_10x_geometry_keeps_starsolo_defaults— a 12-base barcode changes nothingten_x_lengths_without_a_whitelist_keep_starsolo_defaults— no whitelist, no overrideGate: 566 lib + 26 integration tests,
cargo clippy --all-targets -- -D warnings,cargo fmt --check, all green.