solo: --soloOutLayout CellRanger, CellRanger's output shape (changes output paths on 10x, stacked on #176) - #178
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>
7425639 to
816d4e8
Compare
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>
Writes the same numbers where `cellranger count` writes them:
outs/{raw,filtered}_feature_bc_matrix/, gzipped, a -1 GEM-well suffix on
every barcode, one raw column per observed barcode, and no per-feature
subdirectory when a single feature is requested.
The layout implies --soloOutGzip yes, --soloOutRawBarcodes Observed and an
outs/ output directory; each is still overridable on the command line. On
10x geometry it is the default, alongside the five flags from the previous
commit, so it changes where output files are written on those runs.
Counts are untouched: 13 959 entries and 15 439 counts in both layouts on
the 20 000-read fixture, compared entry by entry. Against a real
cellranger count 10.0.0 run on the same fixture the raw barcode sets match
exactly, 200 of 200.
Recorded in DIVERGENCE.md 3.3.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
816d4e8 to
df2253a
Compare
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 #176 → #175 → #174 → #173. Merge in that order.
This changes where output files are written on 10x runs. First line, not last, for the same reason as #176.
Part of #172 (matching CellRanger), and the output half of #105. #174 got the columns of the raw matrix right; this gets the rest of the shape.
The problem
After #165 and #173–#176 the numbers are within 0.03% of CellRanger, and nothing that reads CellRanger output can open them.
scanpy.read_10x_mtx, Seurat'sRead10X, and every in-house loader key onraw_feature_bc_matrix//filtered_feature_bc_matrix/, on gzipped files, and on the-1suffix. We wroteSolo.out/Gene/raw/, uncompressed, bare barcodes. Right answer, unreadable by the tools it is an answer for.What this does
--soloOutLayout CellRangerwrites the same numbers in the shapecellranger countproduces:Solo.out/Gene/raw/outs/raw_feature_bc_matrix/Solo.out/Gene/filtered/outs/filtered_feature_bc_matrix/matrix.mtx,barcodes.tsv,features.tsv.gzAAAGAACGTCGAGCTCAAAGAACGTCGAGCTC-1It implies
--soloOutGzip yes,--soloOutRawBarcodes Observedand anouts/output directory. Each is still overridable, so the layout can be adopted piecemeal:--soloOutLayout CellRanger --soloOutGzip nogives the directory names without the compression.The per-feature level (
Gene/,GeneFull/) is dropped when exactly one feature is requested, which is what puts the matrices directly underouts/. With several features it is kept — CellRanger has no equivalent, and collapsing them would have each feature silently overwrite the last.On 10x geometry it is the default, via the same detection as #176.
Measured
Counts are untouched. On the 20 000-read fixture from #172, compared entry by entry after stripping the
-1suffix and matching on (gene id, barcode):7127 3686400 139597127 200 13959Identical, entry by entry.
Against a real
cellranger count10.0.0 run on the same fixture: same three filenames in both directories, same-1suffix, and the raw barcode sets match exactly, 200 of 200. Before this PR that overlap was 0.What it does not do
raw_feature_bc_matrix.h5,filtered_feature_bc_matrix.h5,metrics_summary.csv,molecule_info.h5,web_summary.html,cloupe.cloupeandanalysis/are not here. The.h5files need an HDF5 dependency, whichCONTRIBUTING.mdsays to raise in an issue before a PR, so I have: #177.metrics_summary.csvis real work and is its own PR. The last three I am not proposing at all:cloupeis an undocumented proprietary format, andanalysis/is a PCA/clustering/UMAP pipeline that does not belong inside an aligner.I also deliberately do not copy CellRanger's
%metadata_json:header line intomatrix.mtx. It carries a"software_version"field, and stamping a version string into a data file makes byte comparison across releases meaningless. Our%comment line stays.Non-STAR flag and a default change, both needing sign-off
--soloOutLayoutdoes not exist in STAR. It is inDIVERGENCE.md§3.3 with the rationale, and it is off by default on everything except 10x geometry.The size of the default change, stated plainly: it broke 7 existing integration tests, all of which assert STARsolo paths on fixtures that use a 16 bp CB and a 12 bp UMI. Each now names
--soloOutLayout STARsolo. That is the honest measure of the blast radius, and it is also a demonstration that the flag escapes the default. If you would rather have the layout opt-in only — one line removed fromCELLRANGER_DEFAULTS, and those 7 tests revert — say so and I will make that change.Verification
test_solo_out_layout_cellranger— end to end from a bare 10x command line:outs/raw_feature_bc_matrix/*.gzexists, noSolo.out/, barcodes carry-1, and the matrix holds the same 2 deduped moleculestest_starsolo_gene_matrixasserts under the STARsolo layoutbarcodes_carry_the_gem_well_suffix_only_under_the_cellranger_layout— both barcode writers, both suffixescellranger_layout_implies_gzip_observed_barcodes_and_outs_diran_explicit_flag_beats_what_the_cellranger_layout_impliesnon_10x_geometry_keeps_the_starsolo_layoutflate2appears under[dev-dependencies]so the integration test can read the gzipped files. It is already a runtime dependency at the same version and with the same features, so nothing new enters the tree.Gate: 570 lib + 27 integration tests,
cargo clippy --all-targets -- -D warnings,cargo fmt --check,cargo +1.89 check --all-targets, all green.