Benchmark --mode fast (and sensitive) against a truth set, not against each other

[wdecoster]

--mode fast (#PR/branch perf/cigar-fast-genotyping) has so far only been compared against the sensitive alignment path on one 30x ONT sample. That measures agreement, not accuracy: it silently assumes the alignment path is right.

Numbers as they stand, 2000 catalog loci on chr1:20-21.4 Mb, single thread:

	sensitive	fast
CPU	1016 s	13 s
peak RSS	643 MB	289 MB

Per-haplotype median read length (MRL) agrees within 2 bp at 95.7% of alleles (median difference 1 bp); 14 of 1760 loci differ by more than 10 bp, nearly all in one region where a ~900 bp insertion sits between four overlapping catalog entries and the two modes attribute it to different ones.

What is needed:

• a truth set — GIAB (HG002) with the tandem repeat benchmark / TR truth VCF, or assembly-derived allele lengths from the HG002 trio assemblies
• run both modes and report per-locus allele-length concordance against truth, stratified by allele length, motif size and locus complexity
• specifically check the case fast is expected to lose on: expansions long enough that reads clip rather than span, where fast drops the read and falls back to sensitive. On the sample tested so far the fallback never triggered, so it is untested against real data.

This also gives the first accuracy number for the sensitive path itself, which has never been benchmarked against a truth set either.

[wdecoster]

The harness is merged as misc/benchmark_truth.py (#28). It runs both modes over a truth set's loci and scores each against the truth allele lengths, stratified by how far the truth allele departs from the reference.

Still open, because the actual benchmark has not been run: nothing GIAB is reachable from the dev machine, so this needs running where HG002 lives.

python misc/benchmark_truth.py \
    --binary ./target/release/STRdust \
    --fasta GRCh38.fa --bam HG002.cram \
    --truth-vcf HG002_TR_benchmark.vcf.gz \
    --confident-bed HG002_TR_benchmark_regions.bed \
    --threads 8 --out-prefix hg002_tr --plot

Worth doing --regions chr1 --max-loci 2000 first to check the truth VCF parses as expected — the script reports how many loci it skipped and why.

Two things the results should settle:

1. Where sensitive and fast disagree, which is closer to the truth. The current numbers (identical MRL at 56% of alleles, within 2 bp at 93%) say nothing about that.
2. Whether the read-dropping rules in fast are set sensibly. The clip rejection and the --support fallback never fired on the sample tested, so the path that is supposed to protect large expansions is effectively untested against real expanded alleles — which is exactly what the >200bp stratum in the summary is there to expose.

[wdecoster]

Benchmarked against the GIAB HG002 tandem-repeat truth set over the adotto v1.2.1 catalog: 10,000 randomly sampled loci, two independent seeds, both modes, ~30x ONT.

fast is the more accurate path: 73.3% exact allele length against sensitive's 53.0% (72.9% vs 52.6% on the second seed).

Contrary to the expectation in this issue, it is better at long expansions rather than worse — 75.0% vs 45.5% exact in the >200bp stratum, with fewer no-calls. Expansion detection is identical in both modes (98% of loci with an expanded allele ≥50bp, 94% at ≥200bp) and neither mode ever reported a real expansion as a small number. The cost is recall, 55% vs 65%, bought back as precision, 61% vs 37%.

Most of the difference traces to the POA consensus in the sensitive path running long: small insertions near the repeat/flank junction get folded into the allele. fast avoids it by reading length off the existing alignment.

Closing — the comparison this issue asked for now exists. The default will probably move to fast; that and the remaining sizing work will be evaluated on actual pathogenic repeat loci before anything changes.