solo: two CellRanger-matching flags are inert, and the count matrix is 8.9% off

[BenjaminDEMAILLE]

Summary

Two of the five flags STAR documents for matching CellRanger are accepted by rustar-aligner and change nothing. Under the full CellRanger-matching flag set, rustar counts 8.9% more UMIs than STAR (16 465 against 15 124), and the repository's existing differential cannot see it because it compares three matrix entries.

Why the existing test does not catch this

test/solo_cellranger_diff.py builds a fixture with two barcodes, two genes and three non-zero matrix entries. Every barcode is an exact whitelist hit, every UMI is unique, every read is unambiguous. None of the paths the CellRanger flags touch is exercised, so the test passes whatever those flags do.

I built a fixture that can fail: 20 000 reads, 200 cells drawn from the real 3M-february-2018 whitelist, 400 yeast genes, and deliberately

• 1 799 barcodes one substitution from a whitelist entry,
• 560 barcodes containing an N,
• 1 421 UMIs one substitution from another UMI on the same cell and gene,
• 720 UMIs shared between two genes in the same cell.

The measurement

STAR 2.7.11b against rustar-aligner at 3dedcb4, same fixture, same index inputs, raw Gene matrix compared entry by entry. One flag added at a time on top of defaults:

variant	identical entries	STAR counts	rustar counts
defaults	14 768 / 14 808	16 453	16 465
--soloUMIdedup 1MM_CR	14 768 / 14 808	16 453	16 465
--soloUMIfiltering MultiGeneUMI_CR	13 749 / 14 806	15 423	16 465
--soloCBmatchWLtype 1MM_multi_Nbase_pseudocounts	14 437 / 14 806	16 120	16 465
--clipAdapterType CellRanger4	14 768 / 14 808	16 453	16 465
--outFilterScoreMin 30	14 768 / 14 808	16 453	16 465
all five together	13 458 / 14 804	15 124	16 465

Read the rustar column: 16 465 in every row. No CellRanger-matching flag changes rustar's output at all, while MultiGeneUMI_CR removes 1 030 counts from STAR and 1MM_multi_Nbase_pseudocounts removes 333.

The barcode sets agree exactly — 200 barcodes on both sides, none unique to either — so this is not barcode detection. It is what happens to reads after the barcode is resolved.

The two defects

1. --soloUMIfiltering MultiGeneUMI_CR is inert. 862 matrix entries that STAR removes are kept, 1 030 counts. filter_multi_gene_umi has a MultiGeneUmiCr arm with a unit test, so the logic exists; something between the flag and the counted matrix is not connected. Worth checking against #152, which touches the neighbouring MultiGeneUMI_All.

2. --soloCBmatchWLtype 1MM_multi_Nbase_pseudocounts is inert. 276 entries, 333 counts. All five CbMatchType variants parse, and resolve_multi_cb implements the count-and-quality posterior, so again the pieces exist. The Nbase variant is the one that decides what happens to the 560 reads whose barcode contains an N, which is exactly the population this fixture inflates.

Both are the same shape as --clipAdapterType CellRanger4 before #148: a flag that validates, runs, and does nothing, so a user asking for CellRanger-matching behaviour gets a successful run and different numbers.

3. A smaller residual at defaults, which is not about CellRanger: 40 entries differ out of 14 808 and rustar counts 12 more. Separate question, listed so the 8.9% is not read as including it.

4. A validation divergence, found while setting this up: STAR refuses --soloUMIfiltering MultiGeneUMI_CR unless --soloUMIdedup 1MM_CR is also given (SOLUTION: rerun with ... --soloUMIdedup 1MM_CR). rustar accepts the combination silently.

Proposed work

One PR per defect, each with a locking test built on this fixture rather than on the three-entry one:

1. Wire MultiGeneUMI_CR through to the counted matrix.
2. Wire the Nbase pseudocount variant through.
3. Adopt STAR's validation pairing for MultiGeneUMI_CR.
4. Replace the three-entry differential with this fixture, so a regression in any of the above fails a test rather than being found by hand.

Reproducing

The generator and comparison are in the branch for item 4. Two notes for anyone rerunning it:

The fixture whitelist is the real 3.7-million-barcode v3 list, taken from a local copy rather than vendored — the lists are 26 MB and this crate publishes to crates.io, so they should stay out of the repository and the test should skip when they are absent.

STAR could only be run in the container. Natively on macOS, --soloFeatures Gene fails with Transcriptome.cpp:18 ... could not open input file /geneInfo.tab: STAR restores sjdbGTFfile from the index, so trInfoDir = P.pGe.sjdbGTFfile=="-" ? P.pGe.gDir : P.sjdbInsert.outDir takes the second branch and sjdbInsert.outDir is empty. Same root cause as the blocker noted on #158.

[BenjaminDEMAILLE]

The two inert flags in #172 turn out to be fixed by two separate PRs, and together they close almost the whole gap.

Attribution

Same 20 000-read 10x fixture, STAR 2.7.11b against each build, raw Gene matrix compared entry by entry.

--soloCBmatchWLtype 1MM_multi_Nbase_pseudocounts alone:

build	identical entries	STAR counts	rustar counts
main	14 437 / 14 806	16 120	16 465
#165	14 518 / 14 530	16 120	16 108

The pseudocounts were already wired on main; what was missing is the cbMinP posterior threshold, so nothing ever rejected a low-confidence correction and the flag could not change an outcome. #165 adds the threshold and the oneExact guard, and the flag starts working.

#165 also fixes the baseline: at default flags it takes rustar from 21 entries STAR does not have down to zero, 14 773 / 14 787 identical.

--soloUMIfiltering MultiGeneUMI_CR alone is #173 (tie goes to nobody): +1 042 counts to −9.

Both together

build	identical entries	STAR counts	rustar counts	gap
main	13 458 / 14 804 (90.9%)	15 124	16 465	+1 341 (+8.9%)
#165 only	13 532 / 14 528	15 124	16 108	+984
#173 only	13 585 / 13 979	15 124	15 414	+290
#165 + #173	13 671 / 13 708 (99.73%)	15 124	15 091	−33 (−0.22%)

Full CellRanger-matching flag set, so this is the configuration the flags exist for.

What is left

37 entries of 13 708: 28 STAR-only, 2 rustar-only, 7 present in both and differing.

One known cause, already stated on #173: STAR requires the winning gene of a multi-gene UMI to hold the top count among uncorrected UMIs as well (SoloFeature_collapseUMIall.cpp:226-232), which needs the pre-correction counts that the current code does not retain. That is the next thing to chase, and it is a restructure rather than a one-line rule change.

The 28 STAR-only entries are the shape that condition would produce — STAR counting a UMI that we drop — so I would look there first.

Note on merge order

Neither PR depends on the other, but they touch src/solo/count.rs and will conflict trivially on CHANGELOG.md. Whichever lands second needs a rebase, and the combined figure above only holds with both.