fix(solo): MultiGeneUMI_CR decides ownership on corrected UMIs (stacked on #174) - #175
fix(solo): MultiGeneUMI_CR decides ownership on corrected UMIs (stacked on #174)#175BenjaminDEMAILLE wants to merge 4 commits into
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>
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. |
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. |
1 participant
Stacked on #174, which is stacked on #173. Merge in that order; this rebases to one commit afterwards.
Finishes the
MultiGeneUMI_CRrule. #173 fixed the tie; this adds the condition I explicitly deferred there, and it needed the restructure I said it would.What STAR actually does
Two things, in this order (
SoloFeature_collapseUMIall.cpp:134-148, then:203-235):UMIs are corrected within each gene first, and the per-gene read totals are recorded twice — once under the raw UMI (
umiGeneMapCount0), once under the corrected one (umiGeneMapCount).Then ownership is decided on the corrected map, subject to two conditions:
The second condition exists because correction moves reads between UMIs. A gene can end up winning only because correction folded a neighbouring UMI onto it, and STAR refuses that win instead of counting it.
Why this could not be a smaller change
The generic path here filters multi-gene UMIs first and corrects afterwards. That order cannot express either condition: by the time correction runs, ownership is already decided. So
MultiGeneUMI_CRnow takes its own path through the per-cell loop, which mirrors STAR — the flag is only valid with--soloUMIdedup 1MM_CR, so no other combination is affected.cellranger_1mm_mapexposes the mappingcellranger_1mmwas already computing and discarding.Measured against CellRanger 10.0.0
Real
cellranger count, not a proxy:cellranger mkrefon the yeast reference, thencellranger counton the fixture from #172. With #165 and #173 also applied:Entries CellRanger has and we do not: 29 → 7.
#165 is a precondition for these numbers and should merge first; without it the same stack sits at 15 439, +2.17%. With it, STAR 2.7.11b on the same fixture with the same flags counts 15 124, so it sits +13 from CellRanger while this sits +5. All three agree to within a fraction of a percent. The full table is in the comment thread below.
Verification
multi_gene_umi_cr_rejects_a_winner_that_only_wins_after_correction— two UMIs one substitution apart, arranged so a gene wins after correction and loses before itmulti_gene_umi_cr_drops_a_tie_entirelyandmulti_gene_umi_cr_keeps_top_genefrom fix(solo): MultiGeneUMI_CR was inert — a tied UMI goes to nobody, not everybody #173 still pass unchangedGate: 562 lib + 26 integration tests,
cargo clippy --all-targets -- -D warnings,cargo fmt --check, all green.What is still different
7 entries CellRanger has that we do not, 19 the other way, of 13 709. I have not chased them and will not guess at a cause here.
The honest limit of this method: anything left is small enough that telling a genuine algorithmic difference from a fixture artefact needs a second, differently-shaped dataset. That belongs with the fixture work in #172 rather than with this rule.