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

crates/navigator-analysis/examples/archaic_classify_dump.rs

@@ -0,0 +1,38 @@
+//! Dump the archaic **diagnostic** sites (position, derived base, lineage class) as TSV.
+//!
+//! Written to test a different observable for the Tier B HMM. The current model counts *all*
+//! private variants per window, and that signal is weak: measured on a real European, archaic
+//! tracts carry only 2.89x the background density while the background itself varies 5.3x between
+//! its 10th and 90th percentile. Restricting the observable to sites where the derived allele is
+//! actually known to be archaic should be far more specific.
+//!
+//! ```sh
+//! cargo run --release -p navigator-analysis --example archaic_classify_dump -- \
+//!   ~/.decodingus/ancestry/archaic_classify_chm13v2.0.bin chr21 chr22 > classify.tsv
+//! ```
+
+use navigator_analysis::archaic::ArchaicClassify;
+
+fn main() -> Result<(), Box<dyn std::error::Error>> {
+    let mut a = std::env::args().skip(1);
+    let path = a.next().expect("usage: archaic_classify_dump <classify.bin> [contig ...]");
+    let want: Vec<String> = a.collect();
+
+    let cls = ArchaicClassify::from_bytes(&std::fs::read(&path)?).map_err(|e| e.to_string())?;
+    println!("contig\tposition\tderived\tclass");
+    for c in &cls.contigs {
+        let name = &c.positions.contig;
+        if !want.is_empty() && !want.contains(name) {
+            continue;
+        }
+        let mut n = 0usize;
+        for (i, p) in c.positions.iter().enumerate() {
+            let d = c.derived.get(i).copied().unwrap_or(b'N') as char;
+            let k = c.classes.get(i).copied().unwrap_or(2);
+            println!("{name}\t{p}\t{d}\t{k}");
+            n += 1;
+        }
+        eprintln!("{name}: {n} diagnostic sites");
+    }
+    Ok(())
+}

crates/navigator-analysis/examples/archaic_outgroup_density.rs

@@ -0,0 +1,45 @@
+//! Per-window counts of African-outgroup segregating sites — a candidate local mutation-rate proxy.
+//!
+//! The Tier B emission model assumes one background rate genome-wide. Measured, the background
+//! private-variant density varies 5.3x between its 10th and 90th percentile and is 14.6x
+//! overdispersed relative to the Poisson it is modelled with, which is larger than the 2.89x
+//! enrichment inside real archaic tracts — so the model calls its own upper tail archaic. hmmix
+//! avoids this with a mutation-rate map; we have no such asset.
+//!
+//! The density of sites segregating in Africans is already in `archaic_outgroup_af_<build>.bin` and
+//! is a direct measure of how variable a region is, for reasons that have nothing to do with
+//! archaic introgression (mutation rate, reference quality, mappability). This dumps it so that
+//! proxy can be tested as a normalizer before an asset is built for the purpose.
+//!
+//! ```sh
+//! cargo run --release -p navigator-analysis --example archaic_outgroup_density -- \
+//!   ~/.decodingus/ancestry/archaic_outgroup_af_chm13v2.0.bin 1000 chr21 chr22 > og_density.tsv
+//! ```
+
+use navigator_analysis::archaic::ArchaicOutgroup;
+
+fn main() -> Result<(), Box<dyn std::error::Error>> {
+    let mut a = std::env::args().skip(1);
+    let path = a.next().expect("usage: archaic_outgroup_density <outgroup.bin> [window_bp] [contig ...]");
+    let window: i64 = a.next().and_then(|s| s.parse().ok()).unwrap_or(1000);
+    let want: Vec<String> = a.collect();
+
+    let og = ArchaicOutgroup::from_bytes(&std::fs::read(&path)?).map_err(|e| e.to_string())?;
+    println!("contig\twindow_start\tn_outgroup_sites");
+    for c in &og.contigs {
+        if !want.is_empty() && !want.contains(&c.contig) {
+            continue;
+        }
+        let mut counts: std::collections::BTreeMap<i64, u32> = Default::default();
+        let mut n = 0u64;
+        for p in c.iter() {
+            *counts.entry(p / window * window).or_insert(0) += 1;
+            n += 1;
+        }
+        eprintln!("{}: {n} outgroup sites in {} non-empty windows", c.contig, counts.len());
+        for (w, k) in counts {
+            println!("{}\t{}\t{}", c.contig, w, k);
+        }
+    }
+    Ok(())
+}

crates/navigator-analysis/examples/archaic_private_dump.rs

@@ -0,0 +1,50 @@
+//! Dump the **private** variant positions the Tier B HMM actually sees — the subject's derived
+//! variants after the African-outgroup strip — so the input to the model can be checked against an
+//! external truth set independently of the model.
+//!
+//! The segment caller is a density model over exactly these positions. If they are not enriched
+//! inside known archaic tracts, no amount of HMM tuning can help, and the fault is upstream in the
+//! variant calls or the outgroup strip rather than in the model. That question is unanswerable from
+//! the caller's own output, which is why this exists.
+//!
+//! ```sh
+//! cargo run --release -p navigator-analysis --example archaic_private_dump -- \
+//!   calls.json ~/.decodingus/ancestry/archaic_outgroup_af_chm13v2.0.bin > private.tsv
+//! ```
+
+use navigator_analysis::archaic::ArchaicOutgroup;
+use navigator_analysis::caller::SiteGenotype;
+
+fn main() -> Result<(), Box<dyn std::error::Error>> {
+    let mut a = std::env::args().skip(1);
+    let calls_path = a.next().expect("usage: archaic_private_dump <calls.json> <outgroup.bin>");
+    let og_path = a.next().expect("usage: archaic_private_dump <calls.json> <outgroup.bin>");
+
+    let calls: Vec<SiteGenotype> = serde_json::from_str(&std::fs::read_to_string(&calls_path)?)?;
+    let og = ArchaicOutgroup::from_bytes(&std::fs::read(&og_path)?).map_err(|e| e.to_string())?;
+
+    // Group by contig, mirroring what the caller does before it strips.
+    let mut by_contig: std::collections::BTreeMap<String, Vec<&SiteGenotype>> = Default::default();
+    for c in &calls {
+        by_contig.entry(c.contig.clone()).or_default().push(c);
+    }
+
+    // Quality columns come out too: whether the background's excess variance is real biology or
+    // this caller's own error rate varying by region is not answerable without them.
+    println!("contig\tposition\tdosage\tgq\tdepth");
+    for (contig, mut sites) in by_contig {
+        sites.sort_by_key(|s| s.position);
+        let carried: Vec<&SiteGenotype> = sites.iter().copied().filter(|s| s.dosage > 0).collect();
+        let positions: Vec<i64> = carried.iter().map(|s| s.position).collect();
+        let keep: std::collections::HashSet<i64> = og.retain_private(&contig, &positions).into_iter().collect();
+        let mut kept = 0usize;
+        for s in &carried {
+            if keep.contains(&s.position) {
+                println!("{contig}\t{}\t{}\t{}\t{}", s.position, s.dosage, s.gq, s.depth);
+                kept += 1;
+            }
+        }
+        eprintln!("{contig}: {kept} private of {} carried variants", carried.len());
+    }
+    Ok(())
+}

documents/design/ArchaicAncestry_Design.md

@@ -1,11 +1,15 @@
 # Archaic Ancestry Report (Neanderthal / Denisovan) — Design
 
 **Status:** **Tier A SHIPPED** (`v0.1.0-alpha.14`). **Tier B GATED OFF** (`ARCHAIC_SEGMENTS_ENABLED
-= false`) after per-individual validation showed the segment caller carries no per-person signal —
-see *Tier B validation* at the end of §10. It shipped enabled in alpha.14 and was withdrawn in the
-next release. Drafted 2026-07-23; plan added 2026-07-26; all three §9 questions resolved.
-**Read *Tier B validation* and *Deviations from the plan* (end of §10) before trusting anything
-below about Tier B, §7's expected percentage, or M3's feature-gate rule.**
+= false`) after per-individual validation showed the segment caller carries no per-person signal, and
+diagnosed as **built on the wrong observable** — §3's choice of a method designed for people who do
+*not* have archaic reference genomes, which we do. It shipped enabled in alpha.14 and was withdrawn
+in the next release. Drafted 2026-07-23; plan added 2026-07-26; all three §9 questions resolved.
+
+> **Read *Tier B validation* and *Why it failed* (both at the end of §10) before anything else in
+> this document about Tier B.** §3's method choice, §5's Tier B pipeline and M3's calibration are all
+> superseded by that diagnosis. *Deviations from the plan* also qualifies §7's expected percentage
+> and M3's feature-gate rule. Tier A (§5 Tier A, M1, M2) stands unaffected.
 **Goal:** Reconstruct a 23andMe-style Neanderthal report — and go beyond it with a Denisovan
 estimate and a true whole-genome introgression map — from public archaic reference genomes and
 recent methods, using the app's existing ancestry/panel/HMM machinery.
@@ -95,6 +99,15 @@ are population-scale; **IBDmix** is reference-based and heavier; **DAIseg** (joi
 Nea/Den, no post-processing) is attractive but an unreviewed 2025 preprint with self-benchmarks — a
 possible future upgrade, not the v1 foundation.
 
+> **THIS CHOICE WAS WRONG, and the one-clause dismissal of IBDmix is the root of it.** hmmix's
+> premise is detecting introgression *without* archaic reference genomes — it infers them from
+> private-mutation density because it assumes you do not have them. We do: all four archaic genomes,
+> and `archaic_classify` with 2,031,406 diagnostic sites derived from them. Measured, the density
+> observable leaves ~1 informative variant per median tract (14.3 % sensitivity at 5 % false
+> positives) where reference-based allele matching leaves ~30 (95.1 %). "Reference-based and
+> heavier" traded the only thing that made the problem tractable for a saving that was never needed.
+> See *Why it failed* at the end of §10.
+
 ---
 
 ## 3a. Marker-list sourcing — RESOLVED: compute our own
@@ -1038,7 +1051,108 @@ intersection with real 23andMe v5 chip content). It remains enabled and is what
 
 ### Re-enabling
 
-Needs a method change, not a threshold sweep — Skov 2020 matches a segment's whole haplotype against
-each archaic genome relative to a background expectation, where ours tests private-variant density
-against pre-classified sites. Re-running the validation harness that produced the numbers above is
-the gate.
+Needs a method change, not a threshold sweep. **The change required is the OBSERVABLE, not the
+model** — see *Why it failed* immediately below. Re-running the validation harness that produced the
+numbers above is the gate.
+
+---
+
+## Why it failed (2026-07-31) — we chose a method for a problem we do not have
+
+The caller was not mis-tuned. It was built on the wrong observable, and no amount of fitting could
+have rescued it. This section records how that was established, because the wrong turn is easy to
+repeat and every intermediate hypothesis here was plausible.
+
+### The chain of diagnosis
+
+Each step was measured on HG00096, chr21+22, against hmmix's own calls for the same individual.
+
+1. **The signal is real and the pipeline upstream is sound.** Private-variant density inside real
+   archaic tracts is **2.89x** background (357.8 vs 123.9/Mb). The outgroup strip and the variant
+   calls both work.
+2. **The hard-coded `archaic_rate_multiple = 6.0` is roughly double the real effect.** At 1 kb
+   windows and 0.124 background variants/window, the Poisson crossover lands at k* ~ 0.35 — *any*
+   window containing a single private variant is called archaic.
+3. **Per-individual Baum-Welch EM (what hmmix does and we do not) makes it worse.** Unconstrained,
+   it learns a 22.4x multiple and a 0.116 switch probability (~9 kb tracts), calling 17.4 Mb in
+   1,452 segments — 7x the truth. It fits the zero-inflation of sparse counts, not spatial structure.
+4. **Oracle parameters do not help either.** Given the *measured* rates and several tract-length
+   priors, precision stays pinned at 4-5% and every configuration sits at or near its
+   random-placement null. **The model class is wrong, not its parameters.**
+5. **The background is not the distribution the emission assumes.** In callable, non-archaic 100 kb
+   bins: p10-p90 spread **5.3x**, **overdispersion 14.6x** against the Poisson's assumed 1.0x, and
+   **5.9%** of pure-background bins already exceed 2.89x the median. The noise is bigger than the
+   signal, so a single-lambda model calls its own upper tail archaic — and background bins outnumber
+   archaic ones ~20:1.
+6. **A mutation-rate map would not have been enough.** African-outgroup site density — the best
+   proxy available without new data — explains only **38%** of the background variance and takes
+   overdispersion 14.6x -> 7.4x, still above the signal. Quality filtering does not help either: it
+   lowers overdispersion only by discarding variants proportionally, and enrichment falls with it.
+7. **It is not our variant calling.** HG00096 is in the 1000G callset, so the same quantity is
+   computable from their calls. Ours are noisier — 2.12x the carried SNVs on chr21 (93,643 vs
+   44,255), 6x the private-variant density, 97% of our private calls unique to us — but the
+   **contrast is identical**: 1.98x for their calls against 2.08x for ours. Noise dilutes counts,
+   not the ratio.
+8. **It is not the truth set or the lift.** Checked in native hg38 with no lifting, using only
+   hmmix's own segments and their own DAV SNP list: their tracts are enriched **1.84x** for their
+   own archaic SNPs against a null of 1.04x (p95 1.39x). The callset is internally consistent and
+   the lift preserved it.
+
+### The actual cause: one bit of evidence per tract
+
+At ~20-120 private variants/Mb, the **median 36 kb tract contains about ONE informative variant**
+(1.3 expected inside a tract against 0.64 in background). That is the entire evidence separating a
+real tract from noise. Detectability, computed directly:
+
+| observable | evidence per 36 kb tract | sensitivity at <=5% FP | at <=1% FP |
+|---|---|---|---|
+| private-variant density (built) | ~1 variant | **14.3 %** | 4.3 % |
+| archaic-allele matching at diagnostic sites | ~30 sites | **95.1 %** | 80.9 % |
+
+To reach 80% sensitivity at 5% false positives, allele matching needs **36 kb** — the real median
+tract. Density does not get there **at 500 kb**, and hmmix's p10 tract is 7 kb. Both observables
+carry the *same* ~3x contrast (2.89x vs 3.04x), so contrast was never the problem: **evidence per
+tract** was.
+
+### The wrong turn, named
+
+§3 chose **Skov 2018 (hmmix)**, whose premise is detecting introgression **without archaic reference
+genomes** — it infers them indirectly from private-mutation density precisely because it assumes you
+do not have them. We adopted it while **holding all four archaic genomes** and already shipping
+`archaic_classify` with 2,031,406 diagnostic sites derived from them. §3 rejected the reference-based
+alternative in one clause — *"IBDmix is reference-based and heavier"* — and that clause is the root
+of every failure above.
+
+An introgressed tract is **a haplotype inherited intact from an archaic ancestor**. The question is
+whether a stretch of the genome *matches an archaic genome*, not whether it is slightly more mutated
+than average. Runs, not counts.
+
+Measured, the reference-based observable gives **39.5%** archaic-allele carrying inside real tracts
+against **13.0%** elsewhere, over ~30 sites per tract instead of ~1 variant.
+
+### What the reframe also fixes
+
+- **The mutation-rate map becomes unnecessary.** Diagnostic sites are the denominator, so their
+  uneven density cancels out — the thing the map existed to correct.
+- **The 5.3x background spread and 14.6x overdispersion stop applying**, because nothing is being
+  modelled as a density any more.
+- **The assets already exist.** `archaic_classify_<build>.bin` carries the sites, derived bases and
+  per-lineage classes; no new offline pipeline is needed to start.
+- **Attribution may come back with it.** The same per-site matching, split by lineage class, *is*
+  Skov's post-hoc annotation — the thing `attribute_lineage` was gated for.
+
+### Open before this becomes a plan
+
+- **The background carrying rate measures 13.0%**, against the 4.3% recorded earlier in this
+  document. Unreconciled; probably our call noise. It sets the contrast, so it matters.
+- **~30 sites per tract assumes adequate coverage at diagnostic sites.** Check against real call
+  rates rather than assuming.
+- **n = 1 for the reframe.** The 3.04x enrichment is HG00096 only. The harness in
+  `scripts/archaic-validation/` runs the cohort cheaply now that the CRAM fix has landed.
+
+### Instruments
+
+The examples these numbers came from, all under `crates/navigator-analysis/examples/`:
+`archaic_private_dump` (the HMM's actual input, with quality columns), `archaic_outgroup_density`
+(the rate-map proxy), `archaic_classify_dump` (diagnostic sites), `archaic_callable_dump` (what the
+caller can see at all), and `cram_query_probe` (the CRAM defect found on the way here).

scripts/archaic-validation/observable/README.md

@@ -0,0 +1,52 @@
+# Why Tier B failed — the diagnostic scripts
+
+These are the instruments behind
+[`ArchaicAncestry_Design.md`](../../../documents/design/ArchaicAncestry_Design.md) § *Why it failed*.
+They establish that the Tier B segment caller was built on the **wrong observable**, not mis-tuned —
+and, just as importantly, they rule out the explanations that looked right first.
+
+Kept because every one of the rejected hypotheses is plausible enough to be proposed again.
+
+## Run order
+
+All expect a scratch directory containing, for one individual (HG00096 in the recorded run):
+
+| file | produced by |
+|---|---|
+| `HG00096.chr{21,22}.calls.json` | `sqlite3` on `analysis_artifact`, kind `diploid_denovo:chrN` |
+| `private.chr{21,22}.tsv` | `cargo run --example archaic_private_dump` |
+| `classify.tsv` | `cargo run --example archaic_classify_dump` |
+| `callable.bed` | `cargo run --example archaic_callable_dump` (threshold 0.5) |
+| `og_density.tsv` | `cargo run --example archaic_outgroup_density` |
+| `truth_<sample>.chm13.bed` | hmmix segments lifted hg38→CHM13 — see [../README.md](../README.md) |
+
+## What each one answers
+
+| script | question | answer it gave |
+|---|---|---|
+| `background_variation.py` | is the background the flat Poisson the emission assumes? | **No** — 5.3× p10–p90 spread, **14.6× overdispersed**, larger than the 2.89× signal |
+| `mutrate_check.py` | would a mutation-rate map fix that? | **Not enough** — the best available proxy explains 38% of the variance, leaving 7.4× |
+| `quality_effect.py` | is the excess variance our caller's artifacts? | **No** — filtering lowers overdispersion only by discarding variants; enrichment falls with it |
+| `callset_compare.py` | is our variant calling diluting the signal? | **No** — 1000G's own calls for the same person give the *same* contrast (1.98× vs 2.08×), though ours are 6× noisier |
+| `truth_selfcheck.py` | is the truth set (or my lift) wrong? | **No** — hmmix's tracts are enriched 1.84× for their own archaic SNPs in native hg38, null 1.04× |
+| `observable_compare.py` | is a more specific observable available? | private ∩ diagnostic is unusable (7 sites); the African strip removes almost all diagnostic sites |
+| `haplotype_match.py` | does archaic-allele matching separate tracts? | **Yes** — 39.5% carrying inside tracts vs 13.0% elsewhere, over ~30 sites per tract |
+| `discriminability.py` | how detectable is ONE tract under each observable? | density **14.3%** vs matching **95.1%** sensitivity at 5% false positives |
+
+`discriminability.py` is the one to read first. It needs no inputs and is the whole argument:
+both observables carry ~3× contrast, so contrast was never the problem — **evidence per tract** was.
+
+## The trap these scripts exist to prevent
+
+Two harness bugs in here produced confident wrong answers before being caught, both worth knowing:
+
+- **`haplotype_match.py` originally conditioned on the subject having a variant call** at a
+  diagnostic site, which samples only sites where he already has a variant and reported an
+  impossible ~80% carrying rate against a known 4.3% background. Every diagnostic site in callable
+  territory belongs in the denominator; no call means hom-reference.
+- **Binning at 100 kb to measure in-tract contrast** dilutes a 36 kb tract with 64 kb of background
+  and reported 1.14× where the correct measurement is 1.98×. Measure contrast at tract boundaries,
+  not in fixed bins.
+
+More generally: sensitivity alone is gameable by calling more sequence, because the
+random-placement null rises with it. Always report the null at the extent actually called.