CellRanger GEX pipeline: coverage map, stage by stage

[BenjaminDEMAILLE]

Read 10x's Gene Expression algorithm page stage by stage and checked each one against the code, so that the remaining CellRanger work is a list of decisions rather than one large PR. Everything in the "where we are" column was verified in the tree at bd/solo-cellranger-h5 (the tip of #173 → #180), not recalled.

Context for the numbers: on the 20 000-read fixture, with the five CellRanger flags, rustar counts 15 439 against CellRanger's 15 111 (+2.17%) and STAR's 15 124 (+0.09%). Correcting that figure is in the comment thread on #175/#176; the gap is open, and this table is partly an attempt to find where it comes from.

Coverage

#	CellRanger stage	Where we are	Gap
1	TSO trimmed from R2 5', poly-A from 3'	Half. --clipAdapterType CellRanger4 does the 3' poly-A trim (#148). src/clip/mod.rs:22 says the 5' TSO clip is out of scope.	5' TSO clip missing. #155 proposed a dependency for it; without one it needs an in-tree aligner for a 30 bp probe
1b	ts:i / pa:i BAM tags recording bases removed	Missing. No occurrence in the tree	Small, self-contained
2	STAR splicing-aware alignment	Done	—
2b	Exonic if ≥50% of the read overlaps an exon, else intronic, else intergenic	Differs. classify_read uses STAR's concordance rule: every aligned block must lie wholly within the gene's exons	A read half in an exon counts for CellRanger and not for us. This is a candidate cause of the count gap and worth measuring before anything else
3	MAPQ 255 for a read with one exonic locus plus non-exonic ones	Missing. No such adjustment	Changes which reads are "confidently mapped", so it moves both counts and metrics_summary.csv
4	Transcriptome alignment; introns included by default since v7.0	Differs by default. --soloFeatures defaults to Gene, exonic only. GeneFull exists and is what CellRanger's default means	A default-change decision, same shape as #176
4b	Antisense reads ignored	Done. Counted for the metrics, not toward genes	—
5	Barcode correction, posterior > 0.975 on Hamming-1 neighbours	Done (#165: cbMinP, oneExact)	—
6	UMI correction: 1MM, lesser-supported UMI corrected to the prevalent one	Done (1MM_CR)	—
7	UMI counting; conflicting gene annotations discarded unless one has maximal support; ties discarded	Done (#173, #175)	—
8	Cell calling, OrdMag: cutoff m/10 where m is the 99th percentile of the top-N; minimise (OrdMag(x)−x)²/x; search 2 to ~45 000 cells	Older variant. knee_cr22 fixes N at 3 000 and takes m/10 with no search over N	The loss-function search is missing. Precisely specified on the page, and self-contained
9	Cell calling, EmptyDrops: candidates > 500 UMIs, SGT ambient profile	Done (#156): umi_min 500, SGT, libc++ sampler	—
10	Multiplet detection across two genomes, 10th percentile in both	Missing	Barnyard experiments only. Needs a two-genome index; I would not do this without a user asking
11	Secondary analysis (PCA, clustering, UMAP, analysis/)	Missing, deliberately	Not an aligner's job. Reimplementing scanpy inside rustar is the wrong shape

What I would do, in order

1. Measure #2b before building anything. If CellRanger's ≥50%-exon rule explains a large share of the 328-count gap, that is the most valuable thing on this list and everything else is noise beside it. It is a measurement, not a PR.
2. OrdMag (feat: SJ insertion into Genome + SA at genomeGenerate #8). Fully specified on the page, no ambiguity, self-contained in called_cells, and it is the step that decides how many cells a user sees. A new --soloCellFilter method, opt-in, defaulting on 10x under the same detection as feat(solo): CellRanger behaviour by default on 10x geometry (changes default output, stacked on #175) #176 if you want it there.
3. Introns by default on 10x (Bump docker/build-push-action from 7.0.0 to 7.1.0 #4). One line in CELLRANGER_DEFAULTS, but it changes counts on every 10x run, so it needs the same sign-off as feat(solo): CellRanger behaviour by default on 10x geometry (changes default output, stacked on #175) #176 and its own measurement.
4. MAPQ 255 (Bump actions/upload-artifact from 7.0.0 to 7.0.1 #3) and ts:i/pa:i (#1b), both small.
5. 5' TSO (Output anndata files #1) only if you want an in-tree aligner for it; the dependency route was already declined.
6. Not doing: multiplet detection, secondary analysis.

Opening the OrdMag one first unless you would rather I start elsewhere. Items 3 and 5 are decisions rather than work, and I would rather have your call than guess.

---

Update: two of these were measured, and one changes the plan

The +2.17% gap is #165, not stage 2b. I suspected CellRanger's ≥50%-exon rule was the cause and said it was the first thing to measure. It is not: our rule is stricter per read (every block wholly inside the gene's exons) and we still counted more. Applying #165's cbMinP posterior threshold closes it — 15 116 against CellRanger's 15 111, +0.03%, 13 676 of 13 709 entries identical. Details and a retraction of my own bad correction are on #175.

So stage 2b drops down this list, and #165 becomes a merge-order dependency for the whole stack rather than a separate concern.

Stage 4 is not the lever I assumed either. GeneFull gives fewer counts than Gene on the fixture, 14 887 against 15 439, not more. Yeast genes are dense and overlapping, so a read that is unambiguously exonic becomes ambiguous against overlapping gene bodies and is dropped. "Introns included by default" would therefore not simply add counts, and the decision needs a dataset with real intronic signal before it is worth taking.

Revised order:

1. Measure 2b — done, ruled out.
2. OrdMag (feat: SJ insertion into Genome + SA at genomeGenerate #8) — done, solo: --soloCellFilter OrdMag, CellRanger's cell call (first item of #181, stacked on #180) #182.
3. MAPQ 255 (Bump actions/upload-artifact from 7.0.0 to 7.0.1 #3) and ts:i/pa:i (#1b), both small and additive.
4. Introns by default (Bump docker/build-push-action from 7.0.0 to 7.1.0 #4) — needs a dataset with intronic signal first, not just a sign-off.
5. 5' TSO (Output anndata files #1) — unchanged: needs an out-of-tree aligner, and the dependency route was declined in Dependency discussion: HDF5, for CellRanger's raw/filtered_feature_bc_matrix.h5 #177.
6. Not doing: multiplet detection, secondary analysis.

One metric definition also got pinned down while chasing this, recorded in DIVERGENCE.md §3.4: 10x define Sequencing Saturation as 1 - n_deduped_reads / n_reads over unique (barcode, UMI, gene) combinations. Taken literally — distinct triples before UMI correction — that gives 0.0% on this fixture, because no two reads there share an exact triple. So their numerator is the corrected molecule count, as ours is, and the 7.4%-against-12.7% disagreement is in which reads count as confidently mapped, not in the formula.

---

Update 2: measured on human data, and stage 4 is now the top item

Ran 10x's pbmc_1k_v3 (20 M read pairs, refdata-gex-GRCh38-2024-A on both sides, cellranger count 10.0.0 against the stack plus #165). Everything above was measured on a 20 000-read yeast fixture; this is the second dataset I kept saying was needed, and it moves two items.

Stage 4 (introns by default) is the whole remaining gap, not a marginal decision.

build	counts	vs CellRanger	entries identical
CellRanger 10.0.0	9 187 379	—	—
rustar --soloFeatures Gene (default)	6 382 961	−30.5%	52.9%
rustar --soloFeatures GeneFull	9 039 161	−1.6%	91.8%

On yeast, GeneFull gave fewer counts than Gene and I concluded the intron decision needed a dataset with real intronic signal before it was worth taking. It had one: on human, intronic reads are 31% of the library and the flag moves us from 30% off CellRanger to 1.6% off. --soloFeatures GeneFull belongs in CELLRANGER_DEFAULTS.

Stage 8 (OrdMag) is measured and does not help. Full numbers on #182: without it, EmptyDrops_CR calls 1 178 cells against CellRanger's 1 178; with it, 1 177. The documented search is one cell worse on the one real dataset available.

Two metric disagreements resolved as fixture artefacts. Sequencing Saturation, the one I was least sure of, reads 42.1% against CellRanger's 42.0% on human — the yeast 12.7%-against-7.4% came from a fixture where no two reads share a UMI. And the Q30 metrics plus Reads Mapped Confidently to Intergenic Regions match to the decimal.

Revised order:

1. Introns by default (stage 4) — measured, large, one line. Needs the same sign-off as feat(solo): CellRanger behaviour by default on 10x geometry (changes default output, stacked on #175) #176.
2. MAPQ 255 (stage 3): CellRanger reports 93.5% confidently mapped against our 88.3%, so this one is worth real numbers too.
3. ts:i/pa:i (stage 1b), small and additive.
4. 5' TSO (stage 1) — unchanged; the dependency route was declined in Dependency discussion: HDF5, for CellRanger's raw/filtered_feature_bc_matrix.h5 #177.
5. OrdMag (stage 8) — done in solo: --soloCellFilter OrdMag, CellRanger's cell call (first item of #181, stacked on #180) #182, and its own measurement argues against it.
6. Not doing: multiplet detection, secondary analysis.