Make subsampling contamination-aware

[jbrestel]

The problem

Subsampling runs before any alignment, and calculateMaxReads() was a pure function of
assayType and genomeSize. So the invariant enforced was "N raw reads" while what
downstream workflows actually need is "N on-target reads."

Those coincide only for a pure sample. With substantial host contamination in a mixed sample,
cutting to N raw reads delivers far fewer than N on-target reads, and every child workflow
gets less coverage than it asked for. The pipeline was asserting a coverage guarantee it had
no information to make, because it has never seen an alignment.

The approach

Measure the on-target fraction here, so the cut is made against a number that means something.

• SKETCH_REFERENCE sketches a user-supplied target-organism FASTA once per run and
  measures genome size.
• MEASURE_SAMPLE draws a pilot of reads per sample, sketches it, and runs sourmash gather
  to get f_unique_weighted — the on-target fraction.
• depthPlan (a pure, unit-tested function) inflates the retained read count:
  rawReads = targetOnTargetReads / onTargetFraction, capped at reads available.

A sample that is 12% target retains ~8.3x more raw reads than a clean one, and both land at
the requested coverage after mapping. Verified arithmetic: a 100M-read run at 12% target
against a 23Mb genome at 60x retains 76,666,667 reads and yields exactly 60.0x on-target.

The duplicated CONCATENATE → SUBSAMPLE → FORMAT → collectFile tail that existed in both
main.nf and retrieve_from_sra.nf is extracted into a shared PREPARE_SAMPLES subworkflow.
The fromSra fork now decides only how reads are obtained.

Interface changes

Added	--referenceFasta (required), --targetCoverage, --minOnTargetFraction, --minPlausibleFraction, --pilotSize
Removed	--genomeSize — now measured from the reference FASTA
New output	sample_metrics.csv sidecar

samplesheet.csv's contract is unchanged at sample,fastq_1,fastq_2,var1 — verified
byte-for-byte. Metrics go to a separate sidecar specifically so downstream consumers are
unaffected.

Behavior changes to be aware of

1. --referenceFasta is required. Runs fail immediately without it, before any download.
2. An SRA run where fasterq-dump emits 3 files now fails instead of succeeding. Default
   --split-3 emits R1, R2 and an orphan file when some reads lack a mate; the vendored
   module's own .diff renames that orphan so it reached CONCATENATE_FASTQ, which treated
   3 files as single-end and concatenated R1+R2+orphans into one blob — silently destroying
   mate pairing. This predates this branch. Anyone who has processed such a run will now see
   a loud failure where they previously got quietly corrupted output.
3. An unrecognized --assayType now fails rather than silently defaulting to DNASeq.

Bugs fixed along the way

• Protein FASTA passed the "is this a nucleotide FASTA" gate. A/C/G/T/N are all valid
  amino-acid codes, so a proteome cleared the old absolute base-count floor and produced a
  fabricated genome size. Now gated on strict-ACGTN as a ratio (>=90%). IUPAC ambiguity
  codes are deliberately excluded from that ratio — they are also amino-acid codes, and
  including them lets protein score ~70%.
• sourmash gather needs --threshold-bp 0. The default 50kbp threshold writes no result
  row at low overlap, so a 0.2%-target sample parsed as 0.0 — indistinguishable from a wrong
  reference genome. Measured: 0.2% reads as nothing by default, 0.0018 with the flag.
• || true on gather swallowed real crashes. An OOM or corrupt sketch reported 0.0
  on-target, which the policy reads as maximal contamination and inflates to the cap. Verified
  gather exits 0 on genuine no-match, so the guard was unnecessary and is gone.
• Zero-read samples crashed with a bare Division by zero in a shared function. Now fails
  naming the sample.
• -stub-run was completely broken — concatenate_fastq.nf's stub referenced an
  undefined bare hasPairedReads. Never caught because nothing here had ever run a stub.
• Same-basename inputs collided. A multi-lane samplesheet (lane1/s.fastq.gz,
  lane2/s.fastq.gz) failed with an input-collision error in the very module meant to
  concatenate them. Fixed with stageAs: "?/*".
• Single-run samples were decompressed and re-gzipped to produce byte-identical output.
  Now symlinked — concat is the pipeline's peak-disk step.

Test plan

nf-test was not installed and had never run here (nextflow test in CLAUDE.md is not a real
command). Bootstrapped, plus 36 tests across 5 files:

• depth_policy (21) — clamps, floors, RNASeq path, flagging, and loud failure on bad input
• sketch_reference (5) — genome measurement, protein rejection, headers-only, gzipped
• measure_sample (3) — known 10% mixture (0.115 vs 0.10 mixing ratio), zero-overlap, reservoir path
• concatenate_fastq (4) — byte-identity passthrough, multi-file concat, mate-collision fall-through
• subsample_fastq (3) — cut, passthrough, paired-end pairing integrity
• Two-sample end-to-end run; published FASTQs verified readable through the symlink chain
• Stub runs, single-end and paired
• Missing-reference error fires before any process submission
• All-flagged run triggers the wrong-genome diagnostic
• Protein FASTA as reference is rejected

Known gaps

• No end-to-end test demonstrates inflation. The 1,000,000-read MIN_TARGET_READS floor
  means any fixture-scale run has rawReads capped at totalReads, so inflation never
  changes the outcome. The arithmetic is unit-tested and the wiring is integration-tested;
  nothing is untested, but no single test exercises the core behavior. Closing it would need
  either gigabyte fixtures or a configurable floor.
• Paired-end coverage runs ~2x high. total_reads counts pairs while coverage credits one
  read length per pair. Pre-existing, deliberately deferred, now documented in README.md.
• README.org is stale — still documents --genomeSize/--maxReads. It predates
  README.md and appears to be an abandoned duplicate. Left untouched pending a decision to
  delete or maintain it.

🤖 Generated with Claude Code