[joshfactorial]

Closes #870. Supersedes #873, which can be closed: this PR rewrites the same docs page and keeps its measured numbers.

Semver: MINOR. This adds a config key, and seeded runs now write FASTQ records in a different order with the same content. Existing configs still load and run. Treat it as MAJOR instead if record order counts as public API; the release notes need an upgrade line either way.

Why

gen-reads wrote reads in genome order, which no sequencer produces. On a 21x T2T-CHM13 simulation, bwa-mem2 peaked at 166 GB on the reads as written and 90 GB on the same reads shuffled, with matching flagstat (#870, jobs 23383724 and 23394381). The old in-memory shuffle was removed in v1.5.0 because it held every read at once.

Change

• New key shuffle_fastq, default true; false restores genome order byte for byte.

• eidolon-core/src/file_tools/fastq_shuffle.rs is a two-pass external shuffle in the combine step:

  1. Scatter each record, both mates together, into one of K bucket files chosen at random. K is sized from the input, at about 256 MB per bucket, capped at 256 to stay under open-file limits. Each temp input is deleted as it is read, so temp disk stays near one copy.
  2. Permute each bucket in memory and write it at the usual compression level. The buckets are then concatenated, since gzip members concatenate.

  Pass 2 runs only as many buckets at once as fit in 4 GiB, judged by the largest bucket. Bucket files use libdeflate's fastest level (BlockGzWriter::with_level).

• Determinism: the scatter draws from one RNG child of the run (SHUFFLE_RNG_CHILD), and each bucket's permutation from a child of that. Thread scheduling only changes which bucket runs when.

• Isolation: reads, the truth VCF and the golden BAM do not depend on the shuffle.

• Tests: determinism.rs::ordinary_reads_are_emitted_in_reference_file_order now pins the shuffle_fastq: false path. That order is also the input order the shuffle's determinism rests on. GenReadsConfig gains shuffle_fastq.

• Docs: the shuffling page is rewritten for the built-in shuffle (the seqkit recipe stays, for older output); the template documents the key; the parallel-processing page links to it.

Evidence

Unit tests (fastq_shuffle.rs, 8):

• 20,000 records give an exact permutation with mates in register.
• Adjacent records within 100 indexes: < 0.03, against ~0.01 expected for random order and 1.0 for input order.
• Byte-identical output at 1 and 8 threads; a different seed gives a different order.
• k = 1 and k = 64 both permute and mix. The tests force k explicitly, because at test sizes bucket_count would always choose 1 and leave the scatter untested.
• Single-ended and empty inputs work; unequal mate files are an error.
• Bucket count scales with input and is capped; pass-2 thread count fits the memory budget.

Integration (tests/fastq_shuffle_output.rs, gen-reads on H1N1):

• Adjacent reads on the same contig within 200 bp: 0.046 shuffled against 0.996 with shuffle_fastq: false.
• The read multiset is identical either way, and mates stay in register.
• The truth VCF and the golden BAM records are byte-identical with the shuffle on and off.
• Shuffled R1 and R2 are byte-identical at 1 and 4 threads.

Mutation checks (each confirmed applied; source restored and verified with cmp):

mutant	result
skip the in-bucket permutation	2 unit + integration fail
R2 written in a different order from R1	2 unit + integration fail (rerun after the storage change: still caught)
bucket RNG seeded from scheduling order	thread-count unit + integration fail
default shuffle_fastq: false	integration fails

Not catchable from output: sending every record to one bucket still yields a correct shuffle and only loses the memory bound. bucket_count and shuffle_threads are unit-tested instead.

Cost (local, CHM13 chr21, 20×, ~3M pairs, 8 threads; off → on):

• wall time 99.4 s → 112.9 s;
• peak memory 0.46 GB → 2.34 GB, down from 2.93 GB before records were stored contiguously;
• R1 gzip 230.8 MB → 372.1 MB, +61%. Genome-ordered reads overlap their neighbors and gzip exploits that; shuffled reads, like real FASTQs, compress less well. This is stated on the docs page.

The workspace suite passes (65 test binaries, --no-fail-fast), along with clippy -D warnings and fmt.

Not verified

• Genome scale. Not run at human or CHM13 scale. The 4 GiB budget and the 256-bucket cap are reasoned, not measured there. At 30× human, buckets would be ~0.9 GB, so 4 run at a time.
• Cancer merges (gen-cancer-reads) concatenate the normal and tumor passes, each now shuffled. The merged FASTQ is therefore unordered within each pass but not interleaved between passes.
• Delta harnesses that align eidolon FASTQs now get shuffled input by default. germline_e2e's SHUFFLE=1 becomes redundant, and nothing there depends on order.

🤖 Generated with Claude Code

[joshfactorial]

Both review points addressed in 7b4e54c.

File descriptors at the bucket cap.

• The scatter no longer holds a writer per bucket per mate. Each bucket accumulates in a 256 KiB buffer; a full buffer is compressed as one gzip member and appended to the bucket file, which is then closed again. The scatter therefore has at most one bucket file open, at any bucket count.
• Pass 2 opens at most a few files per thread, with at most 8 threads.
• With descriptors no longer the constraint, MAX_BUCKETS rises from 256 to 512. Its limit is now scatter buffer memory: 512 × 2 × 256 KiB = 256 MiB.
• Resource test (cap_path_runs_under_a_low_file_descriptor_limit, Linux): the test binary re-runs itself under bash -c 'ulimit -n 48 …' and shuffles paired input through 512 buckets. The child prints a marker, so a filter that matched nothing cannot pass. On the previous code it fails with Too many open files (os error 24) at bucket 22; it passes now.
• Mutant: skipping the final buffer flush fails 6 unit tests.

Thread-count determinism on raw bytes.

• tests/fastq_shuffle_output.rs now also compares std::fs::read() of the 1- and 4-thread FASTQs.
• The unit test compares raw output bytes at 1 and 8 threads, with 16 buckets.
• Mutant: appending an empty gzip member only when threads > 1 passes the decompressed-record assertions and fails both byte assertions.

Cost after the change (CHM13 chr21, 20×, ~3M pairs, 8 threads):

• wall time 117.8 s, against 112.9 s before the change and 99.4 s with the shuffle off. The per-flush open and close costs ~5 s here;
• peak memory 2.08 GB, against 2.34 GB;
• output size identical to the previous implementation: 372,120,594 bytes.

BlockGzWriter::with_level is dropped as unused. The workspace suite (--no-fail-fast), clippy -D warnings and fmt all pass.