Archaic Tier B: diagnose the failure as the wrong observable, and bank the instruments

[JamesKane]

Tier B was gated off in #40 for carrying no per-individual signal. This records why — the caller wasn't mis-tuned, it was built on an observable that leaves about one bit of evidence per tract, and no amount of fitting could have rescued it.

Documentation and dev-only examples. No shipped behaviour changes; Tier B stays gated.

The number that settles it

Detectability of one 36 kb tract (the measured median):

observable	evidence per tract	sensitivity @ ≤5% FP	@ ≤1% FP
private-variant density (what we built)	~1 variant	14.3%	4.3%
archaic-allele matching at diagnostic sites	~30 sites	95.1%	80.9%

Reaching 80% sensitivity at 5% false positives takes 36 kb by allele matching, and is not reached by 500 kb with density. hmmix's p10 tract is 7 kb.

Both observables carry the same ~3× contrast (2.89× vs 3.04×) — so contrast was never the problem. Evidence per tract was.

The wrong turn

§3 chose Skov 2018 (hmmix), whose premise is detecting introgression without archaic reference genomes — it infers them from private-mutation density precisely because it assumes you don't have them. We do: all four archaic genomes, and archaic_classify already ships 2,031,406 diagnostic sites derived from them.

§3 dismissed the reference-based alternative in one clause — "IBDmix is reference-based and heavier" — and that clause is the root of every failure below. An introgressed tract is a haplotype inherited intact; the question is whether a stretch matches an archaic genome, not whether it's slightly more mutated. Runs, not counts.

What was ruled out first — each looked right at the time

• The hard-coded 6.0× multiple is wrong (real enrichment 2.89×), but correcting it changes little
• Per-individual Baum-Welch EM — the thing hmmix does and we don't — makes it worse: unconstrained it learns 22.4× and ~9 kb tracts, calling 7× the truth
• Oracle parameters taken from the truth still sit at their random-placement null, precision pinned at 4–5% → the model class is wrong, not its parameters
• The background is 14.6× overdispersed against the Poisson it's modelled with (5.3× p10–p90 spread) — noise larger than signal
• A mutation-rate map wouldn't have sufficed: the best proxy explains 38% of that variance, leaving 7.4×
• Quality filtering doesn't help — it lowers dispersion only by discarding variants proportionally
• Not our variant calling: 1000G's own calls for the same person give the same contrast (1.98× vs 2.08×), though ours are 6× noisier
• Not the truth set or my lift: in native hg38 with no lifting, hmmix's tracts are enriched 1.84× for their own archaic SNPs against a 1.04× null

What the reframe also fixes

• The mutation-rate map becomes unnecessary — diagnostic sites are the denominator, so their uneven density cancels
• Background dispersion stops applying — nothing is modelled as a density any more
• The assets already exist — no new offline pipeline to start
• Attribution may come back with it — per-site matching split by lineage class is Skov's post-hoc annotation, the thing attribute_lineage was gated for

Still open

• Background carrying rate measures 13.0% here against the 4.3% recorded earlier in the design. Unreconciled, and it sets the contrast.
• ~30 sites per tract assumes adequate coverage at diagnostic sites — needs checking against real call rates.
• The reframe is n=1 so far. The cohort harness runs cheaply now that CRAM region queries decoded the whole chromosome; skip non-overlapping containers #39 landed.
• Separately: our caller emits 2.12× the SNVs 1000G does on the same genome, with 97% of private calls unique to us. Worth its own investigation, well beyond archaic work.

Banked

Three diagnostic examples (archaic_private_dump, archaic_outgroup_density, archaic_classify_dump) and eight analysis scripts under scripts/archaic-validation/observable/, with a README mapping each to the question it answered.

Two harness bugs that produced confident wrong answers are documented there too, since both are easy to repeat: conditioning on "has a call" (reported ~80% carrying against a known 4.3% background), and measuring tract contrast in 100 kb bins (reported 1.14× where the correct figure is 1.98×).

discriminability.py needs no inputs and is the whole argument in twelve lines of output.

🤖 Generated with Claude Code