solo: EmptyDrops_CR on CellRanger's actual statistics (SGT ambient profile, libc++ sampler) - #156
Open
BenjaminDEMAILLE wants to merge 3 commits into
Open
Conversation
…_distribution STARsolo's `EmptyDrops_CR` rescue draws from `std::mt19937`, converts to doubles with `std::generate_canonical<double, 53>`, and picks categories with `std::discrete_distribution`. Two of those three are implementation-defined in the parts that matter: the standard fixes mt19937's output but not how `generate_canonical` consumes it, and says nothing about how `discrete_distribution` maps a uniform onto categories. So porting "the algorithm" is not enough — it has to be libc++'s algorithm, because that is what STAR is built against and where its numbers come from. libc++ accumulates two 32-bit draws in *ascending* significance and divides by 2^64; a most-significant-first accumulation, or one draw scaled to 53 bits, both give perfectly good uniforms and neither reproduces STAR. Every expected value in the tests came out of a C++ program compiled against the real libc++ and run, not from reading its source. `tests/libcxx_oracle.cpp` is that program, kept so the values can be regenerated rather than trusted. `generate_canonical` is compared as bit patterns, since a difference in the last place changes which category a sample lands in. Not yet wired into the EmptyDrops path. `solo::count` samples with a `SplitMix64` stream under a comment calling it "WeightedIndex-equivalent; empirically byte-identical EmptyDrops cell calls" — a claim that cannot hold in general, since two unrelated generators cannot agree on an arbitrary number of draws. It is true of whatever was checked and unknown elsewhere. Replacing it moves cell calls, so it belongs in its own change with the solo differential run against it. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Two approximations in the CellRanger cell-calling path are replaced by what CellRanger and STAR actually compute. Both move cell calls, which is the point: the previous numbers were plausible rather than right. The ambient profile is now smoothed with Simple Good-Turing (Gadsby & Sampson, via Elworthy's implementation, which is what STAR vendors). The ambient counts come from a small sample of empty droplets, so a gene seen twice there is not twice as likely as one seen once, and a gene seen zero times is not impossible — it is one the sample was too small to show. SGT fits the frequency spectrum and reserves mass for the unseen from the singleton rate, then smooths the rest along a log-log line. What was here before had the right shape and the wrong numbers: it reserved mass the same way but distributed the remainder in proportion to raw counts, with no smoothing at all. The Monte-Carlo null is now drawn with libc++'s `std::mt19937` and `std::discrete_distribution`, seeded `19760110 * (isim + 1)` per simulation, as STAR seeds it. The previous sampler was a SplitMix64 stream under a comment calling it "WeightedIndex-equivalent; empirically byte-identical EmptyDrops cell calls" — a claim that cannot hold in general, since two unrelated generators cannot agree over an arbitrary number of draws. The libc++ types were ported and checked against real libc++ in the previous commit on this branch; this wires them in. One generator per simulation, no shared state, so the walks still run in any order on any number of threads and give the same p-values. D17 comes with it: STAR leaves `PZero` uninitialised when the spectrum has fewer than five distinct frequencies and `analyse()` bails, so it reads whatever the stack held. Here it is zero from construction, which is what "no basis for reserving unseen mass" means. Recorded in docs-old/dev/divergences.md.
Section 1.2, in the What STAR does / What rustar-aligner does / Why / Impact / Source format CONTRIBUTING.md asks for, replacing the docs-old file the earlier version of this work carried.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
--soloCellFilter EmptyDrops_CRcomputes what CellRanger computes: a Simple-Good-Turing ambient profile, sampled with libc++'s generator.What changed
The ambient profile is smoothed with Simple Good-Turing (
src/solo/sgt.rs), the method CellRanger uses and STAR vendors as Elworthy'ssgt.h. The ambient counts come from a small sample of empty droplets, so a gene seen twice there is not twice as likely as one seen once, and a gene seen zero times is not impossible — it is one the sample was too small to show. SGT reserves mass for the unseen from the singleton rate and smooths the rest along a fitted log-log line.What this replaces had the right shape and the wrong numbers: it reserved unseen mass the same way, then distributed the remainder in proportion to raw counts with no smoothing at all.
The Monte-Carlo null is drawn with libc++'s
std::mt19937andstd::discrete_distribution(src/solo/libcxx_rng.rs), seeded19760110 * (isim + 1)per simulation as STAR seeds it. The standard fixes mt19937's output but not howgenerate_canonicalconsumes it, nor howdiscrete_distributionmaps a uniform draw onto categories, so reproducing STAR's numbers needs libc++'s algorithm specifically, not a correct categorical sampler.The module is wired into the EmptyDrops path in this same PR — it is not a module landed ahead of its use.
Why
solo::countsampled with aSplitMix64stream under a comment calling it "WeightedIndex-equivalent; empirically byte-identical EmptyDrops cell calls". That claim cannot hold in general: two unrelated generators cannot agree over an arbitrary number of draws. It was true of whatever cases were checked and unknown everywhere else.Divergence
DIVERGENCE.md§1.2. STAR leavesPZerouninitialised when the ambient spectrum has fewer than five distinct frequencies andanalyse()bails, so it reads whatever the stack held. Here it is zero from construction. Degenerate inputs only, but on those an uninitialised read can place arbitrary mass on unseen genes.Verification
Every expected value in the RNG tests was produced by compiling a C++ program against real libc++ (
clang++ -stdlib=libc++) and printing the results, not derived from reading its source.tests/libcxx_oracle.cppis in-tree so they can be regenerated.generate_canonicalis asserted on bit patterns rather than decimals: a difference in the last place changes which category a sample lands in.SGT tests cover a well-formed spectrum summing to one, monotonicity in the observation count, the absent-count case, the no-singletons case, and the D17 guard.
Gate: 570 lib + 26 integration tests,
cargo clippy --all-targets -- -D warnings,cargo fmt --check, MSRV 1.89 — all green.test/solo_diff_docker.sh(one run, Linux container with STAR 2.7.11b) passes:Read that for what it is. The harness exercises the raw matrix, which these changes do not touch, so it confirms nothing was broken on the way past — not that the cell calls are right. Both changes move
filtered/membership by construction, and the harness has no EmptyDrops fixture to compare against. Validating the cell calls themselves needs a dataset with a known CellRangerfiltered/result, which the harness does not currently carry.Split out of #152 following the one-theme rule in CONTRIBUTING.md.