solo: cbMinP posterior threshold, oneExact guard, adapter-anchored geometry (replaces #150)

[BenjaminDEMAILLE]

STARsolo barcode resolution: the cbMinP posterior threshold, the oneExact guard, --soloCBtype, adapter-anchored CB_UMI_Complex geometry, and --soloOutFormatFeaturesGeneField3.

Replaces #150, rebased onto current main. Opening a new PR rather than force-pushing that branch, so the review history there stays readable. Every objection @Psy-Fer raised on #150 is addressed below, including one that was a real correctness bug.

What #150 got wrong, and what changed

The adapter search took the leftmost acceptable match; STAR takes the best-scoring one. @Psy-Fer's diagnosis was exact, including why the tests missed it. localAlignHammingDist (SequenceFuns.cpp:341) scores every offset and keeps the fewest-mismatch one:

uint32 distBest = query.size();
for (ii = 0 .. text.size()-query.size()) {
    dist1 = <mismatches at offset ii>;
    if (dist1 < distBest) { distBest = dist1; pos = ii; }
}

The old code returned the first offset inside the mismatch budget, which is the opposite rule: an early approximate match beat a later exact one. This is not cosmetic. CB_UMI_Complex measures its barcode positions from the adapter, so the wrong offset shifts the whole barcode and changes what the read is counted as.

Reading that function properly surfaced two more details the old code also had wrong:

• Ties go to the leftmost of the best, because STAR replaces its incumbent only on a strictly smaller distance.
• An N in the adapter is never a mismatch. It is a wildcard in the query, not a base that fails to match, and charging it spends budget STAR does not spend.

STAR also seeds distBest with the query length, so an offset mismatching everywhere can never win. That only shows when the budget is as large as the adapter, where it keeps the answer None instead of an arbitrary offset. Reproduced.

find_adapter_takes_the_best_scoring_match_not_the_leftmost is the discriminating case the old tests lacked: a one-mismatch match at offset 0 against an exact one at offset 5. Leftmost gives 0, STAR gives 5. Plus find_adapter_breaks_ties_leftmost and an_n_in_the_adapter_is_not_a_mismatch.

--soloChemistry and --soloCBwhitelistDir are gone. Both objections were right and they were the same mistake: neither flag exists in STAR. They came from the STAR-rs donor, which invented them, and #150 presented them as parity work without saying so. Per CONTRIBUTING.md a non-STAR flag needs explicit sign-off and must never be presented as faithfulness, so rather than relabel them I reverted them. If chemistry presets are wanted they should be their own proposal, argued as an addition.

DIVERGENCE.md had no entry, and needed one. Added §1.2, with the STAR source checked rather than recalled: rustar rejects all four homopolymer UMIs, STAR under CB_UMI_Complex rejects only poly-A. The reason is initialisation order, not a threshold. The four packed constants are built in the SoloReadBarcode constructor (SoloReadBarcode.cpp:16-21) from pSolo.umiL, but that path does not learn its UMI length until it has extracted one (SoloReadBarcode_getCBandUMI.cpp:353-354), so all four constants stay 0. Poly-A packs to 0 and is still caught; poly-C/G/T pack non-zero and match nothing (:139). This is not new behaviour in this PR, rustar always rejected all four; what is new is that it is deliberate, locked by umi_valid_rejects_every_homopolymer, and written down.

The rest of the diff

cbMinP = 0.975 and the oneExact guard. STAR accepts a multi-match CB correction only when the posterior clears 0.975 (SoloReadBarcode_getCBandUMI.cpp, CB_MIN_P), and only when no exact whitelist match was available. The posterior itself already existed in resolve_multi_cb; the threshold and the guard did not, so a low-confidence correction was accepted where STAR drops the read.

--soloCBtype. Accepts Sequence (default) and String.

Adapter-anchored CB_UMI_Complex geometry, positions measured from the adapter match rather than from the read start, which is what --soloAdapterSequence and --soloAdapterMismatchesNmax are for.

--soloOutFormatFeaturesGeneField3, the third column of features.tsv, with - suppressing it.

One thing to flag rather than hide

CONTRIBUTING.md says one theme per PR. Five of the six commits are barcode and UMI resolution; --soloOutFormatFeaturesGeneField3 is features.tsv output, written by the same code path but arguably a different area. I have kept it here rather than split it because it is 51 lines and touches no shared logic, but say the word and I will lift it into its own PR.

Verification

Unit tests in src/solo/mod.rs:

• find_adapter_takes_the_best_scoring_match_not_the_leftmost
• find_adapter_breaks_ties_leftmost
• an_n_in_the_adapter_is_not_a_mismatch

In src/solo/whitelist.rs:

• umi_valid_rejects_every_homopolymer — all four bases, plus a UMI that merely starts with a run, plus N still failing on its own grounds

In src/solo/count.rs:

• cb_posterior_below_min_p_is_rejected
• a_cb_never_seen_exactly_is_refused_without_pseudocounts

Gate on this branch, at the commit pushed: 570 lib + 26 integration tests, cargo clippy --all-targets -- -D warnings, cargo fmt --check, MSRV 1.89 — all green.

No new dependency.

[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.