Read-budget (N) benchmark — 2,000 reads is enough (default cut ~100× from 200k)

[lhqing]

A developer lesson. How many reads does probe → resolve actually need to call a library's
chemistry correctly, and what do size and speed cost for reads beyond that? The default was
DEFAULT_MAX_READS = 200 000, set arbitrarily. This measured whether it should be smaller. It should:
the default is now 2 000 (PROBE_VERSION 2026.7.2, issue #63).

Related: issue #63 (the fingerprint that made this measurable), issue #66 (the CPU/MEM profile).

TL;DR

• Across the 13 C. elegans re-mapping datasets, the resolved chemistry is invariant from N = 1 000
  to N = 200 000 — every dataset calls its production chemistry at every N, exit 0, not one flip.
• New default: N = 2 000 — a deliberate 2× cushion over the observed ≤1 000 floor, ~100×
  smaller and ~faster than the old 200 000, still robust to the caveats below.
• Reads beyond the floor buy nothing for the chemistry decision — the winner is unchanged; only
  the (advisory) score drifts slightly.

Why this matters

N (the bounded read budget) sets three things at once:

1. Correctness — too few reads and the whitelist-hit / geometry signals get noisy, and a v2/v3
   tie or a barcode-absent refusal can appear.
2. Fingerprint size — a Preflight: portable fingerprint package (seqforge preflight) — subsampled FASTQs + info files + pinned content-address #63 package slices N reads per file; the slice is ~linear in N.
3. Speed + memory — probe and resolve are ~linear in N (issue perf: profile of the deterministic compile path — probe + resolve dominate and scale linearly with N #66).

So the smallest N that is still correct is a strict win on the other two axes. 200 000 was never
measured; this was that measurement.

Method

• Datasets: the 13 worm scRNA/snRNA studies under 260612-worm/mapping/<ACC>, one representative
  run (an SRX's FASTQs) each.
• Grid: N ∈ {1 000, 5 000, 10 000, 50 000, 100 000, 200 000}.
• For each (dataset, N): the real resolve_runs(max_reads=N) → top chemistry, advisory score, rung,
  exit code; and preflight(reads=N) → the slice's package size.
• Reference = each dataset's production manifest.yaml (the full-pipeline ground truth), not
  the benchmark itself. A smaller N "matches" if it yields the same top chemistry and exit 0.
• Real shipped whitelists (3M-february-2018, …), run in parallel across 48 cores. No harvest — the
  read budget never touches prose extraction.

Results — the call is N-invariant, and the slice is ~100× smaller

Dataset	ref chemistry (production)	invariant 1k→200k?	200k slice	1k slice	saving
GSE126954	10x-3p-gex-v2	yes	206.6 MB	1.08 MB	~191×
GSE136049	10x-3p-gex-v2	yes	26.4 MB	136 KB	~193×
GSE208154	10x-3p-gex-v3	yes	216.8 MB	1.14 MB	~190×
GSE208229	10x-3p-gex-v3	yes	15.6 MB	81 KB	~193×
GSE226715	10x-3p-gex-v3	yes	58.7 MB	309 KB	~189×
GSE229022	10x-3p-gex-v3	yes	33.8 MB	172 KB	~196×
GSE234962	10x-3p-gex-v3	yes	31.1 MB	159 KB	~196×
GSE256266	10x-3p-gex-v3	yes	47.5 MB	242 KB	~195×
GSE274290	bd-rhapsody-wta	yes	54.5 MB	278 KB	~196×
GSE310667	10x-3p-gex-v3	yes	29.0 MB	147 KB	~197×
PRJNA1027859	10x-3p-gex-v3	yes	13.8 MB	71 KB	~193×
PRJNA1195922	10x-3p-gex-v3	yes	14.6 MB	75 KB	~195×
PRJNA658829	10x-3p-gex-v3	yes	31.0 MB	148 KB	~208×

Families: 10× 10x-3p-gex-v3, 2× 10x-3p-gex-v2 (GSE126954, GSE136049), 1× bd-rhapsody-wta
(GSE274290). Every one is production-consistent at every N. The chosen default, N = 2 000, sits
just above the 1k column — ~2× its size, still ~100× under 200k.

(GSE126954 is 10x-v2 from the bytes despite being a "sci-RNA-seq" study — the production manifest
confirms v2; the benchmark is scored against that manifest, not the paper's wording.)

The speed angle (from #66)

probe and resolve scale ~linearly with N, so the read cut is also a ~100× compute win. One caution
on the raw profiler numbers: pyinstrument (statistical sampling) inflates tight pure-Python loops,
so its per-function self-times over-state the absolute cost — but its ranking of the hotspots is
reliable. Two things came out of it:

• per_cycle_composition (top of the probe tree) was vectorized — 8–13× faster, byte-for-byte
  identical to the old loop (it feeds the observation hash, so identity was the bar). Merged.
• The resolve-side whitelist/anchor/motif scans are the next hotspots, but they decide the chemistry
  — vectorizing them byte-identically is reviewed-follow-up work, not merged.

Net: lowering the default is a free ~100× before any code change; the vectorization compounds it on the
explicit large---max-reads path.

Recommendation → shipped

Default N = 2 000, DEFAULT_MAX_READS, applied globally (probe, resolve, run,
manifest fill, remote probe). A caller who wants to read more of a full-size FASTQ passes a larger
--max-reads (the explicit opt-in; run/manifest fill gained the knob). PROBE_VERSION → 2026.7.2.

Caveats

• Head-slice, not a random sample. We read the first N reads (you cannot random-sample without
  reading the whole file). Fine for chemistry; the theoretical failure mode is a file whose first N
  reads are unrepresentative. Not seen here. (An HQ/quality-ranked sample was considered and dropped —
  it needs 10× the reads for no gain the benchmark could find.)
• One representative run per dataset. A dataset with heterogeneous runs could hide a run that needs
  more reads. The 2× cushion is for exactly this; the real fingerprint package slices all runs.
• This re-hashes pinned manifests, deliberately. Lowering DEFAULT_MAX_READS changes observation
  values → content addresses, so it is a PROBE_VERSION bump and a re-hash. That was the intended,
  accepted step this benchmark justified — not a silent change.

Appendix — reproduce

bench_n.py (representative-run sweep) and analyze_bench.py (the tables) under the seqforge env.
Both operate on real FASTQs bounded by N; neither reads a whole file.