subsample: Add per-tile binning with skip-ahead sampling and threading

[jdidion]

Summary

Adds per-tile binning to the subsample command for normalizing FASTQ files based on physical sequencing location, mitigating batch effects from uneven coverage or flow cell edge effects.

• Reads are binned by lane and tile parsed from Illumina read headers (@<instrument>:<run>:<flowcell>:<lane>:<tile>:<x>:<y>)
• Bins with fewer than the target number of records are discarded
• Records are randomly sampled from each retained bin
• Supports single-end and paired-end reads

New CLI options

Option	Description
--record-count-per-tile N	Explicit per-tile count (exclusive with -n/-p)
--bin-by-tile	Enable tile binning with -n or -p (auto-computes per-tile count)
--fast	Use faster skip-ahead sampling instead of default exact method
--in-memory	Keep tiles in RAM instead of writing temp files
--temp-dir DIR	Directory for temp tile files (default: system temp)
--sampling-threads N	Parallel tile processing (default: 1)
--compression-threads N	Parallel gzip output (default: 1)

Three tile count modes

• --record-count-per-tile N — explicit per-tile count
• -n N --bin-by-tile — per-tile = N / num_bins
• -p P --bin-by-tile — per-tile = floor(P * total / num_bins)

Three storage/sampling strategies

Mode	Storage	Sampling	Exact counts?
Default	temp files + offset index	indexed seek	yes
--fast	temp files	skip-ahead (exponential byte jumps)	approximate
--in-memory	RAM	random index selection	yes

Usage examples

# Exact 10,000 reads per tile (paper's method)
fq subsample --record-count-per-tile 10000 --seed 42 \
  --r1-dst out.R1.fq.gz --r2-dst out.R2.fq.gz in.R1.fq.gz in.R2.fq.gz

# Auto-compute per-tile from target count, parallel sampling
fq subsample -n 100000 --bin-by-tile --sampling-threads 4 --seed 42 \
  --r1-dst out.R1.fq.gz in.R1.fq.gz

# Fast approximate mode with custom temp dir
fq subsample --record-count-per-tile 10000 --fast --temp-dir /scratch --seed 42 \
  --r1-dst out.R1.fq.gz in.R1.fq.gz

# In-memory mode (no temp files, higher RAM usage)
fq subsample --record-count-per-tile 10000 --in-memory --seed 42 \
  --r1-dst out.R1.fq.gz in.R1.fq.gz

Dependencies

• Added tempfile crate for temp directory management

Test plan

• Illumina header parsing (standard, with description, edge cases)
• Exponential distribution sampling
• FASTQ record boundary detection after seeking
• Record offset index building
• Exact tile sampling (single-end and paired-end)
• Skip-ahead tile sampling (single-end)
• In-memory tile sampling (paired-end)
• Auto per-tile count from --record-count
• Empty bins edge case
• All 73 tests pass (63 original + 10 new)
• Test with real Illumina FASTQ data

🤖 Generated with Claude Code

[zaeleus]

Hi,

I suspect this feature may be a bit too niche to include in fq.

I am curious, however, about the motivation behind per-tile binning. With a sufficiently large number of reads, is there a measurable difference in downstream quality compared to a uniform sampling?

[jdidion]

Hi @zaeleus - this paper was the original motivation: https://www.biorxiv.org/content/10.64898/2026.03.08.710357v1. They describe the per-tile binning method, but do not publish any tool to replicate it. It seemed a straight-forward and logical thing to add to a tool like fq that already supports multiple different downsampling methods.

[zaeleus]

fq only has one downsampling method, i.e., uniform random sampling. Choosing between approximate and exact affects the number of records that are emitted.

I don't think an experimental and untested methodology is a good fit here. Have you tested the practicality of the output, and is there an meaningful improvement over random sampling?

[jdidion]

Yes, I've been using the per-tile subsampled output (using the build from this branch) in an active R&D project.

I'll note that it's not straight-forward to select the value for the per-tile target number of reads, so if this PR is accepted, it would be good to have a companion tool that counts the number of reads per tile and outputs the distribution.

It also seems that, at least for NovaSeq X, BCL collates the read by tile (i.e. all the reads from each tile are consecutive in the fastq). So an indexing procedure that notes the offset of the first read in each tile during the per-tile counting pass would enable parallelization of the per-tile downsampling.