Shuffle gen-reads FASTQ output by default (shuffle_fastq, #870)

docs-site/src/gen-reads/fastq-shuffling.md

@@ -1,29 +1,48 @@
 # FASTQ Shuffling
-`eidolon` writes reads in genome order: contig by contig, by position within each contig. A
-real sequencer's FASTQ is in effectively random order. Variant calling does not depend on this,
-since the BAM is coordinate-sorted before calling. Tools that read the FASTQ in batches,
-aligners in particular, do see the difference: a batch can be filled entirely with reads from
-one region, such as a long satellite array. Measured on a 21x simulation of T2T-CHM13 (#870),
-bwa-mem2 peaked at 166 GB on the reads as written and at 90 GB on the same reads shuffled,
-with the same alignments either way. Shuffle before aligning a large genome.
-#870 tracks shuffling in `eidolon` itself.
-
-With `seqkit`, shuffle each mate separately **with the same seed**. seqkit has no paired mode,
-but its permutation depends only on the seed and the number of records, and R1 and R2 have the
-same count. So the two files come out in the same order and mates stay in register:
+`eidolon` writes its FASTQs in **random order** by default (`shuffle_fastq: true`), as a
+sequencer does. Mates stay in register: line *n* of R1 and line *n* of R2 are the two reads of
+one pair.
 
-```bash
-# single-ended
-seqkit shuffle -s 11 sample.fastq.gz -o sample_shuffled.fastq.gz
+Reads are generated chunk by chunk along the genome, so without shuffling the FASTQ comes out in
+genome order: contig by contig, by position within each contig. Variant calling does not depend
+on this, because the BAM is coordinate-sorted before calling. Tools that read the FASTQ in
+batches do see the difference, aligners in particular: a batch can be filled entirely with reads
+from one region, such as a long satellite array. Measured on a 21x simulation of T2T-CHM13
+(#870), bwa-mem2 peaked at 166 GB on reads in genome order and at 90 GB on the same reads
+shuffled, with the same alignments either way.
+
+```yaml
+shuffle_fastq: true    # default; false writes reads in genome order
+```
+
+**What it does and does not change.**
+- Only the order of FASTQ records changes. The reads themselves, the truth VCF and the golden BAM
+  are identical with the shuffle on or off. The shuffle draws from its own branch of the random
+  number generator, so it does not disturb the streams that generate reads.
+- Output is deterministic: the same seed gives byte-identical FASTQs at any `num_threads`.
+- Turning the shuffle off restores genome order exactly as earlier versions wrote it.
+
+**Cost.** The shuffle is a two-pass external sort. Records are scattered at random into bucket
+files in the run's temp directory, then each bucket is shuffled in memory and written out. Buckets
+are shuffled in parallel only as far as they fit in about 4 GB, so memory stays bounded at any
+genome size. Temp disk stays close to one copy of the output, because each input is deleted as it
+is read.
 
-# paired-ended: same -s on both mates
+Measured on T2T-CHM13 chr21 at 20x (about 3M pairs, 8 threads), against `shuffle_fastq: false`:
+- wall time 112 s against 99 s;
+- peak memory 2.3 GB against 0.5 GB;
+- **gzipped FASTQs about 60% larger.** Reads in genome order overlap their neighbors, which gzip
+  exploits; shuffled reads, like a real sequencer's, compress less well.
+
+**Shuffling output from older versions, or from `shuffle_fastq: false`.** With `seqkit`, shuffle
+each mate separately **with the same seed**. seqkit has no paired mode, but its permutation
+depends only on the seed and the number of records, and R1 and R2 have the same count:
+
+```bash
 seqkit shuffle -s 11 sample_r1.fastq.gz -o sample_r1_shuffled.fastq.gz
 seqkit shuffle -s 11 sample_r2.fastq.gz -o sample_r2_shuffled.fastq.gz
 ```
 
-**Memory:** `seqkit shuffle` reads every record into memory, so expect memory on the order of
-the uncompressed FASTQ, per mate. Its low-memory `-2/--two-pass` mode supports FASTA only.
-Do not pass both mates to one `seqkit shuffle` call: they would be read as one input and shuffled
-together into a single output.
-
-(`seqkit` is an open-source FASTQ toolkit: https://github.com/shenwei356/seqkit.)
+`seqkit shuffle` reads every record into memory; its low-memory `-2/--two-pass` mode supports
+FASTA only. Do not pass both mates to one call: they would be shuffled together into a single
+output. (`seqkit`: https://github.com/shenwei356/seqkit.)

docs-site/src/gen-reads/parallel-processing.md

@@ -1,7 +1,7 @@
 # Parallel Processing
 `eidolon gen-reads` processes contigs in parallel by default using rayon's work-stealing thread pool. Each contig is an independent unit of work — variant generation, fragment sampling, and FASTQ/BAM writing all happen concurrently across contigs — so references with many contigs scale well across cores.
 
-Output is **byte-identical regardless of `num_threads`** (the same seed always produces the same reads in the same order), so you can change the thread count freely without affecting results.
+Output is **byte-identical regardless of `num_threads`** (the same seed always produces the same reads in the same order, shuffled or not; see [FASTQ Shuffling](fastq-shuffling.md)), so you can change the thread count freely without affecting results.
 
 A note on scaling: read generation is largely **memory-bandwidth bound**, so on references dominated by one or a few large chromosomes the wall-time gain from extra cores is limited — the bottleneck is memory throughput, not CPU. An experimental sub-contig **chunk size** knob (below) can split a single chromosome across cores, but it is **disabled by default** because in our benchmarks it did not improve wall time and occasionally regressed it.
 

eidolon-core/src/file_tools.rs

@@ -6,6 +6,7 @@ pub mod bam_writer;
 pub mod bed_reader;
 pub mod block_gz;
 pub mod fasta_stream;
+pub mod fastq_shuffle;
 pub mod fastq_tools;
 pub mod file_io;
 pub mod folder_tools;