Align pipeline stops scaling past ~4-8 threads: decode and write are each a single thread

[BenjaminDEMAILLE]

Measurement

Synthetic 20 Mb genome, 2 M single-end 100 bp reads, --outSAMtype BAM Unsorted, macOS aarch64,
release build (with the batch-size fix from #222 already in, so this is the remaining ceiling):

--runThreadN	wall	user CPU	speedup
1	8.81s	11.67s	1.00x
2	4.82s	12.36s	1.83x
4	2.80s	13.49s	3.15x
8	2.43s	15.79s	3.63x
12	2.44s	17.19s	3.61x

Throughput stops improving after ~8 threads while user CPU keeps climbing and sys roughly doubles
(1.96s → 4.45s). On a 32- or 64-core machine the extra cores are being paid for and not used.

Where it goes

Same input, 8 threads, output disabled:

output	wall
BAM Unsorted	2.43s
None	1.60s

So the writer stage alone is ~34% of wall at 8 threads, and it is one thread: SAM record → BAM
encode → BGZF deflate, serially, for the whole run.

Gzipped input costs on top of that (8 threads, BAM out): 2.74s from .fq.gz vs 2.43s from plain
FASTQ, i.e. ~13%, again on a single decode thread.

Why

The align pipelines are decode(1 thread) → align(N rayon workers) → write(1 thread), joined by
bounded channels. The two ends are serial, so Amdahl caps the whole thing regardless of
--runThreadN. The profile agrees: semaphore_wait_trap + __psynch_cvwait were ~16% of all
samples, i.e. workers parked waiting on the ends.

Options, in payoff order

1. Multithreaded BGZF writer. noodles-bgzf 0.50+ runs its multithreaded writer on its own
   local pool with a supported with_worker_count, so it no longer competes with rayon's global
   pool. Biggest single lever, and the version bump is already in deps: noodles 0.113 -> 0.115 and noodles-bgzf 0.49 -> 0.51 (local thread pools, with_worker_count) #207 / chore(deps): bump noodles 0.113 -> 0.115, noodles-bgzf 0.49 -> 0.51 #211. Needs care: BGZF
   block boundaries must stay byte-identical, or the change must be documented.
2. Parallel or chunked decode. Related to add praseq as an option for input #95 (paraseq, built for exactly this) and Parallel inputs ({C,V,,}BQ or parallel readers for fastq) #97
   (parallel readers). For .gz specifically, note that a BGZF-compressed FASTQ can be inflated in
   parallel with the noodles-bgzf we already depend on, with no new dependency.
3. Only then consider reshaping the stage handoff itself (work-stealing over batches rather
   than fixed producer/consumer threads). Do not start here: it is the largest change and the
   smallest measured share.

Note

The numbers above are from a synthetic genome on a laptop; the shape (serial ends capping the
scaling) is architectural, but the exact crossover point should be re-measured on a real genome and
a real machine before sizing any of the fixes. That needs the benchmark harness from the
dependency sweep.

[BenjaminDEMAILLE]

Measured the parallel-gzip option for the decode stage, since it now exists as a proper crate.

rapidgzip 1.2.1 (2026-04-01, MIT OR Apache-2.0) is a safe Rust wrapper over rapidgzip-sys,
which vendors ~20 MB of C++ (rapidgzip + zlib + rpmalloc) and builds it with CMake. API is
ReaderBuilder::new().parallelism(n).open(path) giving a Read + Seek, plus index import/export.
I built it from crates.io on macOS aarch64 (AppleClang, cmake 4.4.2) and benchmarked it against the
flate2/zlib-rs path this project already uses.

Input: the 2 M read synthetic FASTQ from the measurements above, 427 MB decompressed.

gzip -6 (the default anyone's data is in), 73 MB:

decoder	wall	throughput
flate2 + zlib-rs, 1 thread	0.33s	1300 MB/s
rapidgzip parallelism=1	0.78s	548 MB/s
rapidgzip parallelism=2	0.68s	624 MB/s
rapidgzip parallelism=4	0.28s	1524 MB/s
rapidgzip parallelism=8	0.25s	1705 MB/s

gzip -1, 79 MB (favourable case, fewer back-references):

decoder	wall	throughput
flate2 + zlib-rs, 1 thread	0.52s	819 MB/s
rapidgzip parallelism=1	0.80s	537 MB/s
rapidgzip parallelism=4	0.25s	1719 MB/s
rapidgzip parallelism=12	0.33s	1284 MB/s (regresses)

Reading

On default-compressed input it takes 8 rapidgzip threads to beat one zlib-rs thread by 1.3x,
and single-threaded it is 2.4x slower. zlib-rs is simply already fast. For an aligner those 8
threads are exactly the ones we want aligning, so the trade is bad unless cores are idle.

Sizing the prize from the earlier numbers: gzipped input cost +0.31s of 2.74s wall (11%) at
8 threads, and the decode already overlaps alignment on its own thread. So perfect parallel inflate
buys at most ~11% here, whereas the serial writer stage is 34%.

Cost side, for a crate published to crates.io and built on five platforms including Windows:
a CMake and C++17 toolchain requirement for every downstream build, 20 MB of vendored C++, a
second allocator (rpmalloc) inside a process that installs mimalloc globally, and a young crate
(first release 2026-03-25, ~2.3k downloads).

Recommendation

Not now. Cheaper things first, in this order:

1. The multithreaded BGZF writer (deps: noodles 0.113 -> 0.115 and noodles-bgzf 0.49 -> 0.51 (local thread pools, with_worker_count) #207) — 3x the payoff of input decode.
2. If the input is BGZF, inflate it in parallel with noodles-bgzf, already a dependency, no new
   surface at all.
3. --readFilesCommand already lets a user pipe igzip/rapidgzip/pigz in today at zero cost
   to us, which covers the person who really is input-bound.

Caveat on the numbers: synthetic reads (random ACGT, constant quality) compress unlike real FASTQ,
and this is one laptop. The ordering (zlib-rs saturating, writer dominating) is what matters, not
the absolute MB/s.

[BenjaminDEMAILLE]

Two corrections and a pointer, from the COMBINE-lab work published this month.

1. The rapidgzip numbers above were the wrong crate. They measured rapidgzip on crates.io (an FFI wrapper over C++, built with CMake). The pure-Rust rapidgzip-core 0.3 is a different project and much better: 1581 MB/s on one thread (no penalty against flate2's 1552) and 3356 MB/s on eight, with no C toolchain because its inflate backend is zlib-rs, the same one flate2 uses here. #224 and #225 are updated.

2. The conclusion for this issue is unchanged, and now quantified. Decode is not what caps us. We consume decompressed FASTQ at ~140 MB/s while one inflate thread supplies ~1550 MB/s, so the decode stage has roughly 10x headroom. Parallel decode cannot help until the aligner itself is ~10x faster. The serial writer stage, at ~34% of wall, remains the actual ceiling.

3. There is now a crate for the split problem itself, once it does become live: thread-broker 0.1 (BSD-3-Clause, from piscem-rs), described as closed-loop scheduling of a fixed thread budget between a decompressor and the work that consumes it. It solves for the decode share from measured busy time on each side (d* = N·P/(P+C)) instead of fixing a ratio, and is advisory: if the broker dies the run finishes on the last split.

That is directly the shape of option 3 in this issue, and it is worth remembering when the writer work lands, because the same argument applies to a writer budget: today --runThreadN sizes only the align pool, and the decode and write stages get one thread each by construction. If we ever make either of them parallel, the split between the three should be measured, not guessed.