[joshfactorial]

Refs #861 (steps 3 and 4 of its plan). Harness only; no product code changes.

What

• stage_raw_pairs.sbatch accepts local reads. If $WORK/r1.fastq.gz and r2.fastq.gz exist, it skips ENA, does not need SRR, and derives N_PAIRS from r1, checked against r2. Simulated reads then go through the same bwa-mem2 | fixmate | sort | markdup path as the real ones. A resumed real-data staging behaves as before, since its FASTQs were already line-count checked on download.

• telomere_sim.sbatch prepares the matched simulation:

  • fits the fragment model from the real CHM13 BAM, then validates it with validate_frag_model.sh;
  • simulates 2×151 at the real library's depth, rounded because coverage is an integer;
  • draws reads from hs1 minus chrY, because CHM13 is 46,XX and hs1 carries HG002's Y;
  • uses the shipped default quality and error models.

• telomere_ends.sbatch measures every chromosome end, in both libraries:

  • array length in the reference;
  • depth in 100 bp bins over the terminal 20 kb, at all MAPQ and at MAPQ ≥ 20;
  • MAPQ-0 and clip fractions over the array;
  • from the FASTQs, reads with 4× TTAGGG/CCCTAA per million.

  Depth is relative to each library's median depth over interior control windows, placed 75% along each contig. Hard failures:

  • a region with a short depth record;
  • FASTQ read count ≠ BAM primary count;
  • no contigs selected.

Evidence (local, CHM13 chr21 golden BAM at 30x)

• Must-not-fire: the same BAM passed as both libraries gives identical columns.

• Known answers:

  • Arrays: p = 2.7 kb, q = 4.3 kb. The q array has a degenerate stretch at 3.7–3.8 kb (8 of 16 copies), then full repeat to 4.3 kb. The first version stopped at that stretch; the scan now tolerates up to two weak windows in a row.
  • Edge depth 0.52 / 0.62 of interior, consistent with the taper recorded in Telomere realism: compare eidolon against real CHM13 reads at T2T telomeres #861 (0 → 28× over ~400 bp).
  • 1,345 telomeric reads. 7.0 kb of array × 29.8× / 151 bp ≈ 1,380 before the end taper.

• Must-fire mutants:

  edit	result
  drop every other p-array read	array 0.97 → 0.49, flank unchanged
  MAPQ 0 on q-array reads	q20 depth 0.94 → 0.02, MAPQ0 0 → 0.956
  soft-clip p-array reads	clip 0 → 0.885
  FASTQs from a different BAM	FATAL, read counts differ

• Chain: telomere_sim → stage_raw_pairs ran locally on a 2-contig reference:

  • chrY was excluded: all 49,669 pairs came from the other contig;
  • the fragment fit validated against its BAM;
  • depth came out at 5.0×, as requested;
  • the local-reads path derived the pair count with no SRR, and the primary-count check passed.

Not verified

• Not run on Delta or on real CHM13. No CI tests for these three scripts; they are validation harnesses and were checked by hand as above.
• The control windows are not screened for repeats, so a control landing in a satellite would skew normalization. The script prints control MAPQ0 for each library, so that is visible.

🤖 Generated with Claude Code