From 2152cb02dffb05d1cda034b609717828f7eb82ea Mon Sep 17 00:00:00 2001 From: James Kane Date: Fri, 31 Jul 2026 06:48:39 -0500 Subject: [PATCH 1/9] Archaic Tier B: diagnose the failure as the wrong OBSERVABLE, and bank the instruments MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Tier B was gated off (#40) for carrying no per-individual signal. This records WHY, because the caller was not 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. THE NUMBER THAT SETTLES IT. Detectability of one 36 kb tract (the measured median): evidence sens @5% FP @1% FP private-variant density (built) ~1 variant 14.3% 4.3% archaic-allele matching ~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. 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. §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 do not 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 is slightly more mutated. Runs, not counts. WHAT WAS RULED OUT FIRST, each of which looked right at the time: - the 6.0x hard-coded multiple IS wrong (real enrichment 2.89x), but correcting it changes little; - per-individual Baum-Welch EM -- the thing hmmix does and we do not -- makes it WORSE: unconstrained it learns 22.4x and ~9 kb tracts, calling 7x the truth; - oracle parameters taken from the truth still sit at their random-placement null, precision pinned at 4-5%, so the model CLASS is wrong, not its parameters; - the background is 14.6x overdispersed against the Poisson it is modelled with, with a 5.3x p10-p90 spread -- noise larger than the signal; - a mutation-rate map would not have sufficed: the best proxy available explains 38% of that variance, leaving 7.4x; - quality filtering does not help (it lowers dispersion only by discarding variants); - it is NOT our variant calling: 1000G's own calls for the same person give the same contrast (1.98x vs 2.08x), though ours are 6x noisier -- 2.12x the carried SNVs on chr21 and 97% of private calls unique to us, which is worth its own investigation; - it is NOT the truth set or my lift: in native hg38 with no lifting, hmmix's tracts are enriched 1.84x for their own archaic SNPs against a 1.04x null. WHAT THE REFRAME ALSO FIXES: diagnostic sites become the denominator, so their uneven density cancels and the mutation-rate map is unnecessary; the background's dispersion stops applying because nothing is modelled as a density; the assets already exist; and per-site matching split by lineage class IS Skov's post-hoc annotation, so the gated attribute_lineage may be recoverable with it. STILL OPEN: the 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. And the reframe is n=1 so far. Banked: three diagnostic examples (archaic_private_dump, archaic_outgroup_density, archaic_classify_dump) and the eight analysis scripts, with a README mapping each to the question it answered. Two harness bugs that produced confident wrong answers are documented there too -- conditioning on "has a call" (reported ~80% carrying against a known 4.3% background) and measuring tract contrast in 100 kb bins (reported 1.14x where the correct figure is 1.98x). Documentation and dev-only examples; no shipped behaviour changes. Tier B stays gated. Co-Authored-By: Claude Opus 5 (1M context) --- .../examples/archaic_classify_dump.rs | 38 +++++ .../examples/archaic_outgroup_density.rs | 45 ++++++ .../examples/archaic_private_dump.rs | 50 +++++++ documents/design/ArchaicAncestry_Design.md | 132 ++++++++++++++++-- .../archaic-validation/observable/README.md | 52 +++++++ .../observable/background_variation.py | 106 ++++++++++++++ .../observable/callset_compare.py | 100 +++++++++++++ .../observable/discriminability.py | 73 ++++++++++ .../observable/haplotype_match.py | 121 ++++++++++++++++ .../observable/mutrate_check.py | 132 ++++++++++++++++++ .../observable/observable_compare.py | 123 ++++++++++++++++ .../observable/quality_effect.py | 123 ++++++++++++++++ .../observable/truth_selfcheck.py | 120 ++++++++++++++++ 13 files changed, 1206 insertions(+), 9 deletions(-) create mode 100644 crates/navigator-analysis/examples/archaic_classify_dump.rs create mode 100644 crates/navigator-analysis/examples/archaic_outgroup_density.rs create mode 100644 crates/navigator-analysis/examples/archaic_private_dump.rs create mode 100644 scripts/archaic-validation/observable/README.md create mode 100644 scripts/archaic-validation/observable/background_variation.py create mode 100644 scripts/archaic-validation/observable/callset_compare.py create mode 100644 scripts/archaic-validation/observable/discriminability.py create mode 100644 scripts/archaic-validation/observable/haplotype_match.py create mode 100644 scripts/archaic-validation/observable/mutrate_check.py create mode 100644 scripts/archaic-validation/observable/observable_compare.py create mode 100644 scripts/archaic-validation/observable/quality_effect.py create mode 100644 scripts/archaic-validation/observable/truth_selfcheck.py diff --git a/crates/navigator-analysis/examples/archaic_classify_dump.rs b/crates/navigator-analysis/examples/archaic_classify_dump.rs new file mode 100644 index 0000000..2150622 --- /dev/null +++ b/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> { + let mut a = std::env::args().skip(1); + let path = a.next().expect("usage: archaic_classify_dump [contig ...]"); + let want: Vec = 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(()) +} diff --git a/crates/navigator-analysis/examples/archaic_outgroup_density.rs b/crates/navigator-analysis/examples/archaic_outgroup_density.rs new file mode 100644 index 0000000..78e356e --- /dev/null +++ b/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_.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> { + let mut a = std::env::args().skip(1); + let path = a.next().expect("usage: archaic_outgroup_density [window_bp] [contig ...]"); + let window: i64 = a.next().and_then(|s| s.parse().ok()).unwrap_or(1000); + let want: Vec = 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 = 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(()) +} diff --git a/crates/navigator-analysis/examples/archaic_private_dump.rs b/crates/navigator-analysis/examples/archaic_private_dump.rs new file mode 100644 index 0000000..86c573b --- /dev/null +++ b/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> { + let mut a = std::env::args().skip(1); + let calls_path = a.next().expect("usage: archaic_private_dump "); + let og_path = a.next().expect("usage: archaic_private_dump "); + + let calls: Vec = 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> = 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 = carried.iter().map(|s| s.position).collect(); + let keep: std::collections::HashSet = 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(()) +} diff --git a/documents/design/ArchaicAncestry_Design.md b/documents/design/ArchaicAncestry_Design.md index 3174393..851c248 100644 --- a/documents/design/ArchaicAncestry_Design.md +++ b/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_.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). diff --git a/scripts/archaic-validation/observable/README.md b/scripts/archaic-validation/observable/README.md new file mode 100644 index 0000000..c1e635d --- /dev/null +++ b/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_.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. diff --git a/scripts/archaic-validation/observable/background_variation.py b/scripts/archaic-validation/observable/background_variation.py new file mode 100644 index 0000000..214d26c --- /dev/null +++ b/scripts/archaic-validation/observable/background_variation.py @@ -0,0 +1,106 @@ +"""Is the private-variant background uniform, as the caller's emission model assumes? + +The model says: background windows emit Poisson(lambda) with ONE genome-wide lambda, so a window +with several private variants is evidence of an archaic tract. That inference only holds if the +non-archaic background is actually flat. hmmix does not assume this — it requires a mutation-rate +map and scales the expected density per window by it. We have no such asset. + +This measures the spread of background density directly, in regions that are callable and NOT in +hmmix's archaic tracts, and asks how a Poisson model would fare against it. +""" + +import bisect +import collections +import math + +BIN = 100_000 +CONTIGS = ('chr21', 'chr22') + + +def merge(iv, tol=0): + iv = sorted(iv) + out, cs, ce = [], *iv[0] + for s, e in iv[1:]: + if s > ce + tol: + out.append((cs, ce)) + cs, ce = s, e + else: + ce = max(ce, e) + out.append((cs, ce)) + return out + + +def load_bed(path, tol=0): + d = collections.defaultdict(list) + for line in open(path): + c, s, e = line.split()[:3] + d[c].append((int(s), int(e))) + return {c: merge(v, tol) for c, v in d.items()} + + +def main(): + priv = collections.defaultdict(list) + for c in CONTIGS: + with open(f'private.{c}.tsv') as f: + next(f) + for line in f: + p = line.split('\t') + priv[p[0]].append(int(p[1])) + for c in priv: + priv[c].sort() + + callable_r = load_bed('callable.bed') + truth = load_bed('truth_HG00096.chm13.bed', 1000) + + def covered(regions, s, e): + tot = 0 + for s2, e2 in regions: + lo, hi = max(s, s2), min(e, e2) + if hi > lo: + tot += hi - lo + return tot + + dens = [] + for c in CONTIGS: + a = priv[c] + lo = min(s for s, _ in callable_r[c]) + hi = max(e for _, e in callable_r[c]) + for b in range(lo, hi, BIN): + e = b + BIN + cal = covered(callable_r[c], b, e) + if cal < BIN * 0.5: + continue # mostly uncallable: not a background sample + if covered(truth.get(c, []), b, e) > 0: + continue # overlaps a real archaic tract: not background + n = bisect.bisect_left(a, e) - bisect.bisect_left(a, b) + dens.append(n / (cal / 1e6)) + dens.sort() + n = len(dens) + mean = sum(dens) / n + var = sum((d - mean) ** 2 for d in dens) / (n - 1) + q = lambda f: dens[int(f * (n - 1))] + + print(f'BACKGROUND private-variant density, {BIN // 1000} kb bins, callable and NOT archaic') + print(f' n bins {n}') + print(f' mean {mean:.0f}/Mb median {q(.5):.0f} p10 {q(.1):.0f} p90 {q(.9):.0f}' + f' max {dens[-1]:.0f}') + print(f' p90/p10 spread = {q(.9) / max(q(.1), 1e-9):.1f}x') + print() + print(f' enrichment inside real archaic tracts, for comparison: 2.89x') + print(f' fraction of BACKGROUND bins already above 2.89x the median: ' + f'{sum(1 for d in dens if d > 2.89 * q(.5)) / n * 100:.1f}%') + print() + # Poisson would predict variance == mean for the per-bin COUNT; overdispersion is the degree to + # which a single-lambda model is simply the wrong distribution. + counts = [d * BIN / 1e6 for d in dens] + cm = sum(counts) / len(counts) + cv = sum((x - cm) ** 2 for x in counts) / (len(counts) - 1) + print(f' per-bin counts: mean {cm:.1f}, variance {cv:.1f} -> overdispersion {cv / cm:.1f}x') + print(' (Poisson assumes variance = mean, i.e. 1.0x. Anything well above that means a') + print(' single-lambda background will call its own upper tail archaic.)') + print() + print(f' sd/mean of background density = {math.sqrt(var) / mean:.2f}') + + +if __name__ == '__main__': + main() diff --git a/scripts/archaic-validation/observable/callset_compare.py b/scripts/archaic-validation/observable/callset_compare.py new file mode 100644 index 0000000..5c05da7 --- /dev/null +++ b/scripts/archaic-validation/observable/callset_compare.py @@ -0,0 +1,100 @@ +"""Same individual, same outgroup, same build: 1000G's call set versus ours. + +The Tier B observable is "variants this person carries that no African carries". Its contrast +inside real archaic tracts is 2.89x for our calls -- too weak to separate from a background that +varies 5.3x. hmmix's documented example emission rates imply a background near 40/Mb and a ~10x +contrast, and ours is 124/Mb at 2.89x. + +Our outgroup track was just verified complete against its source, so if the observable is diluted +the dilution must come from the variant calls. HG00096 is in the 1000G callset, so the same +quantity can be computed from their calls and ours and compared directly. +""" + +import bisect +import collections + + +def merge(iv, tol=0): + iv = sorted(iv) + out, cs, ce = [], *iv[0] + for s, e in iv[1:]: + if s > ce + tol: + out.append((cs, ce)) + cs, ce = s, e + else: + ce = max(ce, e) + out.append((cs, ce)) + return out + + +def load_bed(path, tol=0): + d = collections.defaultdict(list) + for line in open(path): + c, s, e = line.split()[:3] + d[c].append((int(s), int(e))) + return {c: merge(v, tol) for c, v in d.items()} + + +def covered(regions, s, e): + tot = 0 + for s2, e2 in regions: + lo, hi = max(s, s2), min(e, e2) + if hi > lo: + tot += hi - lo + return tot + + +def stats(positions, callable_r, truth, label): + a = sorted(positions) + cnt = lambda s, e: bisect.bisect_left(a, e) - bisect.bisect_left(a, s) + in_bp = sum(covered(callable_r, s, e) for s, e in truth) + in_n = sum(cnt(s, e) for s, e in truth) + cb = sum(e - s for s, e in callable_r) + cn = sum(cnt(s, e) for s, e in callable_r) + out_n, out_bp = cn - in_n, cb - in_bp + enrich = (in_n / in_bp) / (out_n / out_bp) + print(f'{label:26s} {len(a):7d} private in-tract {in_n / (in_bp / 1e6):6.1f}/Mb ' + f'background {out_n / (out_bp / 1e6):6.1f}/Mb CONTRAST {enrich:5.2f}x') + return enrich + + +def main(): + afr = set(int(x) for x in open('afr_sites.chr21.txt')) + kgp_carried = [int(x) for x in open('kgp_carried.chr21.txt')] + ours_carried, ours_private = [], [] + with open('private.chr21.tsv') as f: + next(f) + for line in f: + p = line.split('\t') + ours_private.append(int(p[1])) + + callable_r = load_bed('callable.bed')['chr21'] + truth = load_bed('truth_HG00096.chm13.bed', 1000)['chr21'] + + kgp_private = [p for p in kgp_carried if p not in afr] + + print(f'HG00096 chr21, carried SNVs: 1000G {len(kgp_carried):7d} ours {93643:7d}' + f' ratio {93643 / len(kgp_carried):.2f}x') + print(f'African outgroup: {len(afr)} segregating sites (652 unrelated individuals)\n') + + print('observable count density in tract / background contrast') + e_kgp = stats(kgp_private, callable_r, truth, "1000G's calls") + e_ours = stats(ours_private, callable_r, truth, 'our calls') + print() + print(f' contrast ratio 1000G/ours = {e_kgp / e_ours:.2f}x') + print() + print(' If 1000G\'s calls give a much higher contrast on the SAME person, the observable is') + print(' diluted by our variant calling, not by the outgroup or the model -- and the fix is') + print(' upstream of the HMM entirely.') + + # Where do the extra calls sit? If ours are mostly at positions 1000G did not call at all, + # they are ours to explain. + kgp_set = set(kgp_carried) + extra = [p for p in ours_private if p not in kgp_set] + shared = len(ours_private) - len(extra) + print(f'\n of our {len(ours_private)} private calls: {shared} also called by 1000G, ' + f'{len(extra)} ours alone ({len(extra) / max(len(ours_private), 1) * 100:.0f}%)') + + +if __name__ == '__main__': + main() diff --git a/scripts/archaic-validation/observable/discriminability.py b/scripts/archaic-validation/observable/discriminability.py new file mode 100644 index 0000000..b548439 --- /dev/null +++ b/scripts/archaic-validation/observable/discriminability.py @@ -0,0 +1,73 @@ +"""How detectable is a single tract under each observable? + +Contrast alone does not decide this. What decides it is how much evidence one tract carries, and the +two observables differ 30-fold there: + + private-variant density : ~1 informative variant per 36 kb tract (Poisson 0.6 vs 1.3) + archaic-allele matching : ~30 diagnostic sites (Binomial 30 x 13% vs 30 x 39.5%) + +Same ~3x contrast, wildly different decidability. This computes, for each, the sensitivity +achievable at a fixed false-positive rate -- i.e. what fraction of real tracts a perfect classifier +could call while keeping background calls rare. +""" + +import math + + +def binom_pmf(k, n, p): + return math.comb(n, k) * p ** k * (1 - p) ** (n - k) + + +def pois_pmf(k, lam): + return math.exp(-lam) * lam ** k / math.factorial(k) + + +def roc(bg, sig, kmax): + """Sensitivity at the threshold where background false-positive rate first drops below 5%/1%.""" + out = {} + for target in (0.05, 0.01): + for t in range(kmax + 1): + fp = sum(bg(k) for k in range(t, kmax + 1)) + if fp <= target: + out[target] = (t, sum(sig(k) for k in range(t, kmax + 1)), fp) + break + else: + out[target] = (None, 0.0, 0.0) + return out + + +print('DETECTABILITY OF ONE TRACT (36 kb, the measured median)\n') + +# --- private-variant density --------------------------------------------------------------- +lam_bg, lam_sig = 0.64, 1.30 +r = roc(lambda k: pois_pmf(k, lam_bg), lambda k: pois_pmf(k, lam_sig), 30) +print(f'private-variant density background Poisson({lam_bg}) tract Poisson({lam_sig})') +for tgt, (t, sens, fp) in r.items(): + print(f' at <= {tgt:.0%} false positives: threshold k >= {t}, sensitivity {sens:6.1%}') + +# --- archaic-allele matching --------------------------------------------------------------- +n, p_bg, p_sig = 30, 0.130, 0.395 +r = roc(lambda k: binom_pmf(k, n, p_bg), lambda k: binom_pmf(k, n, p_sig), n) +print(f'\narchaic-allele matching background Binom({n}, {p_bg}) tract Binom({n}, {p_sig})') +for tgt, (t, sens, fp) in r.items(): + print(f' at <= {tgt:.0%} false positives: threshold k >= {t}, sensitivity {sens:6.1%}') + +# What tract size does each need to become usable? +print('\nTRACT SIZE NEEDED FOR 80% SENSITIVITY AT 5% FALSE POSITIVES') +for label, kind in (('private density', 'pois'), ('allele matching', 'binom')): + for kb in (10, 20, 36, 50, 100, 200, 500): + scale = kb / 36 + if kind == 'pois': + bg, sig = lam_bg * scale, lam_sig * scale + rr = roc(lambda k: pois_pmf(k, bg), lambda k: pois_pmf(k, sig), 60) + else: + nn = max(int(round(n * scale)), 1) + rr = roc(lambda k: binom_pmf(k, nn, p_bg), lambda k: binom_pmf(k, nn, p_sig), nn) + if rr[0.05][1] >= 0.80: + print(f' {label:18s} {kb:4d} kb -> sensitivity {rr[0.05][1]:.0%}') + break + else: + print(f' {label:18s} not reached by 500 kb') + +print('\n hmmix p10 tract is 7 kb and the median 31-36 kb, so a method needing hundreds of kb') +print(' cannot report tracts at the resolution the feature claims.') diff --git a/scripts/archaic-validation/observable/haplotype_match.py b/scripts/archaic-validation/observable/haplotype_match.py new file mode 100644 index 0000000..e6033aa --- /dev/null +++ b/scripts/archaic-validation/observable/haplotype_match.py @@ -0,0 +1,121 @@ +"""Should the observable be archaic-ALLELE MATCHING rather than private-mutation density? + +The current caller counts private variants per window: measured contrast inside real archaic tracts +is ~2x, and at ~20-120 private variants/Mb a median 31 kb tract carries about ONE informative +variant. That is why no parameter setting works. + +But an introgressed tract is a haplotype inherited intact from an archaic ancestor, and we HOLD the +archaic genomes -- 2,031,406 diagnostic sites where archaics carry a derived allele. A real tract +should carry the archaic allele at a large fraction of the diagnostic sites it spans, while a +non-introgressed region carries it only at the background rate (~4%, per the design's own +measurement). + +Diagnostic sites are the denominator here, not megabases, which controls for their uneven density +for free -- the thing a mutation-rate map was going to have to correct. +""" + +import bisect +import collections +import json +import random + +CONTIGS = ('chr21', 'chr22') + + +def merge(iv, tol=0): + iv = sorted(iv) + out, cs, ce = [], *iv[0] + for s, e in iv[1:]: + if s > ce + tol: + out.append((cs, ce)) + cs, ce = s, e + else: + ce = max(ce, e) + out.append((cs, ce)) + return out + + +def load_bed(path, tol=0): + d = collections.defaultdict(list) + for line in open(path): + c, s, e = line.split()[:3] + d[c].append((int(s), int(e))) + return {c: merge(v, tol) for c, v in d.items()} + + +def main(): + # Diagnostic sites: position -> (derived base, class) + diag = collections.defaultdict(dict) + with open('classify.tsv') as f: + next(f) + for line in f: + c, p, d, k = line.rstrip('\n').split('\t') + diag[c][int(p)] = (d, int(k)) + + # The subject's calls, so we can ask whether he carries the archaic allele. + carries = collections.defaultdict(dict) + for c in CONTIGS: + for rec in json.load(open(f'HG00096.{c}.calls.json')): + ref = rec['reference_allele'] + alt = rec['alternate_allele'] + dos = rec['dosage'] + if dos is None or dos < 0: + continue + alleles = set() + if dos < 2: + alleles.add(ref) + if dos > 0: + alleles.add(alt) + carries[rec['contig']][rec['position']] = alleles + + callable_r = load_bed('callable.bed') + truth = load_bed('truth_HG00096.chm13.bed', 1000) + + def in_regions(regions, p): + i = bisect.bisect_right([s for s, _ in regions], p) - 1 + return i >= 0 and regions[i][1] > p + + stats = {'in': [0, 0], 'out': [0, 0]} + for c in CONTIGS: + cal = callable_r[c] + tr = truth.get(c, []) + for p, (derived, _k) in diag[c].items(): + if not in_regions(cal, p): + continue + # EVERY diagnostic site in callable territory is in the denominator. The caller emits + # only variant records, so no record means hom-reference — i.e. NOT carrying the + # archaic allele, given the panel orients derived as ALT. Conditioning on "has a call" + # instead samples only sites where he already has a variant, which is why that version + # reported an impossible ~80% carrying rate against a known 4.3% background. + called = carries[c].get(p) + bucket = 'in' if in_regions(tr, p) else 'out' + stats[bucket][1] += 1 + if called is not None and derived in called: + stats[bucket][0] += 1 + + print('ARCHAIC-ALLELE CARRYING RATE at diagnostic sites the subject has a call for') + for b, label in (('in', "inside hmmix's tracts"), ('out', 'elsewhere')): + hit, tot = stats[b] + rate = hit / tot * 100 if tot else 0 + print(f' {label:26s} {hit:6d} / {tot:6d} = {rate:5.1f}%') + ri = stats['in'][0] / max(stats['in'][1], 1) + ro = stats['out'][0] / max(stats['out'][1], 1) + print(f' ENRICHMENT {ri / ro if ro else 0:.2f}x') + print() + print(' Compare: private-variant DENSITY gives 2.89x with ~1 informative variant per tract.') + print(' Sites here are the denominator, so uneven diagnostic-site density cancels out.') + + # How much evidence does a typical tract actually carry under this observable? + spans = [e - s for c in CONTIGS for s, e in truth.get(c, [])] + spans.sort() + med = spans[len(spans) // 2] + dens = sum(len(diag[c]) for c in CONTIGS) / sum( + e - s for c in CONTIGS for s, e in merge([(min(x for x, _ in callable_r[c]), + max(y for _, y in callable_r[c]))])) + print(f'\n median tract {med / 1000:.0f} kb; diagnostic sites ~{dens * 1e6:.0f}/Mb' + f' -> ~{med * dens:.0f} informative sites per tract') + print(' (the density model gets ~1 private variant for the same tract)') + + +if __name__ == '__main__': + main() diff --git a/scripts/archaic-validation/observable/mutrate_check.py b/scripts/archaic-validation/observable/mutrate_check.py new file mode 100644 index 0000000..05e6105 --- /dev/null +++ b/scripts/archaic-validation/observable/mutrate_check.py @@ -0,0 +1,132 @@ +"""Does African-outgroup site density explain the background variation the emission model ignores? + +If it does, it is a usable mutation-rate normalizer and the fix is an asset build. If it does not, +a different proxy is needed and building this one would waste the effort. + +Two questions, in order: + 1. Correlation between outgroup density and private-variant density in BACKGROUND regions + (callable, non-archaic). This is the normalizer's whole job. + 2. Whether normalizing by it actually flattens the background — the overdispersion should fall + from the measured 14.6x toward 1.0x. Correlation alone does not guarantee that. +""" + +import bisect +import collections +import math + +BIN = 100_000 +CONTIGS = ('chr21', 'chr22') + + +def merge(iv, tol=0): + iv = sorted(iv) + out, cs, ce = [], *iv[0] + for s, e in iv[1:]: + if s > ce + tol: + out.append((cs, ce)) + cs, ce = s, e + else: + ce = max(ce, e) + out.append((cs, ce)) + return out + + +def load_bed(path, tol=0): + d = collections.defaultdict(list) + for line in open(path): + c, s, e = line.split()[:3] + d[c].append((int(s), int(e))) + return {c: merge(v, tol) for c, v in d.items()} + + +def covered(regions, s, e): + tot = 0 + for s2, e2 in regions: + lo, hi = max(s, s2), min(e, e2) + if hi > lo: + tot += hi - lo + return tot + + +def pearson(x, y): + n = len(x) + mx, my = sum(x) / n, sum(y) / n + num = sum((a - mx) * (b - my) for a, b in zip(x, y)) + dx = math.sqrt(sum((a - mx) ** 2 for a in x)) + dy = math.sqrt(sum((b - my) ** 2 for b in y)) + return num / (dx * dy) if dx and dy else 0.0 + + +def main(): + priv = collections.defaultdict(list) + for c in CONTIGS: + with open(f'private.{c}.tsv') as f: + next(f) + for line in f: + p = line.split('\t') + priv[p[0]].append(int(p[1])) + for c in priv: + priv[c].sort() + + og = collections.defaultdict(collections.Counter) + with open('og_density.tsv') as f: + next(f) + for line in f: + c, w, k = line.split('\t') + og[c][int(w)] = int(k) + + callable_r = load_bed('callable.bed') + truth = load_bed('truth_HG00096.chm13.bed', 1000) + + xs, ys, weights = [], [], [] + for c in CONTIGS: + a = priv[c] + lo = min(s for s, _ in callable_r[c]) + hi = max(e for _, e in callable_r[c]) + for b in range(lo, hi, BIN): + e = b + BIN + cal = covered(callable_r[c], b, e) + if cal < BIN * 0.5 or covered(truth.get(c, []), b, e) > 0: + continue + n = bisect.bisect_left(a, e) - bisect.bisect_left(a, b) + ogn = sum(og[c].get(w, 0) for w in range(b, e, 1000)) + if ogn == 0: + continue + xs.append(ogn / (cal / 1e6)) # outgroup sites per callable Mb + ys.append(n / (cal / 1e6)) # our private variants per callable Mb + weights.append(cal) + + r = pearson(xs, ys) + print(f'BACKGROUND bins ({BIN // 1000} kb, callable, non-archaic): n = {len(xs)}') + print(f' outgroup-site density vs private-variant density: Pearson r = {r:+.3f}' + f' (r^2 = {r * r:.2f})') + print(f' -> outgroup density explains {r * r * 100:.0f}% of the background variance') + print() + + # Does dividing by the proxy actually flatten it? + mean_x = sum(xs) / len(xs) + raw = ys + norm = [y / (x / mean_x) for x, y in zip(xs, ys)] + + def stats(v, label): + v = sorted(v) + n = len(v) + m = sum(v) / n + q = lambda f: v[int(f * (n - 1))] + # overdispersion of the implied per-bin count + cnt = [d * BIN / 1e6 for d in v] + cm = sum(cnt) / len(cnt) + cv = sum((z - cm) ** 2 for z in cnt) / (len(cnt) - 1) + print(f' {label:22s} p10 {q(.1):6.0f} median {q(.5):6.0f} p90 {q(.9):6.0f}' + f' p90/p10 {q(.9) / max(q(.1), 1e-9):5.1f}x overdispersion {cv / cm:5.1f}x') + + print('BACKGROUND FLATNESS (the emission model assumes this is flat):') + stats(raw, 'raw density') + stats(norm, 'normalized by proxy') + print() + print(' The archaic signal is 2.89x. For the model to separate it, the background spread') + print(' must be well below that. Overdispersion near 1.0x is what a Poisson emission needs.') + + +if __name__ == '__main__': + main() diff --git a/scripts/archaic-validation/observable/observable_compare.py b/scripts/archaic-validation/observable/observable_compare.py new file mode 100644 index 0000000..e5f7799 --- /dev/null +++ b/scripts/archaic-validation/observable/observable_compare.py @@ -0,0 +1,123 @@ +"""Compare candidate observables for the Tier B HMM on signal-to-noise. + +The current observable is "any private variant". Measured, that gives 2.89x enrichment inside real +archaic tracts against a background that varies 5.3x (p10-p90) and is 14.6x overdispersed — the +noise is bigger than the signal, which is why no parameter setting works. + +A more specific observable should trade count for contrast: fewer observations, but a much higher +ratio inside tracts. What matters is whether the enrichment clears the background spread, because +that ratio is what any two-state model has to separate. +""" + +import bisect +import collections + +BIN = 100_000 +CONTIGS = ('chr21', 'chr22') + + +def merge(iv, tol=0): + iv = sorted(iv) + out, cs, ce = [], *iv[0] + for s, e in iv[1:]: + if s > ce + tol: + out.append((cs, ce)) + cs, ce = s, e + else: + ce = max(ce, e) + out.append((cs, ce)) + return out + + +def load_bed(path, tol=0): + d = collections.defaultdict(list) + for line in open(path): + c, s, e = line.split()[:3] + d[c].append((int(s), int(e))) + return {c: merge(v, tol) for c, v in d.items()} + + +def covered(regions, s, e): + tot = 0 + for s2, e2 in regions: + lo, hi = max(s, s2), min(e, e2) + if hi > lo: + tot += hi - lo + return tot + + +def assess(pos_by_contig, callable_r, truth, label): + in_n = in_bp = out_n = out_bp = 0 + dens = [] + for c in CONTIGS: + a = sorted(pos_by_contig.get(c, [])) + cnt = lambda s, e: bisect.bisect_left(a, e) - bisect.bisect_left(a, s) + in_bp += sum(covered(callable_r[c], s, e) for s, e in truth[c]) + in_n += sum(cnt(s, e) for s, e in truth[c]) + cb = sum(e - s for s, e in callable_r[c]) + out_bp += cb - sum(covered(callable_r[c], s, e) for s, e in truth[c]) + out_n += sum(cnt(s, e) for s, e in callable_r[c]) - sum(cnt(s, e) for s, e in truth[c]) + lo = min(s for s, _ in callable_r[c]) + hi = max(e for _, e in callable_r[c]) + for b in range(lo, hi, BIN): + e = b + BIN + cal = covered(callable_r[c], b, e) + if cal < BIN * 0.5 or covered(truth.get(c, []), b, e) > 0: + continue + dens.append(cnt(b, e) / (cal / 1e6)) + if not in_bp or not out_n: + print(f'{label:38s} (insufficient data)') + return + enrich = (in_n / in_bp) / (out_n / out_bp) + dens.sort() + n = len(dens) + q = lambda f: dens[int(f * (n - 1))] + spread = q(.9) / q(.1) if q(.1) > 0 else float('inf') + print(f'{label:38s} n={in_n + out_n:6d} in-tract {in_n / (in_bp / 1e6):6.1f}/Mb ' + f'bg {out_n / (out_bp / 1e6):6.1f}/Mb ENRICH {enrich:5.2f}x bg p90/p10 {spread:6.1f}x') + + +def main(): + priv = collections.defaultdict(list) + for c in CONTIGS: + with open(f'private.{c}.tsv') as f: + next(f) + for line in f: + p = line.rstrip('\n').split('\t') + priv[p[0]].append((int(p[1]), int(p[2]), int(p[3]), int(p[4]))) + + diag = collections.defaultdict(dict) + with open('classify.tsv') as f: + next(f) + for line in f: + c, p, d, k = line.rstrip('\n').split('\t') + diag[c][int(p)] = (d, int(k)) + + callable_r = load_bed('callable.bed') + truth = load_bed('truth_HG00096.chm13.bed', 1000) + + print('observable count density in/out contrast') + assess({c: [r[0] for r in v] for c, v in priv.items()}, callable_r, truth, + 'all private variants (current)') + + # Private variants that also sit at a known archaic-diagnostic site. + d_all = {c: [r[0] for r in v if r[0] in diag.get(c, {})] for c, v in priv.items()} + assess(d_all, callable_r, truth, 'private AND archaic-diagnostic') + + for k, name in ((0, 'Neanderthal-diagnostic'), (2, 'shared-archaic')): + f = {c: [r[0] for r in v if diag.get(c, {}).get(r[0], (None, -1))[1] == k] + for c, v in priv.items()} + assess(f, callable_r, truth, f' ...of which {name}') + + # Control: diagnostic sites the subject does NOT carry should show no enrichment. + carried = {c: {r[0] for r in v} for c, v in priv.items()} + notc = {c: [p for p in diag.get(c, {}) if p not in carried.get(c, set())] for c in CONTIGS} + assess(notc, callable_r, truth, 'diagnostic sites NOT carried (control)') + + print() + print(' Contrast is what matters: the enrichment has to clear the background spread for a') + print(' two-state model to separate the states. "all private" fails that (2.89x vs 5.3x).') + + +if __name__ == '__main__': + main() diff --git a/scripts/archaic-validation/observable/quality_effect.py b/scripts/archaic-validation/observable/quality_effect.py new file mode 100644 index 0000000..23f008c --- /dev/null +++ b/scripts/archaic-validation/observable/quality_effect.py @@ -0,0 +1,123 @@ +"""Is the background's excess variance real, or is it this caller's error rate varying by region? + +The emission model needs a flat background. Measured raw, it is 14.6x overdispersed; normalizing by +African-outgroup density gets it to 7.4x, still above the 2.89x archaic signal. If the remainder is +low-confidence calls -- artifacts clustering in hard regions -- then filtering should flatten the +background AND raise the archaic enrichment, because artifacts dilute real tracts too. + +Reports both numbers for each filter, since a filter that flattens the background by discarding the +signal has bought nothing. +""" + +import bisect +import collections + +BIN = 100_000 +CONTIGS = ('chr21', 'chr22') + + +def merge(iv, tol=0): + iv = sorted(iv) + out, cs, ce = [], *iv[0] + for s, e in iv[1:]: + if s > ce + tol: + out.append((cs, ce)) + cs, ce = s, e + else: + ce = max(ce, e) + out.append((cs, ce)) + return out + + +def load_bed(path, tol=0): + d = collections.defaultdict(list) + for line in open(path): + c, s, e = line.split()[:3] + d[c].append((int(s), int(e))) + return {c: merge(v, tol) for c, v in d.items()} + + +def covered(regions, s, e): + tot = 0 + for s2, e2 in regions: + lo, hi = max(s, s2), min(e, e2) + if hi > lo: + tot += hi - lo + return tot + + +def load_private(): + rows = collections.defaultdict(list) + for c in CONTIGS: + with open(f'private.{c}.tsv') as f: + next(f) + for line in f: + p = line.rstrip('\n').split('\t') + rows[p[0]].append((int(p[1]), int(p[2]), int(p[3]), int(p[4]))) + return rows + + +def assess(rows, callable_r, truth, label): + pos = {c: sorted(r[0] for r in v) for c, v in rows.items()} + in_n = in_bp = out_n = out_bp = 0 + dens = [] + for c in CONTIGS: + a = pos.get(c, []) + cnt = lambda s, e: bisect.bisect_left(a, e) - bisect.bisect_left(a, s) + tb = sum(covered(callable_r[c], s, e) for s, e in truth[c]) + tn = sum(cnt(s, e) for s, e in truth[c]) + cb = sum(e - s for s, e in callable_r[c]) + cn = sum(cnt(s, e) for s, e in callable_r[c]) + in_n += tn + in_bp += tb + out_n += cn - tn + out_bp += cb - tb + lo = min(s for s, _ in callable_r[c]) + hi = max(e for _, e in callable_r[c]) + for b in range(lo, hi, BIN): + e = b + BIN + cal = covered(callable_r[c], b, e) + if cal < BIN * 0.5 or covered(truth.get(c, []), b, e) > 0: + continue + dens.append(cnt(b, e) / (cal / 1e6)) + enrich = (in_n / in_bp) / (out_n / out_bp) if out_n and in_bp else 0 + dens.sort() + n = len(dens) + q = lambda f: dens[int(f * (n - 1))] + counts = [d * BIN / 1e6 for d in dens] + cm = sum(counts) / len(counts) + cv = sum((z - cm) ** 2 for z in counts) / (len(counts) - 1) + kept = sum(len(v) for v in rows.values()) + print(f'{label:26s} kept {kept:6d} enrich {enrich:4.2f}x ' + f'bg p90/p10 {q(.9) / max(q(.1), 1e-9):5.1f}x overdisp {cv / max(cm, 1e-9):5.1f}x') + + +def main(): + raw = load_private() + callable_r = load_bed('callable.bed') + truth = load_bed('truth_HG00096.chm13.bed', 1000) + + print('filter variants archaic signal background flatness') + assess(raw, callable_r, truth, 'none (current)') + for gq in (20, 30, 50): + f = {c: [r for r in v if r[2] >= gq] for c, v in raw.items()} + assess(f, callable_r, truth, f'GQ >= {gq}') + for dp in (10, 15, 20): + f = {c: [r for r in v if r[3] >= dp] for c, v in raw.items()} + assess(f, callable_r, truth, f'depth >= {dp}') + # Homozygous-derived only: an introgressed tract is usually heterozygous, so this should HURT + # if the signal is real -- a useful control that the enrichment is not an artifact of genotype. + f = {c: [r for r in v if r[1] == 1] for c, v in raw.items()} + assess(f, callable_r, truth, 'het only (control)') + f = {c: [r for r in v if r[1] == 2] for c, v in raw.items()} + assess(f, callable_r, truth, 'hom-alt only (control)') + for gq, dp in ((30, 15), (50, 20)): + f = {c: [r for r in v if r[2] >= gq and r[3] >= dp] for c, v in raw.items()} + assess(f, callable_r, truth, f'GQ>={gq} & depth>={dp}') + print() + print(' Wanted: enrichment UP and overdispersion DOWN together. Overdispersion must fall well') + print(' below the enrichment for a Poisson emission to separate the two states.') + + +if __name__ == '__main__': + main() diff --git a/scripts/archaic-validation/observable/truth_selfcheck.py b/scripts/archaic-validation/observable/truth_selfcheck.py new file mode 100644 index 0000000..97be171 --- /dev/null +++ b/scripts/archaic-validation/observable/truth_selfcheck.py @@ -0,0 +1,120 @@ +"""Is my lifted truth actually where hmmix says it is? + +Everything downstream rests on it: I concluded the segment caller's locations are below chance, and +a shipped feature was gated off on that basis. But measuring private-variant density inside those +tracts gives only ~2x enrichment even in 1000G's OWN calls -- and hmmix found these tracts BY that +density, so it should be far higher. Either their method is weaker than advertised, or my lift put +the tracts in the wrong place. + +This checks the truth in native hg38, with no lifting anywhere, using only hmmix's own two files: +their segment calls and their DAV (derived-in-archaic) SNP list. If their segments are enriched for +their own archaic SNPs in hg38, the callset is self-consistent and any failure is in my lift. +""" + +import bisect +import collections +import os +import random + +SEG = os.path.expanduser('~/.decodingus/ancestry-build/tmp/hmmix_segments_chr21_22.tsv') +SNPS = os.path.expanduser('~/.decodingus/ancestry-build/raw/hmmix/hg38_1000g_SNPS.txt') +SAMPLE = 'HG00096' +CONTIGS = ('chr21', 'chr22') + + +def merge(iv, tol=0): + iv = sorted(iv) + out, cs, ce = [], *iv[0] + for s, e in iv[1:]: + if s > ce + tol: + out.append((cs, ce)) + cs, ce = s, e + else: + ce = max(ce, e) + out.append((cs, ce)) + return out + + +def main(): + # hmmix's own segments for this sample, in hg38, unioned across haplotypes. + seg = collections.defaultdict(list) + seg_all = collections.defaultdict(list) + with open(SEG) as f: + next(f) + for line in f: + p = line.rstrip('\n').split('\t') + if p[4] not in CONTIGS: + continue + seg_all[p[4]].append((int(p[5]), int(p[6]))) + if p[0] == SAMPLE: + seg[p[4]].append((int(p[5]), int(p[6]))) + S = {c: merge(v) for c, v in seg.items()} + + # hmmix's own archaic SNP list, in hg38. + dav = collections.defaultdict(list) + with open(SNPS) as f: + next(f) + for line in f: + p = line.split('\t') + if p[0] in CONTIGS: + dav[p[0]].append(int(p[1])) + for c in dav: + dav[c].sort() + + print(f'hmmix segments for {SAMPLE} (hg38, unioned): ' + f'{sum(len(v) for v in S.values())} tracts, ' + f'{sum(e - s for v in S.values() for s, e in v) / 1e6:.3f} Mb') + print(f'hmmix DAV archaic SNPs on chr21+22: {sum(len(v) for v in dav.values())}') + print() + + # Span over which everyone's tracts fall -- the region the sampling actually covers. + bounds = {c: (min(s for s, _ in v), max(e for _, e in v)) for c, v in seg_all.items()} + + tot_in = tot_bp = 0 + tot_out = tot_out_bp = 0 + for c in CONTIGS: + a = dav[c] + cnt = lambda s, e: bisect.bisect_left(a, e) - bisect.bisect_left(a, s) + lo, hi = bounds[c] + in_n = sum(cnt(s, e) for s, e in S[c]) + in_bp = sum(e - s for s, e in S[c]) + all_n = cnt(lo, hi) + all_bp = hi - lo + tot_in += in_n + tot_bp += in_bp + tot_out += all_n - in_n + tot_out_bp += all_bp - in_bp + + d_in = tot_in / (tot_bp / 1e6) + d_out = tot_out / (tot_out_bp / 1e6) + print('IN NATIVE hg38, NO LIFTING:') + print(f' archaic SNPs inside {SAMPLE}\'s own hmmix tracts : {d_in:7.1f} /Mb ({tot_in} in ' + f'{tot_bp / 1e6:.2f} Mb)') + print(f' archaic SNPs elsewhere : {d_out:7.1f} /Mb') + print(f' ENRICHMENT : {d_in / d_out:6.2f}x') + print() + print(' hmmix called these tracts from archaic-variant density, so a strong enrichment here') + print(' means their callset is internally consistent and my LIFT is what to distrust.') + print(' A weak enrichment here means the truth was never as sharp as assumed.') + + # Null: same tract lengths placed at random in the same span. + rnd = random.Random(0) + vals = [] + for _ in range(200): + n = bp = 0 + for c in CONTIGS: + a = dav[c] + lo, hi = bounds[c] + for s, e in S[c]: + L = e - s + p = rnd.randint(lo, max(lo, hi - L)) + n += bisect.bisect_left(a, p + L) - bisect.bisect_left(a, p) + bp += L + vals.append(n / (bp / 1e6) / d_out) + vals.sort() + print(f' null (same tracts placed at random): mean {sum(vals) / len(vals):.2f}x ' + f'p95 {vals[189]:.2f}x') + + +if __name__ == '__main__': + main() From 7a7f743c0aeba576722ebfa04c2e2901c188ca2e Mon Sep 17 00:00:00 2001 From: James Kane Date: Fri, 31 Jul 2026 07:24:42 -0500 Subject: [PATCH 2/9] Archaic Tier B v2: call tracts by matching the archaic genomes, not counting mutations First increment of the rework diagnosed in the previous commit. New module `archaic_match`, its probe, and the harness fix needed to score it honestly. Not yet wired to the app; Tier B stays gated. THE MODEL. An introgressed tract is a haplotype inherited intact from an archaic ancestor, so it carries the archaic allele at a large share of the diagnostic sites it spans. That is a two-state HMM whose observation is ONE BIT PER DIAGNOSTIC SITE -- carried or not -- with Bernoulli emissions, indexed BY SITE rather than by base pair, and recombination-scaled transitions between consecutive sites. Indexing by site is what makes this robust where density was not: diagnostic sites become the denominator, so their uneven density cancels. The mutation-rate map the density model needed -- and which no available proxy supplied, the best explaining 38% of a 14.6x overdispersion -- is simply not required. MEASURED ON REAL DATA, HG00096 chr21+22, against hmmix's own calls: density (shipped) matching (this) sensitivity 2.1% 38.7% precision 1.5% 19.8% Jaccard 0.009 0.151 tracts hit 4/47 16/47 random-placement null 5.0% 6.5% (p95 12.3, max 20.2) verdict AT OR BELOW CHANCE ABOVE THE NULL'S FULL RANGE 38.7% is 5.9x the null mean and higher than any of 400 random draws. The same data, the same truth, the same harness -- only the observable changed. TWO CHOICES WORTH REVIEWING: - Parameters are MEASURED, not fitted. p_background is estimated directly from the subject's genome-wide carrying rate (dominated by non-archaic sequence, so it estimates the background); p_archaic derives from it by a measured 3.04x ratio. Unconstrained Baum-Welch on the previous caller diverged to a degenerate fit -- 22x emission ratio, 9 kb tracts, 7x the truth -- so EM is deliberately not used here. - Sites where the archaic-derived allele IS the reference base are DROPPED. Every reference-matching genome trivially carries the derived allele there, so the site separates nothing; worse, because the caller emits only variant records, a no-call at such a site means the subject DOES carry it, the opposite of what a no-call means everywhere else. Getting this wrong is not hypothetical: an early version of the analysis conditioned on "has a call" and reported an 80% carrying rate against a known 4.3% background. Five unit tests cover the model without assets or I/O: a real run is found, a background-only contig calls nothing, scattered background-rate carriers do NOT accumulate into a tract (the density caller's defining failure, restated in this observable), reference-derived sites are dropped, and a missing call is a hom-reference non-carrier. HARNESS FIX, without which the numbers above are unreadable: compare_locations.py now computes the random-placement null itself and merges the truth with a 1 kb tolerance. Sensitivity rises with how much sequence you call, so a bare figure means nothing -- the density caller's 2.1% looked like weak performance when it was below chance. The tolerance closes the median 2 bp gaps the lift leaves, which otherwise report 423 shards where there are 47 real tracts. STILL OPEN, and none of it is hidden by the numbers above: we over-call by 1.95x (4.56 Mb against 2.34), precision is 19.8%, and this is n=1. Calibration is now legitimate because the locations are real, but it must be fitted on some individuals and validated on held-out ones -- fitting to a cohort statistic is exactly what produced the last caller. Attribution stays off. Co-Authored-By: Claude Opus 5 (1M context) --- .../examples/archaic_match_probe.rs | 86 +++ .../navigator-analysis/src/archaic_match.rs | 495 ++++++++++++++++++ crates/navigator-analysis/src/lib.rs | 1 + .../archaic-validation/compare_locations.py | 49 +- 4 files changed, 623 insertions(+), 8 deletions(-) create mode 100644 crates/navigator-analysis/examples/archaic_match_probe.rs create mode 100644 crates/navigator-analysis/src/archaic_match.rs diff --git a/crates/navigator-analysis/examples/archaic_match_probe.rs b/crates/navigator-analysis/examples/archaic_match_probe.rs new file mode 100644 index 0000000..84944ee --- /dev/null +++ b/crates/navigator-analysis/examples/archaic_match_probe.rs @@ -0,0 +1,86 @@ +//! Run the reference-based archaic tract caller ([`archaic_match`]) on real cached calls. +//! +//! The unit tests prove the model behaves on synthetic runs; this is what shows whether it finds +//! REAL tracts. Emits segments as JSON for `scripts/archaic-validation/compare_locations.py`, which +//! scores them against an external callset and — critically — against the random-placement null the +//! density caller failed. +//! +//! ```sh +//! cargo run --release -p navigator-analysis --example archaic_match_probe -- \ +//! classify.bin callable.bin chm13v2.0.fa genetic_map.bin calls.chr21.json calls.chr22.json \ +//! > segments.json +//! ``` + +use std::collections::BTreeMap; + +use navigator_analysis::archaic::{ArchaicCallable, ArchaicClassify}; +use navigator_analysis::archaic_match::{call_from_observations, observations_for_contig, MatchConfig, SiteObs}; +use navigator_analysis::caller::SiteGenotype; +use navigator_analysis::ibd::GeneticMap; +use navigator_analysis::reader::read_contig_sequence; + +fn main() -> Result<(), Box> { + let a: Vec = std::env::args().skip(1).collect(); + if a.len() < 5 { + eprintln!( + "usage: archaic_match_probe \ + [calls.json ...]" + ); + std::process::exit(2); + } + let classify = ArchaicClassify::from_bytes(&std::fs::read(&a[0])?).map_err(|e| e.to_string())?; + let callable = ArchaicCallable::from_bytes(&std::fs::read(&a[1])?).map_err(|e| e.to_string())?; + let reference = std::path::PathBuf::from(&a[2]); + + let mut calls: Vec = Vec::new(); + for p in &a[4..] { + let mut v: Vec = serde_json::from_str(&std::fs::read_to_string(p)?)?; + calls.append(&mut v); + } + let mut by_contig: BTreeMap> = BTreeMap::new(); + for c in &calls { + by_contig.entry(c.contig.clone()).or_default().insert(c.position, c); + } + eprintln!("{} calls over {} contig(s)", calls.len(), by_contig.len()); + + let mut observations: BTreeMap> = BTreeMap::new(); + let mut lengths: Vec<(String, i32)> = Vec::new(); + for (contig, pos_map) in &by_contig { + // The reference base decides which diagnostic sites are informative at all, so it is read + // rather than assumed (see `observations_for_contig`). + let seq = read_contig_sequence(&reference, contig)?; + let obs = observations_for_contig( + contig, + &classify, + pos_map, + |p| seq.get((p - 1).max(0) as usize).copied().map(|b| b.to_ascii_uppercase()), + &callable, + 0.5, + ); + let carried = obs.iter().filter(|o| o.carries).count(); + eprintln!( + "{contig}: {} informative diagnostic sites, {carried} carried ({:.1}%)", + obs.len(), + if obs.is_empty() { 0.0 } else { carried as f64 * 100.0 / obs.len() as f64 } + ); + lengths.push((contig.clone(), seq.len() as i32)); + observations.insert(contig.clone(), obs); + } + + let gmap = if a[3] == "-" { + let pairs: Vec<(&str, i32)> = lengths.iter().map(|(c, l)| (c.as_str(), *l)).collect(); + eprintln!("genetic map: uniform 1 cM/Mb"); + GeneticMap::uniform(1.0, &pairs) + } else { + GeneticMap::from_bytes(&std::fs::read(&a[3])?).map_err(|e| e.to_string())? + }; + + let result = call_from_observations(&observations, &gmap, &callable, &MatchConfig::default()); + let s = &result.summary; + eprintln!( + "SEGMENTS {} total {:.3} Mb = {:.2}% of {:.1} Mb callable", + s.n_segments, s.total_mb, s.pct_callable, s.callable_mb + ); + println!("{}", serde_json::to_string(&result)?); + Ok(()) +} diff --git a/crates/navigator-analysis/src/archaic_match.rs b/crates/navigator-analysis/src/archaic_match.rs new file mode 100644 index 0000000..8428fc5 --- /dev/null +++ b/crates/navigator-analysis/src/archaic_match.rs @@ -0,0 +1,495 @@ +//! Tier B, second attempt — archaic tracts by **matching the archaic genomes**, not by counting +//! mutations. +//! +//! # Why this replaces the density caller +//! +//! [`crate::archaic_segments`] follows Skov 2018 (hmmix): strip variants Africans also carry, then +//! look for regions dense in what remains. That method exists for people who do **not** have archaic +//! reference genomes and must infer them indirectly. We have all four, and already ship +//! [`ArchaicClassify`] — 2,031,406 sites where the archaics carry a derived allele. +//! +//! Measured on a real European against hmmix's own calls for the same person, the difference is not +//! subtle. Both observables carry the same ~3x contrast, but they differ 30-fold in how much +//! evidence one tract holds, and that is what decides whether a tract can be called at all: +//! +//! | observable | evidence per 36 kb tract | sensitivity at 5 % false positives | +//! |---|---|---| +//! | private-variant density | ~1 variant | 14.3 % | +//! | archaic-allele matching (this) | ~30 sites | 95.1 % | +//! +//! Density does not reach 80 % sensitivity at **500 kb**; matching reaches 95 % at the real median +//! tract of 36 kb. See `documents/design/ArchaicAncestry_Design.md` § *Why it failed*. +//! +//! # The model +//! +//! An introgressed tract is a haplotype inherited intact from an archaic ancestor, so it carries the +//! archaic allele at a large share of the diagnostic sites it spans; elsewhere the subject carries +//! them only at the background rate. That is a two-state HMM whose observation is one **bit per +//! diagnostic site** — carried or not — with Bernoulli emissions, indexed **by site rather than by +//! base pair**. +//! +//! Indexing by site is what makes this robust where the density model was not. Diagnostic sites +//! become the denominator, so their uneven density cancels out: the mutation-rate map the density +//! model needed (and which no available proxy supplied — the best explained 38 % of a 14.6x +//! overdispersion) is simply not required here. +//! +//! Transitions stay recombination-scaled between consecutive sites, as in [`crate::roh`] and the +//! chromosome painter. + +use std::collections::BTreeMap; + +use serde::{Deserialize, Serialize}; + +use crate::archaic::{ArchaicCallable, ArchaicClassify, DiagnosticClass}; +use crate::archaic_segments::{ArchaicSegment, ArchaicSegmentResult, ArchaicSource, ArchaicSummary}; +use crate::caller::SiteGenotype; +use crate::ibd::GeneticMap; + +/// One diagnostic site, reduced to what the HMM consumes. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct SiteObs { + pub position: i64, + /// Whether the subject carries the archaic-derived allele here. + pub carries: bool, + pub class: DiagnosticClass, +} + +/// Tuning knobs. +#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] +pub struct MatchConfig { + /// Rate at which a **non-introgressed** genome carries the archaic allele at a diagnostic site. + /// `None` estimates it from the subject's own genome-wide rate, which absorbs coverage, call + /// behaviour and ancestry. + /// + /// Estimated rather than fixed because it is the denominator of the whole inference — but + /// estimated *directly*, not by EM. Unconstrained Baum-Welch on the previous caller diverged to + /// a degenerate fit (a 22x emission ratio and 9 kb tracts, calling 7x the truth), so parameters + /// here are measured, not fitted. + pub p_background: Option, + /// Rate inside an introgressed tract. `None` derives it as `p_background * archaic_ratio`. + pub p_archaic: Option, + /// Multiple of the background rate expected inside a tract when `p_archaic` is `None`. + /// Measured at 3.04x (39.5 % inside real tracts against 13.0 % elsewhere). + pub archaic_ratio: f64, + /// Expected state switches per centimorgan. + pub switches_per_cm: f64, + /// Discard tracts whose mean posterior is below this. + pub min_posterior: f64, + /// Minimum diagnostic sites in a tract. A tract resting on one or two sites is exactly the + /// failure mode of the density caller, restated in a new observable. + pub min_sites: usize, + /// Discard tracts shorter than this. + pub min_segment_bp: i64, + /// Minimum callable fraction for a site's window to be used at all. + pub min_callable_fraction: f64, + /// Whether to attempt per-segment Neanderthal/Denisovan attribution. Default `false`, unchanged + /// from the density caller: the lineage signal has not been shown to work, and this module does + /// not by itself change that. + pub attribute_lineage: bool, +} + +impl Default for MatchConfig { + fn default() -> Self { + MatchConfig { + p_background: None, + p_archaic: None, + archaic_ratio: 3.04, + switches_per_cm: 1.0, + min_posterior: 0.80, + min_sites: 8, + min_segment_bp: 5_000, + min_callable_fraction: 0.5, + attribute_lineage: false, + } + } +} + +/// Reduce one contig's diagnostic sites to observations. +/// +/// `ref_base` supplies the reference base at a position; sites where the archaic-derived allele +/// **is** the reference base are dropped. At such a site every reference-matching genome trivially +/// "carries" the derived allele, so it separates nothing and would dilute the contrast — and, +/// because the caller emits only variant records, a no-call there means the subject *does* carry it, +/// the opposite of what a no-call means everywhere else. +/// +/// A site with no variant record is hom-reference, hence **not** carrying. Restricting instead to +/// sites where the subject happens to have a call is the trap that made an early version of this +/// analysis report an 80 % carrying rate against a known 4.3 % background: it samples only sites +/// where a variant already exists. +pub fn observations_for_contig( + contig: &str, + classify: &ArchaicClassify, + calls_by_pos: &BTreeMap, + ref_base: impl Fn(i64) -> Option, + callable: &ArchaicCallable, + min_callable_fraction: f64, +) -> Vec { + let Some(c) = classify.contigs.iter().find(|c| c.positions.contig == contig) else { + return Vec::new(); + }; + let mut out = Vec::new(); + for (i, pos) in c.positions.iter().enumerate() { + let Some(&derived) = c.derived.get(i) else { continue }; + // Uninformative: the reference already carries the archaic allele. + if ref_base(pos) == Some(derived) { + continue; + } + if callable.callable_fraction(contig, pos) < min_callable_fraction { + continue; + } + let carries = calls_by_pos.get(&pos).is_some_and(|g| { + g.dosage > 0 && g.alternate_allele.as_bytes().first() == Some(&derived) + }); + let class = match c.classes.get(i).copied().unwrap_or(2) { + 0 => DiagnosticClass::Neanderthal, + 1 => DiagnosticClass::Denisovan, + _ => DiagnosticClass::SharedArchaic, + }; + out.push(SiteObs { + position: pos, + carries, + class, + }); + } + out +} + +fn ln(x: f64) -> f64 { + x.max(1e-300).ln() +} + +fn ln_sum_exp(a: f64, b: f64) -> f64 { + if a == f64::NEG_INFINITY { + return b; + } + if b == f64::NEG_INFINITY { + return a; + } + let m = a.max(b); + m + ((a - m).exp() + (b - m).exp()).ln() +} + +/// Posterior probability of the archaic state at each observation. +/// +/// Log-space forward/backward with recombination-scaled transitions, as in [`crate::roh`]. Exposed +/// so the decoding can be tested against hand-computed posteriors without constructing assets. +pub fn posteriors(obs: &[SiteObs], contig: &str, gmap: &GeneticMap, p_bg: f64, p_arch: f64, switches_per_cm: f64) -> Vec { + let n = obs.len(); + if n == 0 { + return Vec::new(); + } + let emit = |i: usize| -> [f64; 2] { + if obs[i].carries { + [ln(p_bg), ln(p_arch)] + } else { + [ln(1.0 - p_bg), ln(1.0 - p_arch)] + } + }; + // Switch probability between consecutive sites, from the genetic distance between them. + let sw = |i: usize| -> f64 { + let cm = gmap + .interval_cm(contig, obs[i].position as i32, obs[i + 1].position as i32) + .unwrap_or_else(|| (obs[i + 1].position - obs[i].position).max(0) as f64 / 1_000_000.0); + (1.0 - (-switches_per_cm * cm.max(0.0)).exp()).clamp(1e-9, 0.5) + }; + + // Prior: the stationary share of the archaic state, from the rates themselves rather than a + // tuned constant — with p_arch > p_bg the algebra puts it at a few percent, matching reality. + let prior_arch = ((p_bg - (1.0 - p_arch) * 0.0) / p_arch).clamp(0.001, 0.5) * 0.1; + let mut fwd = vec![[f64::NEG_INFINITY; 2]; n]; + let e0 = emit(0); + fwd[0] = [ln(1.0 - prior_arch) + e0[0], ln(prior_arch) + e0[1]]; + for i in 1..n { + let s = sw(i - 1); + let (stay, go) = (ln(1.0 - s), ln(s)); + let e = emit(i); + for st in 0..2 { + let from0 = fwd[i - 1][0] + if st == 0 { stay } else { go }; + let from1 = fwd[i - 1][1] + if st == 1 { stay } else { go }; + fwd[i][st] = ln_sum_exp(from0, from1) + e[st]; + } + } + let mut bwd = vec![[0.0f64; 2]; n]; + for i in (0..n - 1).rev() { + let s = sw(i); + let (stay, go) = (ln(1.0 - s), ln(s)); + let e = emit(i + 1); + for st in 0..2 { + let to0 = bwd[i + 1][0] + e[0] + if st == 0 { stay } else { go }; + let to1 = bwd[i + 1][1] + e[1] + if st == 1 { stay } else { go }; + bwd[i][st] = ln_sum_exp(to0, to1); + } + } + let total = ln_sum_exp(fwd[n - 1][0], fwd[n - 1][1]); + (0..n) + .map(|i| (fwd[i][1] + bwd[i][1] - total).exp().clamp(0.0, 1.0)) + .collect() +} + +/// Call archaic tracts for one subject by matching the archaic genomes. +/// +/// `observations` is per contig, already reduced by [`observations_for_contig`], so this function +/// does no I/O and no asset decoding — it is the model, and is unit-testable as such. +pub fn call_from_observations( + observations: &BTreeMap>, + gmap: &GeneticMap, + callable: &ArchaicCallable, + cfg: &MatchConfig, +) -> ArchaicSegmentResult { + let (carried, total): (usize, usize) = observations + .values() + .flatten() + .fold((0, 0), |(c, t), o| (c + usize::from(o.carries), t + 1)); + if total == 0 { + return ArchaicSegmentResult { + segments: Vec::new(), + summary: ArchaicSummary { + total_mb: 0.0, + pct_callable: 0.0, + callable_mb: 0.0, + neanderthal_mb: 0.0, + denisovan_mb: 0.0, + unknown_mb: 0.0, + n_segments: 0, + }, + }; + } + // The genome-wide rate is dominated by non-archaic sequence (archaic tracts are a few percent + // of it), so it estimates the background directly. + let p_bg = cfg + .p_background + .unwrap_or((carried as f64 / total as f64).clamp(0.001, 0.5)); + let p_arch = cfg + .p_archaic + .unwrap_or((p_bg * cfg.archaic_ratio).clamp(p_bg * 1.1, 0.95)); + + let mut segments = Vec::new(); + for (contig, obs) in observations { + if obs.len() < cfg.min_sites { + continue; + } + let post = posteriors(obs, contig, gmap, p_bg, p_arch, cfg.switches_per_cm); + let mut i = 0usize; + while i < post.len() { + if post[i] < cfg.min_posterior { + i += 1; + continue; + } + let start = i; + while i < post.len() && post[i] >= cfg.min_posterior { + i += 1; + } + let end = i - 1; + let n_sites = end - start + 1; + let span = obs[end].position - obs[start].position; + if n_sites < cfg.min_sites || span < cfg.min_segment_bp { + continue; + } + let mean_post = post[start..=end].iter().sum::() / n_sites as f64; + let (mut nea, mut den) = (0usize, 0usize); + for o in &obs[start..=end] { + if !o.carries { + continue; + } + match o.class { + DiagnosticClass::Neanderthal => nea += 1, + DiagnosticClass::Denisovan => den += 1, + DiagnosticClass::SharedArchaic => {} + } + } + segments.push(ArchaicSegment { + contig: contig.clone(), + start: obs[start].position, + end: obs[end].position, + posterior: mean_post, + n_private: obs[start..=end].iter().filter(|o| o.carries).count(), + // Attribution stays off by default; the lineage signal is a separate question this + // module does not answer (see `MatchConfig::attribute_lineage`). + source: ArchaicSource::Unknown, + neanderthal_matches: nea, + denisovan_matches: den, + }); + } + } + + let callable_mb: f64 = callable + .contigs + .iter() + .filter(|c| observations.contains_key(&c.contig)) + .flat_map(|c| c.callable_bp.iter()) + .map(|&b| b as f64) + .sum::() + / 1_000_000.0; + let total_mb: f64 = segments.iter().map(|s| s.length_mb()).sum(); + let summary = ArchaicSummary { + total_mb, + pct_callable: if callable_mb > 0.0 { total_mb * 100.0 / callable_mb } else { 0.0 }, + callable_mb, + neanderthal_mb: 0.0, + denisovan_mb: 0.0, + unknown_mb: total_mb, + n_segments: segments.len(), + }; + ArchaicSegmentResult { segments, summary } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::archaic::{CallableContig, ClassifyContig, PositionStream}; + + fn gmap(contig: &str, len: i32) -> GeneticMap { + GeneticMap::uniform(1.0, &[(contig, len)]) + } + + fn callable(contig: &str, windows: usize) -> ArchaicCallable { + ArchaicCallable { + build: "chm13v2.0".into(), + window_bp: 1_000, + contigs: vec![CallableContig { + contig: contig.into(), + start: 0, + callable_bp: vec![1_000u16; windows], + }], + } + } + + fn obs(positions: &[(i64, bool)]) -> Vec { + positions + .iter() + .map(|&(position, carries)| SiteObs { + position, + carries, + class: DiagnosticClass::SharedArchaic, + }) + .collect() + } + + /// A run of carried sites against a background of non-carried ones is what a real tract looks + /// like, and is the thing this model exists to find. + #[test] + fn finds_a_run_of_carried_sites() { + let mut sites: Vec<(i64, bool)> = (0..60).map(|i| (10_000 + i * 500, false)).collect(); + for s in sites.iter_mut().skip(20).take(20) { + s.1 = true; // a 10 kb tract, 20 diagnostic sites, all carried + } + let mut m = BTreeMap::new(); + m.insert("chr21".to_string(), obs(&sites)); + let cfg = MatchConfig { + p_background: Some(0.13), + p_archaic: Some(0.40), + min_sites: 5, + min_segment_bp: 1_000, + ..Default::default() + }; + let r = call_from_observations(&m, &gmap("chr21", 60_000), &callable("chr21", 60), &cfg); + assert_eq!(r.segments.len(), 1, "one tract expected, got {:?}", r.segments); + let seg = &r.segments[0]; + assert!(seg.start >= 19_000 && seg.start <= 21_000, "start {} off", seg.start); + assert!(seg.end >= 29_000 && seg.end <= 31_000, "end {} off", seg.end); + } + + /// The failure that gated the density caller was calling tracts out of background noise. With + /// no carried sites at all there is nothing to call, and the model must say so. + #[test] + fn calls_nothing_on_a_background_only_contig() { + let sites: Vec<(i64, bool)> = (0..200).map(|i| (10_000 + i * 500, false)).collect(); + let mut m = BTreeMap::new(); + m.insert("chr21".to_string(), obs(&sites)); + let r = call_from_observations( + &m, + &gmap("chr21", 200_000), + &callable("chr21", 200), + &MatchConfig { + p_background: Some(0.13), + p_archaic: Some(0.40), + ..Default::default() + }, + ); + assert!(r.segments.is_empty(), "background should call nothing, got {:?}", r.segments); + } + + /// Scattered carried sites at the background rate must not accumulate into a tract — the + /// density caller's defining failure, restated in this observable. + #[test] + fn scattered_background_carriers_do_not_form_a_tract() { + // 13 % carried, evenly spread: exactly the background rate, no run. + let sites: Vec<(i64, bool)> = (0..300).map(|i| (10_000 + i * 500, i % 8 == 0)).collect(); + let mut m = BTreeMap::new(); + m.insert("chr21".to_string(), obs(&sites)); + let r = call_from_observations( + &m, + &gmap("chr21", 300_000), + &callable("chr21", 300), + &MatchConfig { + p_background: Some(0.13), + p_archaic: Some(0.40), + ..Default::default() + }, + ); + assert!(r.segments.is_empty(), "background-rate carriers formed {:?}", r.segments); + } + + /// A site whose derived allele IS the reference base separates nothing, and a no-call there + /// means the opposite of what it means elsewhere. Such sites must be dropped, not counted. + #[test] + fn observations_drop_sites_where_reference_is_derived() { + let classify = ArchaicClassify { + build: "chm13v2.0".into(), + contigs: vec![ClassifyContig { + positions: PositionStream::encode("chr21", &[1_000, 2_000, 3_000]), + derived: vec![b'A', b'C', b'G'], + classes: vec![0, 1, 2], + }], + }; + let calls: BTreeMap = BTreeMap::new(); + // The reference carries the derived allele at 2_000 only. + let out = observations_for_contig( + "chr21", + &classify, + &calls, + |p| if p == 2_000 { Some(b'C') } else { Some(b'T') }, + &callable("chr21", 10), + 0.5, + ); + assert_eq!(out.len(), 2, "the reference-derived site must be dropped"); + assert!(out.iter().all(|o| o.position != 2_000)); + assert!(out.iter().all(|o| !o.carries), "no calls means nothing carried"); + } + + /// A no-call is hom-reference, i.e. NOT carrying. Conditioning on "has a call" instead is what + /// made an early version of this analysis report ~80 % carrying against a 4.3 % background. + #[test] + fn a_missing_call_is_not_a_carrier() { + let classify = ArchaicClassify { + build: "chm13v2.0".into(), + contigs: vec![ClassifyContig { + positions: PositionStream::encode("chr21", &[1_000, 2_000]), + derived: vec![b'A', b'A'], + classes: vec![0, 0], + }], + }; + let carried = SiteGenotype { + name: String::new(), + contig: "chr21".into(), + position: 1_000, + reference_allele: "T".into(), + alternate_allele: "A".into(), + ploidy: 2, + dosage: 1, + gq: 60, + depth: 30, + ref_depth: 15, + alt_depth: 15, + pls: Vec::new(), + gt: None, + allele_depths: None, + }; + let mut calls: BTreeMap = BTreeMap::new(); + calls.insert(1_000, &carried); + let out = observations_for_contig("chr21", &classify, &calls, |_| Some(b'T'), &callable("chr21", 10), 0.5); + assert_eq!(out.len(), 2); + assert!(out[0].carries, "a called derived allele carries"); + assert!(!out[1].carries, "an absent call is hom-reference, not a carrier"); + } +} diff --git a/crates/navigator-analysis/src/lib.rs b/crates/navigator-analysis/src/lib.rs index 87ea507..62e49a8 100644 --- a/crates/navigator-analysis/src/lib.rs +++ b/crates/navigator-analysis/src/lib.rs @@ -11,6 +11,7 @@ pub mod ancestry; pub mod archaic; +pub mod archaic_match; pub mod archaic_segments; pub mod caller; pub mod callset; diff --git a/scripts/archaic-validation/compare_locations.py b/scripts/archaic-validation/compare_locations.py index 98e08c1..4a25f83 100644 --- a/scripts/archaic-validation/compare_locations.py +++ b/scripts/archaic-validation/compare_locations.py @@ -17,7 +17,7 @@ from collections import defaultdict -def load_bed(path, keep=None): +def load_bed(path, keep=None, tol=0): iv = defaultdict(list) for line in open(path): p = line.rstrip("\n").split("\t") @@ -26,14 +26,14 @@ def load_bed(path, keep=None): if keep and p[0] not in keep: continue iv[p[0]].append((int(p[1]), int(p[2]))) - return {c: union(v) for c, v in iv.items()} + return {c: union(v, tol) for c, v in iv.items()} -def union(iv): +def union(iv, tol=0): iv = sorted(iv) out, cs, ce = [], *iv[0] for s, e in iv[1:]: - if s > ce: + if s > ce + tol: out.append((cs, ce)) cs, ce = s, e else: @@ -65,6 +65,30 @@ def intersect(a, b): return out +def null_sensitivity(ours, truth, truth_bp, draws=400): + """Sensitivity a caller of OUR extent would score by placing its segments at random. + + Sensitivity rises with how much sequence you call, so the raw number is meaningless without + this. The density caller scored 2.1% against a 5.0% null -- below chance -- and reporting only + the 2.1% would have looked like weak performance rather than none. + """ + import random + rnd = random.Random(0) + lens = {c: [e - s for s, e in v] for c, v in ours.items()} + vals = [] + for _ in range(draws): + r = {} + for c, L in lens.items(): + if c not in truth: + continue + lo = min(s for s, _ in truth[c]) + hi = max(e for _, e in truth[c]) + r[c] = union([(p, p + l) for l in L for p in [rnd.randint(lo, max(lo, hi - l))]]) + vals.append(intersect(r, truth) / truth_bp * 100) # both in bp + vals.sort() + return sum(vals) / len(vals), vals[int(0.95 * len(vals)) - 1], vals[-1] + + def main(): ours_path, truth_path = sys.argv[1], sys.argv[2] contigs = set(sys.argv[3:]) or None @@ -78,7 +102,9 @@ def main(): continue ours_iv[c].append((int(s["start"]), int(s["end"]))) ours = {c: union(v) for c, v in ours_iv.items()} if ours_iv else {} - truth = load_bed(truth_path, contigs) + # tol=1000: the lift splits tracts at median 2 bp gaps; a strict union reports 423 + # shards where there are 48 real tracts, which makes per-tract recovery meaningless. + truth = load_bed(truth_path, contigs, tol=1000) o_mb, t_mb = total(ours) / 1e6, total(truth) / 1e6 inter = intersect(ours, truth) / 1e6 @@ -108,9 +134,16 @@ def main(): miss += 1 print(f" their tracts hit : {hit}/{hit + miss} ({hit / (hit + miss) * 100:.1f}%)" if hit + miss else " n/a") print() - print(" A random caller with our extent would score sensitivity ~= our_Mb / callable_Mb,") - print(" i.e. a few percent. Sensitivity near that floor means we match the AMOUNT but not") - print(" the LOCATION -- which would mean the 1.01x cohort agreement was luck.") + if t_mb and ours: + mean, p95, mx = null_sensitivity(ours, truth, total(truth)) + sens = inter / t_mb * 100 + print() + print("4. AGAINST THE NULL (sensitivity alone rises with how much you call)") + print(f" random placement : mean {mean:5.1f}% p95 {p95:5.1f}% max {mx:5.1f}%") + verdict = ("ABOVE the null's full range" if sens > mx + else "above p95" if sens > p95 + else "AT OR BELOW CHANCE") + print(f" observed {sens:5.1f}% : {verdict}") if __name__ == "__main__": From e9957bd0f33d596149eff02222fe7e2837de1a19 Mon Sep 17 00:00:00 2001 From: James Kane Date: Fri, 31 Jul 2026 09:05:51 -0500 Subject: [PATCH 3/9] Archaic Tier B v2: calibrate on 30 Europeans, report on 30 held-out ones Thresholds fitted by grid search on a TRAIN split and reported on a TEST split that never touched the fit. The split exists because the density caller was tuned until a cohort statistic matched and the statistic was then reported as evidence; this makes that failure visible instead of absorbing it. HELD-OUT TEST, 30 Europeans, chr21+22, against hmmix's own calls: density defaults calibrated base-level F1 -- 27.9% 32.5% precision 1.5% 20.2% 28.8% extent ratio ours/theirs 1.45 2.23 1.30 per-individual extent r -0.018 +0.520 +0.642 (p=0.94) (p=0.0001) Locations were already decisive on the earlier 20-sample cohort: 20/20 individuals above their OWN random-placement null, computed at their own called extent, mean sensitivity 45.3% against a 7.1% null. The density caller scored 2.1% against 5.0% -- below chance. Objective was base-level F1, not sensitivity: sensitivity is bought by calling more sequence, and the uncalibrated caller over-called 2.2x while still scoring 45%. F1 makes over-calling cost something. THE ARGMAX WAS NOT TAKEN, deliberately. The grid's best TRAIN F1 wanted min_segment_bp = 40 kb, which discards 61% of real tracts by construction -- hmmix's median tract is 31-36 kb and its p10 is 7 kb. It bought 0.1 F1 points over the 5 kb floor and was WORSE on both things the number exists for: sensitivity (37.4% vs 38.2%) and per-individual extent correlation (+0.642 vs +0.658). The design doc already records this trap once, at 50 kb. Chosen: min_posterior 0.95, min_sites 24, min_segment_bp 5 kb. One reporting bug fixed while writing this: the harness flagged ANY train-test gap over 5 points as overfitting, including the case where TEST scores higher. Test above train is split-to-split variation at n=30, not the grid finding the split. Overfitting is train >> test, and only that now trips the warning. STILL NOT GOOD ENOUGH TO RE-ENABLE, and the module says so: precision is 28.8% and the extent ratio 1.30, so the caller over-calls by about a third. ARCHAIC_SEGMENTS_ENABLED stays false until that closes and the result reproduces outside Europe -- the cohort here is European only, on two chromosomes, and `archaic_ratio` was held at its measured value rather than swept because it changes the emissions and needs a re-run per value. Banks calibrate.py (the split + grid) and cohort_score.py (per-individual scoring against each individual's own null) alongside the existing harness. Co-Authored-By: Claude Opus 5 (1M context) --- .../navigator-analysis/src/archaic_match.rs | 35 ++- scripts/archaic-validation/calibrate.py | 200 ++++++++++++++++++ scripts/archaic-validation/cohort_score.py | 166 +++++++++++++++ 3 files changed, 399 insertions(+), 2 deletions(-) create mode 100644 scripts/archaic-validation/calibrate.py create mode 100644 scripts/archaic-validation/cohort_score.py diff --git a/crates/navigator-analysis/src/archaic_match.rs b/crates/navigator-analysis/src/archaic_match.rs index 8428fc5..72ce308 100644 --- a/crates/navigator-analysis/src/archaic_match.rs +++ b/crates/navigator-analysis/src/archaic_match.rs @@ -35,6 +35,28 @@ //! //! Transitions stay recombination-scaled between consecutive sites, as in [`crate::roh`] and the //! chromosome painter. +//! +//! # Validation +//! +//! Scored against hmmix's own calls for the same individuals, 60 Europeans on chr21+22, split 30 +//! **train** / 30 **test** on a fixed seed. Thresholds were fitted on train only; every figure below +//! is the held-out half. The split exists because the previous caller was tuned until a cohort +//! statistic matched and the statistic was then reported as evidence. +//! +//! | | density caller | this, defaults | this, calibrated | +//! |---|---|---|---| +//! | base-level F1 | — | 27.9 % | **32.5 %** | +//! | precision | 1.5 % | 20.2 % | **28.8 %** | +//! | extent ratio ours/theirs | 1.45 | 2.23 | **1.30** | +//! | per-individual extent `r` | −0.018 (p = 0.94) | +0.520 | **+0.642 (p = 0.0001)** | +//! +//! On locations, all 20 individuals of an earlier cohort scored above their own random-placement +//! null (mean 45.3 % sensitivity against a 7.1 % null); the density caller scored 2.1 % against a +//! 5.0 % null, i.e. below chance. +//! +//! **Not yet good enough to re-enable.** Precision is 28.8 % and the extent ratio 1.30, so the +//! caller still over-calls by about a third. `ARCHAIC_SEGMENTS_ENABLED` stays `false` until that is +//! closed and the result reproduces outside Europe. use std::collections::BTreeMap; @@ -93,10 +115,19 @@ impl Default for MatchConfig { MatchConfig { p_background: None, p_archaic: None, + // Measured (39.5 % carrying inside real tracts against 13.0 % elsewhere), not fitted. archaic_ratio: 3.04, switches_per_cm: 1.0, - min_posterior: 0.80, - min_sites: 8, + // CALIBRATED on 30 Europeans and reported on 30 held-out ones (see the module docs). + // Fitted by grid search on base-level F1, because sensitivity alone is bought by calling + // more sequence — the uncalibrated caller over-called 2.2x and still scored 45 %. + min_posterior: 0.95, + min_sites: 24, + // 5 kb, NOT the 40 kb the grid's argmax preferred. That floor scored 0.1 F1 points + // higher while discarding 61 % of real tracts by construction — hmmix's median tract is + // 31-36 kb and its p10 is 7 kb — and it scored *worse* on the two things the number is + // for: sensitivity (37.4 % vs 38.2 %) and per-individual extent correlation (+0.642 vs + // +0.658). The design already recorded this exact trap once, at 50 kb. min_segment_bp: 5_000, min_callable_fraction: 0.5, attribute_lineage: false, diff --git a/scripts/archaic-validation/calibrate.py b/scripts/archaic-validation/calibrate.py new file mode 100644 index 0000000..8768e7d --- /dev/null +++ b/scripts/archaic-validation/calibrate.py @@ -0,0 +1,200 @@ +"""Calibrate the reference-based archaic caller on a TRAIN split, report on a held-out TEST split. + +The density caller looked validated because it was tuned until a cohort statistic matched, and the +statistic was then reported as evidence. The split exists so that cannot happen again: every number +quoted as performance comes from individuals whose data never touched the fit. + +Objective is base-level F1 against the external callset. F1 rather than sensitivity because +sensitivity alone is bought by calling more sequence -- the current caller over-calls 2.2x and still +scores 45% sensitivity. F1 makes over-calling cost something. + +Reading the output: TRAIN F1 is the fit, TEST F1 is the claim. A large gap between them means the +grid found the split rather than the signal. +""" + +import collections +import glob +import json +import math +import os +import random +import subprocess +import sys + +SC = os.path.dirname(os.path.abspath(__file__)) +V2 = os.path.join(SC, 'v2') +CONTIGS = ('chr21', 'chr22') + + +def union(iv, tol=0): + iv = sorted(iv) + out, cs, ce = [], *iv[0] + for s, e in iv[1:]: + if s > ce + tol: + out.append((cs, ce)) + cs, ce = s, e + else: + ce = max(ce, e) + out.append((cs, ce)) + return out + + +def load_bed(path, tol=0): + d = collections.defaultdict(list) + for line in open(path): + p = line.split() + if len(p) >= 3: + d[p[0]].append((int(p[1]), int(p[2]))) + return {c: union(v, tol) for c, v in d.items()} + + +def bp(d): + return sum(e - s for v in d.values() for s, e in v) + + +def intersect(a, b): + tot = 0 + for c, av in a.items(): + for s, e in av: + for s2, e2 in b.get(c, []): + lo, hi = max(s, s2), min(e, e2) + if hi > lo: + tot += hi - lo + return tot + + +def pearson(x, y): + n = len(x) + if n < 3: + return 0.0 + mx, my = sum(x) / n, sum(y) / n + num = sum((a - mx) * (b - my) for a, b in zip(x, y)) + dx = math.sqrt(sum((a - mx) ** 2 for a in x)) + dy = math.sqrt(sum((b - my) ** 2 for b in y)) + return num / (dx * dy) if dx and dy else 0.0 + + +def segments_from(doc, min_post, min_sites, min_bp): + """Re-filter a caller run at new thresholds. + + The HMM posterior is fixed per run, so the three post-hoc thresholds can be swept without + re-running the caller. `archaic_ratio` cannot -- it changes the emissions -- so it is swept by + re-running the probe, outside this function. + """ + oi = collections.defaultdict(list) + for s in doc['segments']: + if s['contig'] not in CONTIGS: + continue + if s['posterior'] < min_post or s['n_private'] < min_sites: + continue + if s['end'] - s['start'] < min_bp: + continue + oi[s['contig']].append((s['start'], s['end'])) + return {c: union(v) for c, v in oi.items()} if oi else {} + + +def score(samples, runs, truths, min_post, min_sites, min_bp): + sens_n = sens_d = prec_d = 0 + t_mb, o_mb = [], [] + for s in samples: + ours = segments_from(runs[s], min_post, min_sites, min_bp) + truth = truths[s] + tb, ob = bp(truth), bp(ours) + ov = intersect(ours, truth) if ours else 0 + sens_n += ov + sens_d += tb + prec_d += ob + t_mb.append(tb / 1e6) + o_mb.append(ob / 1e6) + sens = sens_n / sens_d if sens_d else 0.0 + prec = sens_n / prec_d if prec_d else 0.0 + f1 = 2 * sens * prec / (sens + prec) if (sens + prec) else 0.0 + return { + 'sens': sens * 100, 'prec': prec * 100, 'f1': f1 * 100, + 'ratio': (sum(o_mb) / sum(t_mb)) if sum(t_mb) else 0.0, + 'r': pearson(t_mb, o_mb), + } + + +def main(): + runs, truths = {}, {} + for f in sorted(glob.glob(f'{V2}/*.match.json')): + s = os.path.basename(f)[:-len('.match.json')] + t = f'{V2}/truth_{s}.bed' + if os.path.exists(t): + runs[s] = json.load(open(f)) + truths[s] = load_bed(t, tol=1000) + samples = sorted(runs) + if len(samples) < 20: + print(f'only {len(samples)} scored runs — need the cohort first') + return + + rnd = random.Random(20260731) + shuffled = samples[:] + rnd.shuffle(shuffled) + half = len(shuffled) // 2 + train, test = sorted(shuffled[:half]), sorted(shuffled[half:]) + print(f'n = {len(samples)} train {len(train)} test {len(test)} (seed fixed)\n') + + grid = [] + for mp in (0.80, 0.85, 0.90, 0.95, 0.98): + for ms in (8, 12, 16, 24, 32): + for mb in (5_000, 10_000, 20_000, 40_000): + grid.append((mp, ms, mb)) + + scored = [(score(train, runs, truths, *g), g) for g in grid] + scored.sort(key=lambda x: -x[0]['f1']) + + print('TOP 5 ON TRAIN (fitted)') + print(f" {'post':>5s} {'sites':>5s} {'minbp':>7s} {'F1':>6s} {'sens':>6s} {'prec':>6s} " + f"{'ratio':>6s} {'r':>6s}") + for sc, g in scored[:5]: + print(f' {g[0]:5.2f} {g[1]:5d} {g[2]:7d} {sc["f1"]:5.1f}% {sc["sens"]:5.1f}% ' + f'{sc["prec"]:5.1f}% {sc["ratio"]:6.2f} {sc["r"]:+6.3f}') + + best = scored[0][1] + base = (0.80, 8, 5_000) + print(f'\nCHOSEN (best TRAIN F1): min_posterior {best[0]}, min_sites {best[1]}, ' + f'min_segment_bp {best[2]}') + print('\n F1 sens prec ratio r') + for label, samp in (('TRAIN (fitted)', train), ('TEST (held out)', test)): + b = score(samp, runs, truths, *base) + c = score(samp, runs, truths, *best) + print(f' {label:22s}') + print(f' before (defaults) {b["f1"]:5.1f}% {b["sens"]:6.1f}% {b["prec"]:6.1f}% ' + f'{b["ratio"]:6.2f} {b["r"]:+7.3f}') + print(f' after (calibrated){c["f1"]:5.1f}% {c["sens"]:6.1f}% {c["prec"]:6.1f}% ' + f'{c["ratio"]:6.2f} {c["r"]:+7.3f}') + tr = score(train, runs, truths, *best)['f1'] + te = score(test, runs, truths, *best)['f1'] + gap = tr - te + # Overfitting is train >> test. Test scoring HIGHER is split-to-split variation, not the grid + # having found the split -- reading any large gap as overfitting would be its own error. + if gap > 5: + verdict = 'OVERFIT — the grid found the split, not the signal' + elif gap < -5: + verdict = 'test split is easier; not overfitting, but the gap is split variation at n=30' + else: + verdict = 'generalises' + print(f'\n train F1 {tr:.1f}% test F1 {te:.1f}% gap {gap:+.1f} points ({verdict})') + + # Is the held-out extent correlation real, or n=30 noise? + t_mb, o_mb = [], [] + for s in test: + ours = segments_from(runs[s], *best) + t_mb.append(bp(truths[s]) / 1e6) + o_mb.append(bp(ours) / 1e6 if ours else 0.0) + r = pearson(t_mb, o_mb) + rnd = random.Random(0) + yy = list(o_mb) + hits = 0 + for _ in range(20000): + rnd.shuffle(yy) + if abs(pearson(t_mb, yy)) >= abs(r): + hits += 1 + print(f' held-out extent correlation r = {r:+.3f}, permutation p = {(hits + 1) / 20001:.4f}' + f' (density caller: r = -0.018, p = 0.94)') + + +if __name__ == '__main__': + main() diff --git a/scripts/archaic-validation/cohort_score.py b/scripts/archaic-validation/cohort_score.py new file mode 100644 index 0000000..9a2f8f3 --- /dev/null +++ b/scripts/archaic-validation/cohort_score.py @@ -0,0 +1,166 @@ +"""Score the reference-based caller across the cohort, on the tests the density caller failed. + +Two questions, both of which the density caller answered badly: + + LOCATIONS -- is base-level overlap above the random-placement null, per individual? + (density: 2.1% against a 5.0% null, i.e. below chance) + AMOUNTS -- does extent track the individual across people? + (density: Pearson r = -0.018, p = 0.94, against a 2.5x range of true values) + +Reports both, plus precision, because over-calling inflates sensitivity and the pair has to be +read together. +""" + +import collections +import glob +import json +import math +import os +import random + +SC = os.path.dirname(os.path.abspath(__file__)) +V2 = os.path.join(SC, 'v2') +CONTIGS = ('chr21', 'chr22') + + +def union(iv, tol=0): + iv = sorted(iv) + out, cs, ce = [], *iv[0] + for s, e in iv[1:]: + if s > ce + tol: + out.append((cs, ce)) + cs, ce = s, e + else: + ce = max(ce, e) + out.append((cs, ce)) + return out + + +def load_bed(path, tol=0): + d = collections.defaultdict(list) + for line in open(path): + p = line.split() + if len(p) < 3: + continue + d[p[0]].append((int(p[1]), int(p[2]))) + return {c: union(v, tol) for c, v in d.items()} + + +def bp(d): + return sum(e - s for v in d.values() for s, e in v) + + +def intersect(a, b): + tot = 0 + for c, av in a.items(): + for s, e in av: + for s2, e2 in b.get(c, []): + lo, hi = max(s, s2), min(e, e2) + if hi > lo: + tot += hi - lo + return tot + + +def null_mean_p95(ours, truth, tbp, draws=200, seed=0): + rnd = random.Random(seed) + lens = {c: [e - s for s, e in v] for c, v in ours.items()} + vals = [] + for _ in range(draws): + r = {} + for c, L in lens.items(): + if c not in truth: + continue + lo = min(s for s, _ in truth[c]) + hi = max(e for _, e in truth[c]) + r[c] = union([(p, p + l) for l in L for p in [rnd.randint(lo, max(lo, hi - l))]]) + vals.append(intersect(r, truth) / tbp * 100) + vals.sort() + return sum(vals) / len(vals), vals[int(0.95 * len(vals)) - 1] + + +def pearson(x, y): + n = len(x) + mx, my = sum(x) / n, sum(y) / n + num = sum((a - mx) * (b - my) for a, b in zip(x, y)) + dx = math.sqrt(sum((a - mx) ** 2 for a in x)) + dy = math.sqrt(sum((b - my) ** 2 for b in y)) + return num / (dx * dy) if dx and dy else 0.0 + + +def rank(v): + order = sorted(range(len(v)), key=lambda i: v[i]) + r = [0.0] * len(v) + for pos, i in enumerate(order): + r[i] = pos + 1 + return r + + +def perm_p(x, y, draws=20000): + obs = abs(pearson(x, y)) + rnd = random.Random(0) + yy = list(y) + hits = 0 + for _ in range(draws): + rnd.shuffle(yy) + if abs(pearson(x, yy)) >= obs: + hits += 1 + return (hits + 1) / (draws + 1) + + +def main(): + rows = [] + for f in sorted(glob.glob(f'{V2}/*.match.json')): + s = os.path.basename(f)[:-len('.match.json')] + tpath = f'{V2}/truth_{s}.bed' + if not os.path.exists(tpath): + continue + doc = json.load(open(f)) + oi = collections.defaultdict(list) + for seg in doc['segments']: + if seg['contig'] in CONTIGS: + oi[seg['contig']].append((seg['start'], seg['end'])) + if not oi: + continue + ours = {c: union(v) for c, v in oi.items()} + truth = load_bed(tpath, tol=1000) + tbp, obp = bp(truth), bp(ours) + ov = intersect(ours, truth) + mean, p95 = null_mean_p95(ours, truth, tbp) + rows.append({ + 'sample': s, 'truth_mb': tbp / 1e6, 'ours_mb': obp / 1e6, + 'sens': ov / tbp * 100, 'prec': ov / obp * 100, + 'null': mean, 'p95': p95, 'nseg': sum(len(v) for v in ours.values()), + }) + + rows.sort(key=lambda r: r['truth_mb']) + print(f'n = {len(rows)} Europeans, chr21+22, reference-based caller\n') + print(f"{'sample':10s} {'truth':>7s} {'ours':>7s} {'segs':>5s} {'sens':>7s} {'prec':>6s} " + f"{'null':>6s} {'p95':>6s} verdict") + beat = 0 + for r in rows: + ok = r['sens'] > r['p95'] + beat += ok + print(f"{r['sample']:10s} {r['truth_mb']:7.3f} {r['ours_mb']:7.3f} {r['nseg']:5d} " + f"{r['sens']:6.1f}% {r['prec']:5.1f}% {r['null']:5.1f}% {r['p95']:5.1f}% " + f"{'above null' if ok else 'AT/BELOW'}") + + print(f'\nLOCATIONS: {beat}/{len(rows)} individuals score above their own p95 null') + print(f' mean sensitivity {sum(r["sens"] for r in rows) / len(rows):.1f}% ' + f'mean null {sum(r["null"] for r in rows) / len(rows):.1f}% ' + f'mean precision {sum(r["prec"] for r in rows) / len(rows):.1f}%') + print(' (density caller: 2.1% against a 5.0% null -- below chance)') + + t = [r['truth_mb'] for r in rows] + o = [r['ours_mb'] for r in rows] + r_ = pearson(t, o) + rho = pearson(rank(t), rank(o)) + print(f'\nAMOUNTS: per-individual extent correlation') + print(f' Pearson r = {r_:+.3f} permutation p = {perm_p(t, o):.4f}') + print(f' Spearman rho = {rho:+.3f} permutation p = {perm_p(rank(t), rank(o)):.4f}') + print(f' truth range {min(t):.2f}-{max(t):.2f} Mb ours {min(o):.2f}-{max(o):.2f} Mb') + print(f' mean ratio ours/theirs = {sum(o) / sum(t):.3f}') + print(' (density caller: r = -0.018, p = 0.94)') + + +if __name__ == '__main__': + main() From d279a5c99070a7a97eb49fb97ccb74847073cfa8 Mon Sep 17 00:00:00 2001 From: James Kane Date: Fri, 31 Jul 2026 09:31:22 -0500 Subject: [PATCH 4/9] Archaic Tier B v2: sweep the emission ratio; over-calling is gone The previous calibration swept only the three post-hoc thresholds, because archaic_ratio changes the emissions and cannot be re-filtered from a finished run. Sweeping it is what removed the systematic over-calling. The probe now sweeps it IN PROCESS (`ARCHAIC_RATIOS=2.0,2.5,...`): the expensive work -- reading the reference, walking the diagnostic sites -- is per sample, not per ratio, so eight values cost one pass rather than eight. 60 samples x 8 ratios runs in well under a second each. HELD-OUT TEST, 30 Europeans that never touched the fit: density uncalibrated thresholds + ratio base-level F1 -- 27.9% 33.1% 34.5% precision 1.5% 20.2% 31.4% 34.9% extent ratio ours/theirs 1.45 2.23 1.11 0.98 per-individual extent r -0.018 +0.520 +0.628 +0.710 (p=0.94) (p<0.0001) The extent ratio is the line that matters: 0.98 means the caller is no longer systematically over-calling, which no threshold could fix because the emissions were wrong. archaic_ratio is now FITTED, not measured, and the docs say so. The observed enrichment inside real tracts is 3.04x but the model separates best at 4.5x -- not a contradiction: 3.04x is the average over an external tract set that is itself weakly supported, while the emission ratio is what makes the HMM selective enough to place boundaries. p_background remains estimated per individual from the subject's own data, unfitted. THE ARGMAX WAS AGAIN NOT TAKEN. The grid preferred min_segment_bp = 10 kb; 5 kb is chosen. Within the plateau the two differ by 0.1 F1, 5 kb is slightly BETTER on per-individual correlation (+0.710 vs +0.706), and it discards half as many real tracts (8% of the truth under 5 kb against 16% under 10 kb). An earlier sweep wanted 40 kb, which discards 61%; the design records the same trap at 50 kb. Structural exclusion of real tracts is not worth a tenth of a point. One unit test needed pinning rather than fixing: `finds_a_run_of_carried_sites` broke when min_posterior rose to 0.98 and trimmed the run's edges. That is correct behaviour and the wrong thing for a model test to be sensitive to, so it now sets its own thresholds instead of inheriting the calibrated ones. STILL GATED, with specific limits rather than general unease: precision 34.9% means two thirds of called sequence is not in the reference callset; the cohort is European only and chr21+22 only; and the reference callset is itself weakly supported (hmmix's tracts are enriched just 1.84x for their own archaic SNPs), so agreement with it caps well below 100% even for a correct caller. East Asians are the sharp next test -- the truth predicts ~1.18x more archaic sequence there, which a caller merely tracking European structure would miss. Co-Authored-By: Claude Opus 5 (1M context) --- .../examples/archaic_match_probe.rs | 23 +++ .../navigator-analysis/src/archaic_match.rs | 58 +++--- scripts/archaic-validation/calibrate_full.py | 170 ++++++++++++++++++ 3 files changed, 231 insertions(+), 20 deletions(-) create mode 100644 scripts/archaic-validation/calibrate_full.py diff --git a/crates/navigator-analysis/examples/archaic_match_probe.rs b/crates/navigator-analysis/examples/archaic_match_probe.rs index 84944ee..45c1213 100644 --- a/crates/navigator-analysis/examples/archaic_match_probe.rs +++ b/crates/navigator-analysis/examples/archaic_match_probe.rs @@ -75,6 +75,29 @@ fn main() -> Result<(), Box> { GeneticMap::from_bytes(&std::fs::read(&a[3])?).map_err(|e| e.to_string())? }; + // `ARCHAIC_RATIOS=2.0,2.5,3.04` sweeps the emission ratio in one process. It cannot be swept + // post-hoc like the three thresholds — it changes the emissions, so the HMM must be re-decoded — + // but the expensive part (reading the reference, walking the diagnostic sites) is per sample, + // not per ratio, so doing it here costs one pass instead of one per value. + if let Ok(spec) = std::env::var("ARCHAIC_RATIOS") { + let mut out = serde_json::Map::new(); + for tok in spec.split(',').filter(|t| !t.trim().is_empty()) { + let ratio: f64 = tok.trim().parse()?; + let cfg = MatchConfig { + archaic_ratio: ratio, + ..Default::default() + }; + let r = call_from_observations(&observations, &gmap, &callable, &cfg); + eprintln!( + " ratio {ratio:>5.2} -> {} segments, {:.3} Mb", + r.summary.n_segments, r.summary.total_mb + ); + out.insert(tok.trim().to_string(), serde_json::to_value(&r)?); + } + println!("{}", serde_json::to_string(&out)?); + return Ok(()); + } + let result = call_from_observations(&observations, &gmap, &callable, &MatchConfig::default()); let s = &result.summary; eprintln!( diff --git a/crates/navigator-analysis/src/archaic_match.rs b/crates/navigator-analysis/src/archaic_match.rs index 72ce308..8856a91 100644 --- a/crates/navigator-analysis/src/archaic_match.rs +++ b/crates/navigator-analysis/src/archaic_match.rs @@ -43,20 +43,27 @@ //! is the held-out half. The split exists because the previous caller was tuned until a cohort //! statistic matched and the statistic was then reported as evidence. //! -//! | | density caller | this, defaults | this, calibrated | +//! | | density caller | this, uncalibrated | this, calibrated | //! |---|---|---|---| -//! | base-level F1 | — | 27.9 % | **32.5 %** | -//! | precision | 1.5 % | 20.2 % | **28.8 %** | -//! | extent ratio ours/theirs | 1.45 | 2.23 | **1.30** | -//! | per-individual extent `r` | −0.018 (p = 0.94) | +0.520 | **+0.642 (p = 0.0001)** | +//! | base-level F1 | — | 27.9 % | **34.5 %** | +//! | precision | 1.5 % | 20.2 % | **34.9 %** | +//! | extent ratio ours/theirs | 1.45 | 2.23 | **0.98** | +//! | per-individual extent `r` | −0.018 (p = 0.94) | +0.520 | **+0.710 (p < 0.0001)** | +//! +//! The extent ratio of 0.98 is the one to notice: the caller is no longer systematically +//! over-calling, which the emission-ratio sweep is what fixed. //! //! On locations, all 20 individuals of an earlier cohort scored above their own random-placement //! null (mean 45.3 % sensitivity against a 7.1 % null); the density caller scored 2.1 % against a //! 5.0 % null, i.e. below chance. //! -//! **Not yet good enough to re-enable.** Precision is 28.8 % and the extent ratio 1.30, so the -//! caller still over-calls by about a third. `ARCHAIC_SEGMENTS_ENABLED` stays `false` until that is -//! closed and the result reproduces outside Europe. +//! **Still not enough to re-enable**, and the limits are specific rather than general unease: +//! precision is 34.9 %, so two thirds of called sequence is not in the reference callset; the cohort +//! is **European only** and **chr21+22 only**; and the reference callset is itself weakly supported +//! (hmmix's own tracts are enriched only 1.84x for their own archaic SNPs), so agreement with it +//! caps out well below 100 % even for a correct caller. `ARCHAIC_SEGMENTS_ENABLED` stays `false` +//! until this reproduces outside Europe — East Asians are the sharp test, since the truth predicts +//! ~1.18x more archaic sequence there and a caller merely tracking European structure would miss it. use std::collections::BTreeMap; @@ -115,19 +122,25 @@ impl Default for MatchConfig { MatchConfig { p_background: None, p_archaic: None, - // Measured (39.5 % carrying inside real tracts against 13.0 % elsewhere), not fitted. - archaic_ratio: 3.04, + // FITTED (not measured): the observed enrichment inside real tracts is 3.04x, but the + // model separates best at 4.5x. That is not a contradiction — 3.04x is the *average* + // over an external tract set that is itself only weakly supported, while the emission + // ratio is what makes the HMM selective enough to place boundaries. Fitted on 30 + // Europeans, reported on 30 held-out ones; it is the parameter that removed the + // over-calling (extent ratio 2.23 -> 0.98). + archaic_ratio: 4.5, switches_per_cm: 1.0, - // CALIBRATED on 30 Europeans and reported on 30 held-out ones (see the module docs). - // Fitted by grid search on base-level F1, because sensitivity alone is bought by calling - // more sequence — the uncalibrated caller over-called 2.2x and still scored 45 %. - min_posterior: 0.95, - min_sites: 24, - // 5 kb, NOT the 40 kb the grid's argmax preferred. That floor scored 0.1 F1 points - // higher while discarding 61 % of real tracts by construction — hmmix's median tract is - // 31-36 kb and its p10 is 7 kb — and it scored *worse* on the two things the number is - // for: sensitivity (37.4 % vs 38.2 %) and per-individual extent correlation (+0.642 vs - // +0.658). The design already recorded this exact trap once, at 50 kb. + // All three CALIBRATED on train, reported on held-out test (see the module docs). + // Objective was base-level F1: sensitivity alone is bought by calling more sequence, and + // the uncalibrated caller over-called 2.2x while still scoring 45 %. + min_posterior: 0.98, + min_sites: 16, + // 5 kb, though the grid's argmax preferred 10 kb. Within the plateau the two differ by + // 0.1 F1 points, the 5 kb floor is slightly BETTER on per-individual extent correlation + // (+0.710 vs +0.706), and it discards half as many real tracts (8 % of the truth under + // 5 kb against 16 % under 10 kb). An earlier sweep wanted 40 kb, which would have + // discarded 61 %; the design records the same trap once before at 50 kb. Structural + // exclusion of real tracts is not worth a tenth of a point. min_segment_bp: 5_000, min_callable_fraction: 0.5, attribute_lineage: false, @@ -406,9 +419,14 @@ mod tests { } let mut m = BTreeMap::new(); m.insert("chr21".to_string(), obs(&sites)); + // Thresholds pinned rather than inherited: this test is about whether the model finds a + // run at all, and should not move when the calibrated defaults do. (It broke once when + // `min_posterior` rose to 0.98 and trimmed the run's edges — correct behaviour, wrong + // thing for this test to be sensitive to.) let cfg = MatchConfig { p_background: Some(0.13), p_archaic: Some(0.40), + min_posterior: 0.80, min_sites: 5, min_segment_bp: 1_000, ..Default::default() diff --git a/scripts/archaic-validation/calibrate_full.py b/scripts/archaic-validation/calibrate_full.py new file mode 100644 index 0000000..dbdf757 --- /dev/null +++ b/scripts/archaic-validation/calibrate_full.py @@ -0,0 +1,170 @@ +"""Full calibration: sweep the emission ratio AND the three post-hoc thresholds, train/test split. + +`archaic_ratio` is the one parameter that cannot be re-filtered from a finished run -- it changes +the emissions, so the HMM has to be re-decoded. The probe sweeps it in-process (one reference read +per sample, not per value) and writes `{ratio: result}`; this searches the joint grid. + +Fitted on TRAIN only. Every number reported as performance comes from TEST. +""" + +import collections +import glob +import json +import math +import os +import random + +SC = os.path.dirname(os.path.abspath(__file__)) +SWEEP = os.path.join(SC, 'v2', 'sweep') +V2 = os.path.join(SC, 'v2') +CONTIGS = ('chr21', 'chr22') + + +def union(iv, tol=0): + iv = sorted(iv) + out, cs, ce = [], *iv[0] + for s, e in iv[1:]: + if s > ce + tol: + out.append((cs, ce)) + cs, ce = s, e + else: + ce = max(ce, e) + out.append((cs, ce)) + return out + + +def load_bed(path, tol=0): + d = collections.defaultdict(list) + for line in open(path): + p = line.split() + if len(p) >= 3: + d[p[0]].append((int(p[1]), int(p[2]))) + return {c: union(v, tol) for c, v in d.items()} + + +def bp(d): + return sum(e - s for v in d.values() for s, e in v) + + +def intersect(a, b): + tot = 0 + for c, av in a.items(): + for s, e in av: + for s2, e2 in b.get(c, []): + lo, hi = max(s, s2), min(e, e2) + if hi > lo: + tot += hi - lo + return tot + + +def pearson(x, y): + n = len(x) + if n < 3: + return 0.0 + mx, my = sum(x) / n, sum(y) / n + num = sum((a - mx) * (b - my) for a, b in zip(x, y)) + dx = math.sqrt(sum((a - mx) ** 2 for a in x)) + dy = math.sqrt(sum((b - my) ** 2 for b in y)) + return num / (dx * dy) if dx and dy else 0.0 + + +def segs(doc, mp, ms, mb): + oi = collections.defaultdict(list) + for s in doc['segments']: + if s['contig'] not in CONTIGS or s['posterior'] < mp or s['n_private'] < ms: + continue + if s['end'] - s['start'] < mb: + continue + oi[s['contig']].append((s['start'], s['end'])) + return {c: union(v) for c, v in oi.items()} if oi else {} + + +def score(samples, sweeps, truths, ratio, mp, ms, mb): + n = d_t = d_o = 0 + t_mb, o_mb = [], [] + for s in samples: + doc = sweeps[s].get(ratio) + if doc is None: + continue + ours = segs(doc, mp, ms, mb) + tb, ob = bp(truths[s]), (bp(ours) if ours else 0) + n += intersect(ours, truths[s]) if ours else 0 + d_t += tb + d_o += ob + t_mb.append(tb / 1e6) + o_mb.append(ob / 1e6) + sens = n / d_t if d_t else 0.0 + prec = n / d_o if d_o else 0.0 + f1 = 2 * sens * prec / (sens + prec) if (sens + prec) else 0.0 + return {'f1': f1 * 100, 'sens': sens * 100, 'prec': prec * 100, + 'ratio_mb': (sum(o_mb) / sum(t_mb)) if sum(t_mb) else 0.0, + 'r': pearson(t_mb, o_mb), 't': t_mb, 'o': o_mb} + + +def perm_p(x, y, draws=20000): + obs = abs(pearson(x, y)) + rnd = random.Random(0) + yy = list(y) + hits = sum(1 for _ in range(draws) if (rnd.shuffle(yy), abs(pearson(x, yy)) >= obs)[1]) + return (hits + 1) / (draws + 1) + + +def main(): + sweeps, truths = {}, {} + for f in sorted(glob.glob(f'{SWEEP}/*.json')): + s = os.path.basename(f)[:-5] + t = f'{V2}/truth_{s}.bed' + if os.path.exists(t): + sweeps[s] = json.load(open(f)) + truths[s] = load_bed(t, tol=1000) + samples = sorted(sweeps) + if len(samples) < 40: + print(f'only {len(samples)} swept runs so far') + return + ratios = sorted(next(iter(sweeps.values())).keys(), key=float) + + rnd = random.Random(20260731) + sh = samples[:] + rnd.shuffle(sh) + h = len(sh) // 2 + train, test = sorted(sh[:h]), sorted(sh[h:]) + print(f'n = {len(samples)} train {len(train)} test {len(test)} ratios {ratios}\n') + + grid = [(r, mp, ms, mb) + for r in ratios + for mp in (0.90, 0.95, 0.98) + for ms in (8, 16, 24, 32) + for mb in (5_000, 10_000)] + scored = sorted(((score(train, sweeps, truths, *g), g) for g in grid), key=lambda x: -x[0]['f1']) + + print('TOP 8 ON TRAIN (fitted)') + print(f" {'ratio':>5s} {'post':>5s} {'sites':>5s} {'minbp':>6s} {'F1':>6s} {'sens':>6s} " + f"{'prec':>6s} {'ext':>5s} {'r':>7s}") + for sc, g in scored[:8]: + print(f' {g[0]:>5s} {g[1]:5.2f} {g[2]:5d} {g[3]:6d} {sc["f1"]:5.1f}% {sc["sens"]:5.1f}% ' + f'{sc["prec"]:5.1f}% {sc["ratio_mb"]:5.2f} {sc["r"]:+7.3f}') + + best = scored[0][1] + print(f'\nCHOSEN (best TRAIN F1): archaic_ratio {best[0]}, min_posterior {best[1]}, ' + f'min_sites {best[2]}, min_segment_bp {best[3]}') + prev = ('3.04', 0.95, 24, 5_000) + print('\n F1 sens prec ext r') + for label, samp in (('TRAIN (fitted)', train), ('TEST (held out)', test)): + print(f' {label}') + for tag, g in (('previous (ratio 3.04)', prev), ('new (swept) ', best)): + sc = score(samp, sweeps, truths, *g) + print(f' {tag} {sc["f1"]:5.1f}% {sc["sens"]:6.1f}% {sc["prec"]:6.1f}% ' + f'{sc["ratio_mb"]:6.2f} {sc["r"]:+8.3f}') + tr = score(train, sweeps, truths, *best)['f1'] + te_s = score(test, sweeps, truths, *best) + gap = tr - te_s['f1'] + verdict = ('OVERFIT — the grid found the split' if gap > 5 + else 'test split is easier; split variation, not overfitting' if gap < -5 + else 'generalises') + print(f'\n train F1 {tr:.1f}% test F1 {te_s["f1"]:.1f}% gap {gap:+.1f} ({verdict})') + print(f' held-out extent r = {te_s["r"]:+.3f}, permutation p = ' + f'{perm_p(te_s["t"], te_s["o"]):.4f}') + + +if __name__ == '__main__': + main() From e69f913416ddd639385152fe710c8c2c51f91c33 Mon Sep 17 00:00:00 2001 From: James Kane Date: Fri, 31 Jul 2026 10:46:59 -0500 Subject: [PATCH 5/9] Archaic Tier B v2: detection transfers to East Asians; the reported number does not Ran 30 East Asians with the parameters FROZEN at the European fit -- nothing refitted -- because a caller tuned on one population reproducing that population's statistics proves nothing. DETECTION TRANSFERS: Europe (fitted) East Asia (new) above own null 60/60 30/30 sensitivity 31.6% 31.6% precision 32.2% 41.9% per-individual r +0.620 +0.545 Identical sensitivity and BETTER precision on a population the thresholds never saw. The calibration learned archaic structure, not European structure. THE REPORTED EXTENT ORDERS THE POPULATIONS BACKWARDS, and this is now the single reason the feature stays gated. Truth puts East Asian archaic extent at 1.217x Europe's; we call 0.937x. A user would be told an East Asian carries LESS archaic ancestry than a European. Cause: reported extent is true positives plus false positives, and the false-positive load is population-dependent (precision 32.2% vs 41.9%), so Europeans accumulate more spurious extent. A RESULT I ALMOST REPORTED AS A SUCCESS, recorded so it is not quoted later: "detected sequence reproduces 1.22x" is CIRCULAR. detected = sensitivity x truth, and sensitivity is equal across the two populations, so that ratio matches by construction. It restates the invariance; it does not test the ordering. The honest read is that per-individual detection transfers and the population-level number does not. THREE CAUSES RULED OUT, each measured rather than argued: - background contamination of p_background -- carrying rates 11.9% (EUR) vs 12.2% (EAS), and both emission states scale together so the contrast is preserved; - tract length -- median 29 kb in both populations; East Asians simply have MORE tracts (54 vs 46 per person), consistent with the 1.22x truth; - panel ascertainment -- in-tract contrast 2.99x (EUR) vs 3.04x (EAS), ratio 1.014, so the diagnostic panel is equally informative in both. This was the most plausible candidate and it is not the cause. Banks cross_population.py (the frozen-parameter transfer test, with a permutation test on the population difference) and observable/ascertainment.py (the per-population contrast measurement). Co-Authored-By: Claude Opus 5 (1M context) --- .../navigator-analysis/src/archaic_match.rs | 41 +++- .../archaic-validation/cross_population.py | 178 ++++++++++++++++++ .../observable/ascertainment.py | 147 +++++++++++++++ 3 files changed, 359 insertions(+), 7 deletions(-) create mode 100644 scripts/archaic-validation/cross_population.py create mode 100644 scripts/archaic-validation/observable/ascertainment.py diff --git a/crates/navigator-analysis/src/archaic_match.rs b/crates/navigator-analysis/src/archaic_match.rs index 8856a91..6d4b221 100644 --- a/crates/navigator-analysis/src/archaic_match.rs +++ b/crates/navigator-analysis/src/archaic_match.rs @@ -57,13 +57,40 @@ //! null (mean 45.3 % sensitivity against a 7.1 % null); the density caller scored 2.1 % against a //! 5.0 % null, i.e. below chance. //! -//! **Still not enough to re-enable**, and the limits are specific rather than general unease: -//! precision is 34.9 %, so two thirds of called sequence is not in the reference callset; the cohort -//! is **European only** and **chr21+22 only**; and the reference callset is itself weakly supported -//! (hmmix's own tracts are enriched only 1.84x for their own archaic SNPs), so agreement with it -//! caps out well below 100 % even for a correct caller. `ARCHAIC_SEGMENTS_ENABLED` stays `false` -//! until this reproduces outside Europe — East Asians are the sharp test, since the truth predicts -//! ~1.18x more archaic sequence there and a caller merely tracking European structure would miss it. +//! ## Cross-population: transfers per individual, but the reported number does not +//! +//! Run on 30 East Asians with the parameters **frozen** at the European fit, nothing refitted: +//! +//! | | Europe (fitted) | East Asia (new) | +//! |---|---|---| +//! | above own random-placement null | 60/60 | **30/30** | +//! | sensitivity | 31.6 % | **31.6 %** | +//! | precision | 32.2 % | **41.9 %** | +//! | per-individual extent `r` | +0.620 | **+0.545** | +//! +//! Detection transfers: identical sensitivity and *better* precision on a population the thresholds +//! never saw, so the calibration learned archaic structure rather than European structure. +//! +//! **But the reported extent orders the populations backwards.** The truth puts East Asian archaic +//! extent at **1.217x** Europe's; our called extent is **0.937x**. A user would be told an East +//! Asian carries *less* archaic ancestry than a European, which is the wrong way round and is the +//! single reason this is still gated. +//! +//! The cause is that reported extent is true positives *plus* false positives, and the false-positive +//! load is population-dependent (precision 32.2 % against 41.9 %), so Europeans accumulate more +//! spurious extent. Note that "detected sequence reproduces 1.22x" is **not** evidence to the +//! contrary: detected = sensitivity x truth, and sensitivity is equal across the two populations, so +//! that ratio matches by construction. It restates the invariance, it does not test the ordering. +//! +//! Ruled out as causes, each measured rather than argued: background contamination of `p_background` +//! (carrying rates 11.9 % vs 12.2 %, and both states scale together), tract length (median 29 kb in +//! both; East Asians simply have more tracts, 54 vs 46 per person), and panel ascertainment +//! (in-tract contrast 2.99x vs 3.04x, ratio 1.014 — the panel is equally informative in both). +//! +//! **Still not enough to re-enable.** Beyond the ordering: precision is 34.9 % on held-out +//! Europeans, the cohort is **chr21+22 only**, and the reference callset is itself weakly supported +//! (hmmix's own tracts are enriched just 1.84x for their own archaic SNPs), so agreement with it +//! caps well below 100 % even for a correct caller — F1 alone cannot say when this is done. use std::collections::BTreeMap; diff --git a/scripts/archaic-validation/cross_population.py b/scripts/archaic-validation/cross_population.py new file mode 100644 index 0000000..3cc58cd --- /dev/null +++ b/scripts/archaic-validation/cross_population.py @@ -0,0 +1,178 @@ +"""Does the caller transfer to East Asians, using parameters fitted only on Europeans? + +This is the sharp test, and it asks something the calibration could not have bought. Two things +have to hold: + + 1. TRANSFER -- per-individual locations and extent hold up on a population the thresholds were + never fitted to. If performance collapses, the calibration learned European structure. + + 2. THE POPULATION PREDICTION -- hmmix's own data puts East Asian archaic extent at 2.45 Mb + against Europe's 2.09, a ratio of ~1.18. Reproducing that ORDERING is a prediction the caller + was never shown: nothing in the fit knows which population a sample comes from. A caller that + merely reproduces whatever it was tuned on would return the same number for both. + +Parameters are frozen at the European-fitted values; nothing here is refitted. +""" + +import collections +import glob +import json +import math +import os +import random + +SC = os.path.dirname(os.path.abspath(__file__)) +V2 = os.path.join(SC, 'v2') +CONTIGS = ('chr21', 'chr22') +# European-fitted, frozen. +RATIO, MIN_POST, MIN_SITES, MIN_BP = '4.5', 0.98, 16, 5_000 + + +def union(iv, tol=0): + iv = sorted(iv) + out, cs, ce = [], *iv[0] + for s, e in iv[1:]: + if s > ce + tol: + out.append((cs, ce)) + cs, ce = s, e + else: + ce = max(ce, e) + out.append((cs, ce)) + return out + + +def load_bed(path, tol=0): + d = collections.defaultdict(list) + for line in open(path): + p = line.split() + if len(p) >= 3: + d[p[0]].append((int(p[1]), int(p[2]))) + return {c: union(v, tol) for c, v in d.items()} + + +def bp(d): + return sum(e - s for v in d.values() for s, e in v) + + +def intersect(a, b): + tot = 0 + for c, av in a.items(): + for s, e in av: + for s2, e2 in b.get(c, []): + lo, hi = max(s, s2), min(e, e2) + if hi > lo: + tot += hi - lo + return tot + + +def segs(doc): + oi = collections.defaultdict(list) + for s in doc['segments']: + if s['contig'] not in CONTIGS or s['posterior'] < MIN_POST or s['n_private'] < MIN_SITES: + continue + if s['end'] - s['start'] < MIN_BP: + continue + oi[s['contig']].append((s['start'], s['end'])) + return {c: union(v) for c, v in oi.items()} if oi else {} + + +def pearson(x, y): + n = len(x) + if n < 3: + return 0.0 + mx, my = sum(x) / n, sum(y) / n + num = sum((a - mx) * (b - my) for a, b in zip(x, y)) + dx = math.sqrt(sum((a - mx) ** 2 for a in x)) + dy = math.sqrt(sum((b - my) ** 2 for b in y)) + return num / (dx * dy) if dx and dy else 0.0 + + +def null_p95(ours, truth, tbp, draws=200, seed=0): + rnd = random.Random(seed) + lens = {c: [e - s for s, e in v] for c, v in ours.items()} + vals = [] + for _ in range(draws): + r = {} + for c, L in lens.items(): + if c not in truth: + continue + lo = min(s for s, _ in truth[c]) + hi = max(e for _, e in truth[c]) + r[c] = union([(p, p + l) for l in L for p in [rnd.randint(lo, max(lo, hi - l))]]) + vals.append(intersect(r, truth) / tbp * 100) + vals.sort() + return sum(vals) / len(vals), vals[int(0.95 * len(vals)) - 1] + + +def group(samples, source): + rows = [] + for s in samples: + f = f'{source}/{s}.json' + if not os.path.exists(f): + continue + doc = json.load(open(f)) + doc = doc[RATIO] if RATIO in doc else doc + ours = segs(doc) + truth = load_bed(f'{V2}/truth_{s}.bed', tol=1000) + if not truth: + continue + tb, ob = bp(truth), (bp(ours) if ours else 0) + ov = intersect(ours, truth) if ours else 0 + mean, p95 = null_p95(ours, truth, tb) if ours else (0, 0) + rows.append({'s': s, 't': tb / 1e6, 'o': ob / 1e6, 'ov': ov, + 'sens': ov / tb * 100, 'prec': (ov / ob * 100) if ob else 0, + 'null': mean, 'p95': p95}) + return rows + + +def summarize(rows, label): + if not rows: + print(f'{label}: no data') + return None + sens = sum(r['ov'] for r in rows) / sum(r['t'] * 1e6 for r in rows) * 100 + prec = sum(r['ov'] for r in rows) / max(sum(r['o'] * 1e6 for r in rows), 1) * 100 + f1 = 2 * sens * prec / (sens + prec) if sens + prec else 0 + beat = sum(1 for r in rows if r['sens'] > r['p95']) + t = [r['t'] for r in rows] + o = [r['o'] for r in rows] + print(f'{label:16s} n={len(rows):3d} F1 {f1:5.1f}% sens {sens:5.1f}% prec {prec:5.1f}% ' + f'ext {sum(o) / sum(t):5.2f} r {pearson(t, o):+6.3f} above-null {beat}/{len(rows)}') + return {'truth_mean': sum(t) / len(t), 'ours_mean': sum(o) / len(o), 'r': pearson(t, o), + 'n': len(rows), 't': t, 'o': o} + + +def main(): + eur = [l.strip() for l in open(f'{SC}/eur60.txt') if l.strip()] + eas = [l.strip() for l in open(f'{SC}/eas30.txt') if l.strip()] + print(f'parameters FROZEN at the European fit: ratio {RATIO}, posterior {MIN_POST}, ' + f'sites {MIN_SITES}, min_bp {MIN_BP}\n') + + e = summarize(group(eur, f'{V2}/sweep'), 'EUROPE (fitted)') + a = summarize(group(eas, f'{V2}/sweep_eas'), 'EAST ASIA (new)') + if not (e and a): + return + + print('\nTHE POPULATION PREDICTION (nothing in the fit knows a sample\'s population)') + print(f' hmmix truth EAS {a["truth_mean"]:.3f} Mb / EUR {e["truth_mean"]:.3f} Mb' + f' = {a["truth_mean"] / e["truth_mean"]:.3f}x') + print(f' our calls EAS {a["ours_mean"]:.3f} Mb / EUR {e["ours_mean"]:.3f} Mb' + f' = {a["ours_mean"] / e["ours_mean"]:.3f}x') + print(' published expectation: East Asians carry ~1.2x the Neanderthal ancestry of Europeans.') + + # Is the elevation real, or within noise? Permutation over the pooled samples. + pooled = e['o'] + a['o'] + labels = [0] * len(e['o']) + [1] * len(a['o']) + obs = (sum(a['o']) / len(a['o'])) - (sum(e['o']) / len(e['o'])) + rnd = random.Random(0) + hits = 0 + for _ in range(20000): + rnd.shuffle(labels) + ea = [v for v, l in zip(pooled, labels) if l == 1] + eu = [v for v, l in zip(pooled, labels) if l == 0] + if (sum(ea) / len(ea)) - (sum(eu) / len(eu)) >= obs: + hits += 1 + print(f' our EAS-EUR difference {obs:+.3f} Mb, permutation p = {(hits + 1) / 20001:.4f}') + + +if __name__ == '__main__': + main() diff --git a/scripts/archaic-validation/observable/ascertainment.py b/scripts/archaic-validation/observable/ascertainment.py new file mode 100644 index 0000000..171c5c4 --- /dev/null +++ b/scripts/archaic-validation/observable/ascertainment.py @@ -0,0 +1,147 @@ +"""Is the diagnostic-site panel less informative inside EAST ASIAN archaic tracts? + +The caller transfers per-individual (30/30 above null, precision 41.9%) but inverts the population +ordering: the truth puts East Asian archaic extent at 1.22x Europe's, and we call 0.94x. Ruled out +already: background contamination (carrying rates 11.9% vs 12.2%, and both states scale together) +and tract length (median 29 kb in both; EAS simply have MORE tracts). + +What remains is the observable itself. The model detects a tract by an elevated rate of carrying the +archaic allele at panel sites. If East Asian tracts carry those sites at a LOWER rate -- because the +panel's sites were ascertained on data that better represents the haplotypes introgressed into +Europeans -- then the same tract is less visible in an East Asian, and the caller under-calls exactly +where the truth says there is more. + +This measures the in-tract and background carrying rates per population. The contrast, not the +absolute rate, is what the model separates on. +""" + +import bisect +import collections +import json +import os +import subprocess +import sys + +SC = os.path.dirname(os.path.abspath(__file__)) +V2 = f'{SC}/v2' +DB = os.path.expanduser('~/.decodingus/navigator-rs.db') +CONTIGS = ('chr21', 'chr22') + + +def merge(iv, tol=0): + iv = sorted(iv) + out, cs, ce = [], *iv[0] + for s, e in iv[1:]: + if s > ce + tol: + out.append((cs, ce)) + cs, ce = s, e + else: + ce = max(ce, e) + out.append((cs, ce)) + return out + + +def load_bed(path, tol=0): + d = collections.defaultdict(list) + for line in open(path): + p = line.split() + if len(p) >= 3: + d[p[0]].append((int(p[1]), int(p[2]))) + return {c: merge(v, tol) for c, v in d.items()} + + +def sql(q): + return subprocess.run(['sqlite3', DB, q], capture_output=True, text=True).stdout.strip() + + +def calls_for(sample): + aln = sql(f"SELECT a.id FROM biosample b JOIN sequence_run r ON r.biosample_guid=b.guid " + f"JOIN alignment a ON a.sequence_run_id=r.id WHERE b.donor_identifier='{sample}' " + f"AND a.reference_build='chm13v2.0' LIMIT 1;") + if not aln: + return None + out = {} + for c in CONTIGS: + payload = sql(f"SELECT payload FROM analysis_artifact WHERE alignment_id={aln} " + f"AND kind='diploid_denovo:{c}';") + if not payload: + return None + for rec in json.loads(payload): + d = rec['dosage'] + if d is None or d < 0: + continue + alleles = set() + if d < 2: + alleles.add(rec['reference_allele']) + if d > 0: + alleles.add(rec['alternate_allele']) + out[(rec['contig'], rec['position'])] = alleles + return out + + +def main(): + diag = collections.defaultdict(dict) + with open(f'{SC}/classify.tsv') as f: + next(f) + for line in f: + c, p, d, k = line.rstrip('\n').split('\t') + diag[c][int(p)] = d + + callable_r = load_bed(f'{SC}/callable.bed') + + def in_reg(regs, p): + i = bisect.bisect_right([s for s, _ in regs], p) - 1 + return i >= 0 and regs[i][1] > p + + print(f"{'pop':4s} {'sample':10s} {'in-tract':>9s} {'background':>11s} {'contrast':>9s}") + agg = collections.defaultdict(lambda: [0, 0, 0, 0]) + for grp, listfile, n in (('EUR', 'eur60.txt', 8), ('EAS', 'eas30.txt', 8)): + for s in [x.strip() for x in open(f'{SC}/{listfile}') if x.strip()][:n]: + tpath = f'{V2}/truth_{s}.bed' + if not os.path.exists(tpath): + continue + carries = calls_for(s) + if carries is None: + continue + truth = load_bed(tpath, tol=1000) + hits = [0, 0, 0, 0] # in_carried, in_total, bg_carried, bg_total + for c in CONTIGS: + cal = callable_r.get(c, []) + tr = truth.get(c, []) + for p, derived in diag[c].items(): + if not cal or not in_reg(cal, p): + continue + got = carries.get((c, p)) + carried = got is not None and derived in got + if tr and in_reg(tr, p): + hits[1] += 1 + hits[0] += carried + else: + hits[3] += 1 + hits[2] += carried + if hits[1] == 0 or hits[3] == 0: + continue + r_in = hits[0] / hits[1] * 100 + r_bg = hits[2] / hits[3] * 100 + print(f'{grp:4s} {s:10s} {r_in:8.1f}% {r_bg:10.1f}% {r_in / r_bg:8.2f}x') + for i in range(4): + agg[grp][i] += hits[i] + + print() + for grp in ('EUR', 'EAS'): + a = agg[grp] + if a[1] and a[3]: + r_in, r_bg = a[0] / a[1] * 100, a[2] / a[3] * 100 + print(f'{grp} POOLED in-tract {r_in:.1f}% background {r_bg:.1f}% ' + f'CONTRAST {r_in / r_bg:.2f}x') + e, a = agg['EUR'], agg['EAS'] + if e[1] and a[1]: + ce = (e[0] / e[1]) / (e[2] / e[3]) + ca = (a[0] / a[1]) / (a[2] / a[3]) + print(f'\ncontrast EAS/EUR = {ca / ce:.3f}') + print(' < 1 means the panel is less informative inside East Asian tracts, which would') + print(' explain under-calling exactly where the truth says there is MORE archaic DNA.') + + +if __name__ == '__main__': + main() From c4091ac8be060d852c7ce7836e03a8f4c079243c Mon Sep 17 00:00:00 2001 From: James Kane Date: Fri, 31 Jul 2026 11:23:08 -0500 Subject: [PATCH 6/9] Archaic Tier B v2: an arbiter that does not ask another caller's opinion Precision has been scored against hmmix, but a call they did not make is not automatically wrong -- their callset is incomplete by an unknown amount (their own tracts are enriched just 1.84x for their own archaic SNPs). Settling that needed a referee independent of both callers. The Tier A panel records, per site, which of the four archaic genomes carries the derived allele. The Tier B caller NEVER SEES THIS -- it reads only a derived base and a lineage class from ArchaicClassify -- so per-genome concordance is evidence it cannot have been fitted to. Of the sites where a given archaic genome is derived, what fraction does the subject carry: true positive FALSE positive background Europe 93.6% 81.3% 59.0% East Asia 93.5% 72.9% 45.5% Our "false positives" sit 64% (EUR) and 57% (EAS) of the way from background to true positive. They are a MIXTURE: real tracts hmmix missed, genuine noise, and calls that are correctly placed but over-extended. So precision against hmmix understates this caller -- though not enough to dismiss it, and F1 stays a usable objective, which means the calibration was not invalidated. THREE ATTEMPTS, TWO OF THEM NON-ANSWERS, all kept because each looked convincing: - carrying rate inside segments: CIRCULAR. The caller selects on carrying rate, so every segment it emits has a high one. It reported false positives at 98% of the way to true positive, which is meaningless. - overlap with other individuals' tracts: SATURATED. Across 90 individuals the hmmix tracts cover 67% of callable territory, so a random region hits them ~60% of the time and the observed excess was 2-7 points. - the first version of the arbiter itself: VACUOUS. It conditioned on the subject already carrying the derived allele, but at a discordant site at least one genome is derived by construction, so best-genome concordance was ~100% EVERYWHERE INCLUDING BACKGROUND. Fixed by conditioning on the GENOME instead -- of the sites where genome g is derived, how many does the subject carry -- at which point background correctly lands near the genome-wide carrying rate. A LEAD ON THE ORDERING INVERSION: background concordance differs by population, 59.0% in Europeans against 45.5% in East Asians. Europeans carry archaic-derived alleles more often OUTSIDE tracts, which is a plausible mechanism for the population-varying false-positive load and hence for the reported extent ordering the populations backwards. Adds archaic_panel_dump (the per-genome calls as TSV) and banks all three attempts under scripts/archaic-validation/, since the two failures are the more instructive half. Co-Authored-By: Claude Opus 5 (1M context) --- .../examples/archaic_panel_dump.rs | 64 ++++++ .../navigator-analysis/src/archaic_match.rs | 25 ++ scripts/archaic-validation/arbiter.py | 217 ++++++++++++++++++ .../observable/fp_character.py | 176 ++++++++++++++ .../observable/fp_shared_haplotypes.py | 121 ++++++++++ 5 files changed, 603 insertions(+) create mode 100644 crates/navigator-analysis/examples/archaic_panel_dump.rs create mode 100644 scripts/archaic-validation/arbiter.py create mode 100644 scripts/archaic-validation/observable/fp_character.py create mode 100644 scripts/archaic-validation/observable/fp_shared_haplotypes.py diff --git a/crates/navigator-analysis/examples/archaic_panel_dump.rs b/crates/navigator-analysis/examples/archaic_panel_dump.rs new file mode 100644 index 0000000..c9ea5d4 --- /dev/null +++ b/crates/navigator-analysis/examples/archaic_panel_dump.rs @@ -0,0 +1,64 @@ +//! Dump the Tier A marker panel with its **per-archaic-genome** calls, as TSV. +//! +//! This is the independent evidence for arbitrating Tier B calls. The segment caller +//! ([`navigator_analysis::archaic_match`]) reads only `ArchaicClassify` — a derived base and a +//! lineage class per site — and never sees which archaic genome carries what. So the per-genome +//! pattern is information the caller cannot have fitted to, which is what makes it usable as a +//! referee. +//! +//! Why a referee is needed: precision has been measured against hmmix's callset, but a call absent +//! from hmmix is not necessarily wrong — hmmix's own tracts are enriched only 1.84x for their own +//! archaic SNPs, so that callset is incomplete by an unknown amount. Scoring a segment against the +//! archaic genomes directly asks whether it looks like an inherited archaic haplotype, without +//! asking another caller's opinion. +//! +//! ```sh +//! cargo run --release -p navigator-analysis --example archaic_panel_dump -- \ +//! ~/.decodingus/ancestry/archaic_markers_chm13v2.0.bin chr21 chr22 > panel.tsv +//! ``` + +use navigator_analysis::archaic::{ArchaicMarkerPanel, ARCHAIC_GENOMES}; + +fn main() -> Result<(), Box> { + let mut a = std::env::args().skip(1); + let path = a.next().expect("usage: archaic_panel_dump [contig ...]"); + let want: Vec = a.collect(); + + let panel = ArchaicMarkerPanel::from_bytes(&std::fs::read(&path)?).map_err(|e| e.to_string())?; + eprintln!("panel: {} sites, build {}", panel.sites.len(), panel.build); + + // One column per archaic genome: D = carries the derived allele, A = positively called + // homozygous-ancestral, . = no call. The A/. distinction is load-bearing — treating a no-call as + // ancestral is the error that produced ~19 % Denisovan for a European in an earlier pass. + println!("contig\tposition\tderived\tclass\t{}", ARCHAIC_GENOMES.join("\t")); + let mut n = 0usize; + for s in &panel.sites { + if !want.is_empty() && !want.contains(&s.contig) { + continue; + } + let calls: Vec<&str> = s + .calls + .iter() + .map(|c| { + if c.carries_derived() { + "D" + } else if matches!(c, navigator_analysis::archaic::ArchaicCall::HomAncestral) { + "A" + } else { + "." + } + }) + .collect(); + println!( + "{}\t{}\t{}\t{:?}\t{}", + s.contig, + s.position, + s.archaic_derived_allele, + s.diagnostic_class, + calls.join("\t") + ); + n += 1; + } + eprintln!("wrote {n} sites"); + Ok(()) +} diff --git a/crates/navigator-analysis/src/archaic_match.rs b/crates/navigator-analysis/src/archaic_match.rs index 6d4b221..6fe91ea 100644 --- a/crates/navigator-analysis/src/archaic_match.rs +++ b/crates/navigator-analysis/src/archaic_match.rs @@ -87,6 +87,31 @@ //! both; East Asians simply have more tracts, 54 vs 46 per person), and panel ascertainment //! (in-tract contrast 2.99x vs 3.04x, ratio 1.014 — the panel is equally informative in both). //! +//! ## How much of the "false positive" rate is really ours +//! +//! Precision is measured against hmmix, but a call they did not make is not automatically wrong. +//! An independent arbiter settles this without asking another caller: the Tier A panel records, per +//! site, which of the four archaic genomes carries the derived allele, and **this caller never sees +//! that** — it reads only a derived base and a lineage class. So per-genome concordance is evidence +//! it cannot have been fitted to. +//! +//! Of the sites where a given archaic genome is derived, what fraction does the subject carry +//! (best-matching genome): +//! +//! | | true positive | false positive | background | +//! |---|---|---|---| +//! | Europe | 93.6 % | **81.3 %** | 59.0 % | +//! | East Asia | 93.5 % | **72.9 %** | 45.5 % | +//! +//! Our "false positives" sit **64 % / 57 %** of the way from background to true positive. They are a +//! mixture: real tracts hmmix missed, plus genuine noise, plus calls that are correctly placed but +//! over-extended. So precision against hmmix **understates** this caller — though not enough to +//! dismiss it, and F1 remains a usable objective. +//! +//! Note the background rates differ by population (59.0 % against 45.5 %): Europeans carry +//! archaic-derived alleles more often *outside* tracts. That is a candidate mechanism for the +//! population-varying false-positive load, and hence for the ordering inversion above. +//! //! **Still not enough to re-enable.** Beyond the ordering: precision is 34.9 % on held-out //! Europeans, the cohort is **chr21+22 only**, and the reference callset is itself weakly supported //! (hmmix's own tracts are enriched just 1.84x for their own archaic SNPs), so agreement with it diff --git a/scripts/archaic-validation/arbiter.py b/scripts/archaic-validation/arbiter.py new file mode 100644 index 0000000..82838f3 --- /dev/null +++ b/scripts/archaic-validation/arbiter.py @@ -0,0 +1,217 @@ +"""An arbiter for Tier B calls that does not ask another caller's opinion. + +Precision has been measured against hmmix, but a call absent from hmmix is not necessarily wrong: +their callset is incomplete by an unknown amount (their own tracts are enriched just 1.84x for their +own archaic SNPs). Two attempts to settle whether our extra calls are real both failed -- +carrying rate is circular (the caller selects on it), and shared-haplotype overlap is saturated +(hmmix tracts cover 67% of callable territory, so the null is ~60%). + +This asks the archaic genomes instead. A real introgressed tract was inherited from ONE archaic +individual, so the derived alleles the subject carries inside it should concentrate on the genomes +that share that haplotype -- and in particular should be present at sites where ONE genome is +derived and another is positively called ancestral. Those discordant sites are the informative ones: +a site where all four archaics are derived says nothing about which haplotype was inherited. + +Crucially, the caller never sees per-genome calls -- it reads only a derived base and a lineage +class from `ArchaicClassify` -- so this is evidence it cannot have fitted to. + +Reports, for true positives, false positives and background: + * DISCORDANT-SITE concordance -- of the sites where the subject carries the derived allele and the + archaic genomes DISAGREE, what fraction match the best-matching genome. Random carriage gives + the base rate; an inherited haplotype gives much more. +""" + +import bisect +import collections +import json +import os +import subprocess + +SC = os.path.dirname(os.path.abspath(__file__)) +V2 = f'{SC}/v2' +DB = os.path.expanduser('~/.decodingus/navigator-rs.db') +CONTIGS = ('chr21', 'chr22') +RATIO, MIN_POST, MIN_SITES, MIN_BP = '4.5', 0.98, 16, 5_000 +GENOMES = ['AltaiNeanderthal', 'Vindija33.19', 'Chagyrskaya8', 'Denisova3'] + + +def merge(iv, tol=0): + iv = sorted(iv) + out, cs, ce = [], *iv[0] + for s, e in iv[1:]: + if s > ce + tol: + out.append((cs, ce)) + cs, ce = s, e + else: + ce = max(ce, e) + out.append((cs, ce)) + return out + + +def load_bed(path, tol=0): + d = collections.defaultdict(list) + for line in open(path): + p = line.split() + if len(p) >= 3: + d[p[0]].append((int(p[1]), int(p[2]))) + return {c: merge(v, tol) for c, v in d.items()} + + +def sql(q): + return subprocess.run(['sqlite3', DB, q], capture_output=True, text=True).stdout.strip() + + +def carries_for(sample): + aln = sql(f"SELECT a.id FROM biosample b JOIN sequence_run r ON r.biosample_guid=b.guid " + f"JOIN alignment a ON a.sequence_run_id=r.id WHERE b.donor_identifier='{sample}' " + f"AND a.reference_build='chm13v2.0' LIMIT 1;") + if not aln: + return None + out = {} + for c in CONTIGS: + payload = sql(f"SELECT payload FROM analysis_artifact WHERE alignment_id={aln} " + f"AND kind='diploid_denovo:{c}';") + if not payload: + return None + for rec in json.loads(payload): + d = rec['dosage'] + if d is None or d < 0: + continue + al = set() + if d < 2: + al.add(rec['reference_allele']) + if d > 0: + al.add(rec['alternate_allele']) + out[(rec['contig'], rec['position'])] = al + return out + + +def load_panel(): + """Sites where the archaic genomes DISAGREE — the only ones that identify a haplotype.""" + by_contig = collections.defaultdict(list) + with open(f'{SC}/panel.tsv') as f: + hdr = next(f).rstrip('\n').split('\t') + gi = [hdr.index(g) for g in GENOMES] + for line in f: + p = line.rstrip('\n').split('\t') + calls = [p[i] for i in gi] + if 'D' in calls and 'A' in calls: # informative: some derived, some POSITIVELY ancestral + by_contig[p[0]].append((int(p[1]), p[2], calls)) + return {c: sorted(v) for c, v in by_contig.items()} + + +def overlaps(regs, s, e): + for a, b in regs: + if min(e, b) > max(s, a): + return True + return False + + +def concordance(panel, contig, s, e, carries): + """For each archaic genome g: of the sites where g is DERIVED, how many does the subject carry? + + Conditioning on the GENOME, not on the subject. An earlier version conditioned on the subject + already carrying the derived allele, which is vacuous: at a discordant site at least one genome + is derived by construction, so the best-matching genome scored ~100 % everywhere including + background, and the statistic separated nothing. + + Read this way the quantity is the subject's sensitivity to a specific archaic haplotype: + background sits at the genome-wide carrying rate (~13 %), while a tract inherited from that + lineage should be far higher. + """ + v = panel.get(contig, []) + keys = [x[0] for x in v] + lo, hi = bisect.bisect_left(keys, s), bisect.bisect_left(keys, e) + hits = [0] * len(GENOMES) + dens = [0] * len(GENOMES) + for pos, derived, calls in v[lo:hi]: + got = carries.get((contig, pos)) + subject_has = bool(got) and derived in got + for i, c in enumerate(calls): + if c != 'D': + continue + dens[i] += 1 + if subject_has: + hits[i] += 1 + # Best genome by rate, requiring a few sites so a 1/1 does not win. + best_rate, best_hits, best_den = 0.0, 0, 0 + for h, d in zip(hits, dens): + if d >= 3 and (h / d) > best_rate: + best_rate, best_hits, best_den = h / d, h, d + return best_hits, best_den + + +def main(): + panel = load_panel() + n_inf = sum(len(v) for v in panel.values()) + print(f'informative (discordant) panel sites on chr21+22: {n_inf}\n') + + callable_r = load_bed(f'{SC}/callable.bed') + agg = collections.defaultdict(lambda: [0, 0, 0]) # best_hits, total, n_segments + + for grp, listfile, source, n in (('EUR', 'eur60.txt', f'{V2}/sweep', 15), + ('EAS', 'eas30.txt', f'{V2}/sweep_eas', 15)): + for s in [x.strip() for x in open(f'{SC}/{listfile}') if x.strip()][:n]: + f = f'{source}/{s}.json' + tp = f'{V2}/truth_{s}.bed' + if not (os.path.exists(f) and os.path.exists(tp)): + continue + carries = carries_for(s) + if carries is None: + continue + doc = json.load(open(f)) + doc = doc[RATIO] if RATIO in doc else doc + truth = load_bed(tp, tol=1000) + called = collections.defaultdict(list) + for seg in doc['segments']: + c = seg['contig'] + if c not in CONTIGS or seg['posterior'] < MIN_POST or seg['n_private'] < MIN_SITES: + continue + if seg['end'] - seg['start'] < MIN_BP: + continue + called[c].append((seg['start'], seg['end'])) + cls = 'true positive' if overlaps(truth.get(c, []), seg['start'], seg['end']) \ + else 'FALSE positive' + h, d = concordance(panel, c, seg['start'], seg['end'], carries) + if d >= 3: + a = agg[(grp, cls)] + a[0] += h + a[1] += d + a[2] += 1 + # background: callable windows we did NOT call and hmmix did not either + for c in CONTIGS: + cal = callable_r.get(c, []) + mine = merge(called[c]) if called[c] else [] + for a0, b0 in cal[:600]: + if overlaps(mine, a0, b0) or overlaps(truth.get(c, []), a0, b0): + continue + h, d = concordance(panel, c, a0, b0, carries) + if d >= 3: + a = agg[(grp, 'background')] + a[0] += h + a[1] += d + a[2] += 1 + + print(f"{'pop':4s} {'class':16s} {'regions':>8s} {'sites':>7s} {'best-genome concordance':>24s}") + for grp in ('EUR', 'EAS'): + for cls in ('true positive', 'FALSE positive', 'background'): + a = agg[(grp, cls)] + if a[1]: + print(f'{grp:4s} {cls:16s} {a[2]:8d} {a[1]:7d} {a[0] / a[1] * 100:23.1f}%') + print() + for grp in ('EUR', 'EAS'): + tp, fp, bg = agg[(grp, 'true positive')], agg[(grp, 'FALSE positive')], agg[(grp, 'background')] + if tp[1] and fp[1] and bg[1]: + r = lambda a: a[0] / a[1] + span = r(tp) - r(bg) + pos = (r(fp) - r(bg)) / span if span > 0 else 0 + print(f'{grp}: false positives sit {pos * 100:.0f}% of the way from background to true ' + f'positive on ARCHAIC-GENOME concordance') + print('\n This statistic is not one the caller optimises: it reads per-genome calls, which the') + print(' caller never sees. Near 100% means our extra calls carry the same haplotype signature') + print(' as confirmed tracts, and precision against hmmix understates us. Near 0% means they') + print(' are noise and the precision figure is fair.') + + +if __name__ == '__main__': + main() diff --git a/scripts/archaic-validation/observable/fp_character.py b/scripts/archaic-validation/observable/fp_character.py new file mode 100644 index 0000000..d04bec7 --- /dev/null +++ b/scripts/archaic-validation/observable/fp_character.py @@ -0,0 +1,176 @@ +"""Are our "false positives" archaic-looking, or noise? And does that differ by population? + +The blocker is that reported extent orders the populations backwards (0.937x where the truth says +1.217x), driven by precision differing between them (32.2% EUR vs 41.9% EAS) at identical +thresholds and identical in-tract contrast. Two very different explanations fit that: + + (a) OUR FAULT -- we emit more spurious calls in Europeans. Those segments should look like + background: carrying rate near 13%. + + (b) THE TRUTH'S FAULT -- hmmix is less sensitive in Europeans, so real tracts we find are + scored as false positives. Those segments should look like true positives: carrying rate + near 40%. + +The archaic-allele carrying rate inside each segment class separates the two, and it is measured +independently of hmmix -- it uses only our calls and the diagnostic panel, so it does not assume +the reference callset is right. +""" + +import bisect +import collections +import json +import os +import subprocess + +SC = os.path.dirname(os.path.abspath(__file__)) +V2 = f'{SC}/v2' +DB = os.path.expanduser('~/.decodingus/navigator-rs.db') +CONTIGS = ('chr21', 'chr22') +RATIO, MIN_POST, MIN_SITES, MIN_BP = '4.5', 0.98, 16, 5_000 + + +def merge(iv, tol=0): + iv = sorted(iv) + out, cs, ce = [], *iv[0] + for s, e in iv[1:]: + if s > ce + tol: + out.append((cs, ce)) + cs, ce = s, e + else: + ce = max(ce, e) + out.append((cs, ce)) + return out + + +def load_bed(path, tol=0): + d = collections.defaultdict(list) + for line in open(path): + p = line.split() + if len(p) >= 3: + d[p[0]].append((int(p[1]), int(p[2]))) + return {c: merge(v, tol) for c, v in d.items()} + + +def sql(q): + return subprocess.run(['sqlite3', DB, q], capture_output=True, text=True).stdout.strip() + + +def carries_for(sample): + aln = sql(f"SELECT a.id FROM biosample b JOIN sequence_run r ON r.biosample_guid=b.guid " + f"JOIN alignment a ON a.sequence_run_id=r.id WHERE b.donor_identifier='{sample}' " + f"AND a.reference_build='chm13v2.0' LIMIT 1;") + if not aln: + return None + out = {} + for c in CONTIGS: + payload = sql(f"SELECT payload FROM analysis_artifact WHERE alignment_id={aln} " + f"AND kind='diploid_denovo:{c}';") + if not payload: + return None + for rec in json.loads(payload): + d = rec['dosage'] + if d is None or d < 0: + continue + al = set() + if d < 2: + al.add(rec['reference_allele']) + if d > 0: + al.add(rec['alternate_allele']) + out[(rec['contig'], rec['position'])] = al + return out + + +def overlaps(regs, s, e): + for a, b in regs: + if min(e, b) > max(s, a): + return True + return False + + +def main(): + diag = collections.defaultdict(dict) + with open(f'{SC}/classify.tsv') as f: + next(f) + for line in f: + c, p, d, _k = line.rstrip('\n').split('\t') + diag[c][int(p)] = d + diag_sorted = {c: sorted(v) for c, v in diag.items()} + callable_r = load_bed(f'{SC}/callable.bed') + + def rate(contig, s, e, carries): + keys = diag_sorted.get(contig, []) + lo = bisect.bisect_left(keys, s) + hi = bisect.bisect_left(keys, e) + tot = hit = 0 + for p in keys[lo:hi]: + tot += 1 + got = carries.get((contig, p)) + if got and diag[contig][p] in got: + hit += 1 + return hit, tot + + print(f"{'pop':4s} {'class':18s} {'segments':>9s} {'Mb':>7s} {'carrying rate':>14s}") + agg = collections.defaultdict(lambda: [0, 0, 0, 0.0]) + for grp, listfile, source, n in (('EUR', 'eur60.txt', f'{V2}/sweep', 12), + ('EAS', 'eas30.txt', f'{V2}/sweep_eas', 12)): + for s in [x.strip() for x in open(f'{SC}/{listfile}') if x.strip()][:n]: + f = f'{source}/{s}.json' + tp = f'{V2}/truth_{s}.bed' + if not (os.path.exists(f) and os.path.exists(tp)): + continue + doc = json.load(open(f)) + doc = doc[RATIO] if RATIO in doc else doc + truth = load_bed(tp, tol=1000) + carries = carries_for(s) + if carries is None: + continue + for seg in doc['segments']: + c = seg['contig'] + if c not in CONTIGS or seg['posterior'] < MIN_POST or seg['n_private'] < MIN_SITES: + continue + if seg['end'] - seg['start'] < MIN_BP: + continue + cls = 'true positive' if overlaps(truth.get(c, []), seg['start'], seg['end']) \ + else 'FALSE positive' + hit, tot = rate(c, seg['start'], seg['end'], carries) + a = agg[(grp, cls)] + a[0] += 1 + a[1] += hit + a[2] += tot + a[3] += (seg['end'] - seg['start']) / 1e6 + # background: the callable territory outside both our calls and the truth + ours = merge([(x['start'], x['end']) for x in doc['segments'] + if x['contig'] == 'chr21' and x['posterior'] >= MIN_POST + and x['n_private'] >= MIN_SITES] or [(0, 0)]) + for a0, b0 in callable_r.get('chr21', [])[:400]: + if overlaps(ours, a0, b0) or overlaps(truth.get('chr21', []), a0, b0): + continue + hit, tot = rate('chr21', a0, b0, carries) + a = agg[(grp, 'background')] + a[1] += hit + a[2] += tot + + for grp in ('EUR', 'EAS'): + for cls in ('true positive', 'FALSE positive', 'background'): + a = agg[(grp, cls)] + if not a[2]: + continue + print(f'{grp:4s} {cls:18s} {a[0]:9d} {a[3]:7.2f} {a[1] / a[2] * 100:13.1f}%') + print() + for grp in ('EUR', 'EAS'): + tp = agg[(grp, 'true positive')] + fp = agg[(grp, 'FALSE positive')] + bg = agg[(grp, 'background')] + if tp[2] and fp[2] and bg[2]: + r_tp, r_fp, r_bg = tp[1] / tp[2], fp[1] / fp[2], bg[1] / bg[2] + # 1.0 = our false positives look exactly like real tracts; 0.0 = like background. + pos = (r_fp - r_bg) / (r_tp - r_bg) if r_tp > r_bg else 0 + print(f'{grp}: false positives sit {pos * 100:.0f}% of the way from background to ' + f'true positive') + print('\n Near 0% => they are noise and the precision gap is ours.') + print(' Near 100% => they are archaic-looking and hmmix simply did not call them, which makes') + print(' "precision" against this reference a measure of the reference as much as of us.') + + +if __name__ == '__main__': + main() diff --git a/scripts/archaic-validation/observable/fp_shared_haplotypes.py b/scripts/archaic-validation/observable/fp_shared_haplotypes.py new file mode 100644 index 0000000..a51f927 --- /dev/null +++ b/scripts/archaic-validation/observable/fp_shared_haplotypes.py @@ -0,0 +1,121 @@ +"""Do our "false positives" land where OTHER people have archaic tracts? + +The carrying-rate version of this question is circular: the caller selects regions of high carrying +rate, so every segment it emits has one. This test does not use our caller's own evidence at all. + +Introgressed haplotypes are SHARED -- the same archaic haplotype segregates across a population. So +if a segment we call and hmmix did not is really archaic, it should coincide with tracts hmmix calls +in OTHER individuals. If instead it is our noise, it should fall where hmmix never calls anything, +i.e. no better than a random region of the same size. + +The individual's own truth is excluded from the union, so a segment cannot vindicate itself. +""" + +import collections +import json +import os +import random + +SC = os.path.dirname(os.path.abspath(__file__)) +V2 = f'{SC}/v2' +CONTIGS = ('chr21', 'chr22') +RATIO, MIN_POST, MIN_SITES, MIN_BP = '4.5', 0.98, 16, 5_000 + + +def merge(iv, tol=0): + iv = sorted(iv) + out, cs, ce = [], *iv[0] + for s, e in iv[1:]: + if s > ce + tol: + out.append((cs, ce)) + cs, ce = s, e + else: + ce = max(ce, e) + out.append((cs, ce)) + return out + + +def load_bed(path, tol=0): + d = collections.defaultdict(list) + for line in open(path): + p = line.split() + if len(p) >= 3: + d[p[0]].append((int(p[1]), int(p[2]))) + return {c: merge(v, tol) for c, v in d.items()} + + +def overlaps(regs, s, e): + for a, b in regs: + if min(e, b) > max(s, a): + return True + return False + + +def main(): + eur = [x.strip() for x in open(f'{SC}/eur60.txt') if x.strip()] + eas = [x.strip() for x in open(f'{SC}/eas30.txt') if x.strip()] + truths = {} + for s in eur + eas: + p = f'{V2}/truth_{s}.bed' + if os.path.exists(p): + truths[s] = load_bed(p, tol=1000) + + # The population-wide map of where archaic tracts occur at all, per contig. + pool = collections.defaultdict(list) + for s, t in truths.items(): + for c, v in t.items(): + pool[c].extend(v) + print(f'population archaic map from {len(truths)} individuals: ' + ', '.join( + f'{c} {sum(e - s for s, e in merge(v)) / 1e6:.1f} Mb' for c, v in sorted(pool.items()))) + + callable_r = load_bed(f'{SC}/callable.bed') + rnd = random.Random(0) + print(f"\n{'pop':4s} {'class':16s} {'n':>5s} {'in others truth':>16s} {'random null':>12s}") + for grp, names, source in (('EUR', eur[:25], f'{V2}/sweep'), ('EAS', eas[:25], f'{V2}/sweep_eas')): + counts = collections.Counter() + tot = collections.Counter() + null_hit = null_tot = 0 + for s in names: + f = f'{source}/{s}.json' + if not (os.path.exists(f) and s in truths): + continue + doc = json.load(open(f)) + doc = doc[RATIO] if RATIO in doc else doc + # union of everyone ELSE's tracts + others = {} + for c in CONTIGS: + iv = [x for o, t in truths.items() if o != s for x in t.get(c, [])] + others[c] = merge(iv) if iv else [] + for seg in doc['segments']: + c = seg['contig'] + if c not in CONTIGS or seg['posterior'] < MIN_POST or seg['n_private'] < MIN_SITES: + continue + if seg['end'] - seg['start'] < MIN_BP: + continue + cls = 'true positive' if overlaps(truths[s].get(c, []), seg['start'], seg['end']) \ + else 'FALSE positive' + tot[cls] += 1 + if overlaps(others[c], seg['start'], seg['end']): + counts[cls] += 1 + if cls == 'FALSE positive': + # a same-size region placed at random in callable territory + L = seg['end'] - seg['start'] + cal = callable_r.get(c, []) + if cal: + lo = min(a for a, _ in cal) + hi = max(b for _, b in cal) + p = rnd.randint(lo, max(lo, hi - L)) + null_tot += 1 + null_hit += overlaps(others[c], p, p + L) + for cls in ('true positive', 'FALSE positive'): + if tot[cls]: + extra = (f'{null_hit / null_tot * 100:11.1f}%' + if cls == 'FALSE positive' and null_tot else ' ' * 12) + print(f'{grp:4s} {cls:16s} {tot[cls]:5d} {counts[cls] / tot[cls] * 100:15.1f}%{extra}') + print('\n If FALSE positives hit other people\'s tracts far above the random null, they are') + print(' real archaic haplotypes hmmix missed in THIS individual -- and precision against') + print(' hmmix is then measuring their sensitivity, not our specificity.') + + +if __name__ == '__main__': + main() From 3d6a9719e5f0f3bc751c01ea65b429f27d6310e3 Mon Sep 17 00:00:00 2001 From: James Kane Date: Fri, 31 Jul 2026 11:42:34 -0500 Subject: [PATCH 7/9] Archaic Tier B v2: concordance filter takes precision to 90%, and finds a harder limit Scoring each called segment against the archaic genomes and dropping poor matches raises PRECISION FROM 54% TO 90%. The filter is sound on a held-out check: with Denisova excluded from it entirely, kept segments score 74.9% on Denisova concordance against 21.5% for dropped ones -- a 3.5x separation on a genome the filter never saw. IT DOES NOT FIX THE POPULATION ORDERING, and tightening it makes the ordering worse (0.927 unfiltered -> 0.895 at 90% precision -> 0.831 at the strictest threshold, against a 1.098 target on this subset). At 90% precision the reported extent is mostly true positives, so false positives are no longer the explanation. What remains is recovery: ~46% of European truth against ~38% of East Asian. THE CAUSE, and it is not a threshold: EUR true positives Neanderthal-best 89.2% Denisova 69.3% EAS true positives Neanderthal-best 83.4% Denisova 77.4% Denisova beats every Neanderthal: EUR 11.2% EAS 32.2% That 2.9x is the Denisovan ancestry East Asians carry and Europeans essentially lack -- the data reproduces known biology, which is reassuring about the observable. But it also means our four sequenced archaic genomes UNDER-REPRESENT EAST ASIAN ARCHAIC DIVERSITY, so any reference-based filter under-calls East Asians. Fixing it would take archaic genomes closer to the populations that introgressed into East Asia, which do not exist. So a cross-population comparable number is not achievable this way. The caller is defensible WITHIN a population and not BETWEEN them, and the module says so. HOLDING DENISOVA OUT WAS THE WRONG CHOICE, for a reason worth recording: it is not a neutral held-out genome. Excluding it biases the filter against exactly the population that carries Denisovan ancestry, and the best-of-three ordering (0.824) is worse than best-of-four (0.895). A held-out validator has to be neutral with respect to the thing being measured, and this one was not. Also records why p_background cannot fix the ordering: the model's estimate is flat across populations (14.17% EUR vs 14.04% EAS, ratio 0.991) while the false-positive load is not (ratio 0.877), so the parameter is blind to the driver. Banks concordance_filter.py (the filter with the held-out design) and observable/ordering_cause.py (the p_background diagnosis). Co-Authored-By: Claude Opus 5 (1M context) --- .../navigator-analysis/src/archaic_match.rs | 33 ++- .../archaic-validation/concordance_filter.py | 195 ++++++++++++++++++ .../observable/ordering_cause.py | 169 +++++++++++++++ 3 files changed, 393 insertions(+), 4 deletions(-) create mode 100644 scripts/archaic-validation/concordance_filter.py create mode 100644 scripts/archaic-validation/observable/ordering_cause.py diff --git a/crates/navigator-analysis/src/archaic_match.rs b/crates/navigator-analysis/src/archaic_match.rs index 6fe91ea..6d70fc7 100644 --- a/crates/navigator-analysis/src/archaic_match.rs +++ b/crates/navigator-analysis/src/archaic_match.rs @@ -112,10 +112,35 @@ //! archaic-derived alleles more often *outside* tracts. That is a candidate mechanism for the //! population-varying false-positive load, and hence for the ordering inversion above. //! -//! **Still not enough to re-enable.** Beyond the ordering: precision is 34.9 % on held-out -//! Europeans, the cohort is **chr21+22 only**, and the reference callset is itself weakly supported -//! (hmmix's own tracts are enriched just 1.84x for their own archaic SNPs), so agreement with it -//! caps well below 100 % even for a correct caller — F1 alone cannot say when this is done. +//! ## A concordance filter fixes precision, and exposes a harder limit +//! +//! Scoring each called segment against the archaic genomes and dropping the poor matches raises +//! **precision from 54 % to 90 %**. The filter is sound: with Denisova held out of it entirely, kept +//! segments score 74.9 % on Denisova concordance against 21.5 % for dropped ones — a 3.5x separation +//! on a genome the filter never saw. +//! +//! It does **not** fix the population ordering, and tightening it makes the ordering worse. At 90 % +//! precision the reported extent is mostly true positives, and it still orders the populations +//! backwards, so the cause is no longer false positives. What remains is a difference in *recovery*: +//! roughly 46 % of European truth against 38 % of East Asian. +//! +//! The reason is visible in the concordance itself. East Asian tracts match our archaic genomes less +//! well than European ones (83.4 % against 89.2 %), and Denisova is the best match for **32.2 % of +//! East Asian tracts against 11.2 % of European** ones. That 2.9x is the known Denisovan ancestry +//! East Asians carry and Europeans essentially lack — the data reproduces it — but it also means our +//! four sequenced archaic genomes **under-represent East Asian archaic diversity**. Any +//! reference-based filter therefore under-calls East Asians, and holding Denisova out (the first +//! design here) makes it markedly worse. +//! +//! That is a limit of the approach, not a threshold to tune: it would take archaic genomes closer to +//! the populations that introgressed into East Asia, which do not exist. **A cross-population +//! comparable number is therefore not currently achievable this way** — the caller is defensible +//! within a population and not between them. +//! +//! **Still not enough to re-enable.** Beyond the ordering: precision is 34.9 % unfiltered on +//! held-out Europeans, the cohort is **chr21+22 only**, and the reference callset is itself weakly +//! supported (hmmix's own tracts are enriched just 1.84x for their own archaic SNPs), so agreement +//! with it caps well below 100 % even for a correct caller — F1 alone cannot say when this is done. use std::collections::BTreeMap; diff --git a/scripts/archaic-validation/concordance_filter.py b/scripts/archaic-validation/concordance_filter.py new file mode 100644 index 0000000..5652b39 --- /dev/null +++ b/scripts/archaic-validation/concordance_filter.py @@ -0,0 +1,195 @@ +"""Filter called segments by archaic-genome concordance — fitted on the three Neanderthals, +validated on Denisova, which the filter never sees. + +Why: the reported extent orders the populations backwards, and the decomposition says why. The +true-positive component reproduces the truth ordering almost exactly (1.219 against 1.217); the +false positives, which are about twice the true positives, run the other way (0.877). The headline +number is dominated by noise whose population ordering is inverted. So the fix is precision. + +The arbiter discriminates -- false positives score 81.3 %/72.9 % against true positives' ~93.6 % -- +and it uses information the caller does not: which archaic genome carries what. Using it as a filter +is therefore a real gain rather than re-tuning something already fitted. + +The cost is that filtering on the arbiter spends it as an independent referee. Holding out Denisova +keeps one: the filter sees only Altai, Vindija and Chagyrskaya, so Denisova concordance remains an +untouched check on whether the kept segments are genuinely archaic rather than merely +Neanderthal-shaped. +""" + +import bisect +import collections +import json +import os +import subprocess + +SC = os.path.dirname(os.path.abspath(__file__)) +V2 = f'{SC}/v2' +DB = os.path.expanduser('~/.decodingus/navigator-rs.db') +CONTIGS = ('chr21', 'chr22') +RATIO, MIN_POST, MIN_SITES, MIN_BP = '4.5', 0.98, 16, 5_000 +NEANDERTHALS = ['AltaiNeanderthal', 'Vindija33.19', 'Chagyrskaya8'] +HELD_OUT = 'Denisova3' + + +def merge(iv, tol=0): + iv = sorted(iv) + out, cs, ce = [], *iv[0] + for s, e in iv[1:]: + if s > ce + tol: + out.append((cs, ce)) + cs, ce = s, e + else: + ce = max(ce, e) + out.append((cs, ce)) + return out + + +def load_bed(path, tol=0): + d = collections.defaultdict(list) + for line in open(path): + p = line.split() + if len(p) >= 3: + d[p[0]].append((int(p[1]), int(p[2]))) + return {c: merge(v, tol) for c, v in d.items()} + + +def sql(q): + return subprocess.run(['sqlite3', DB, q], capture_output=True, text=True).stdout.strip() + + +def carries_for(sample): + aln = sql(f"SELECT a.id FROM biosample b JOIN sequence_run r ON r.biosample_guid=b.guid " + f"JOIN alignment a ON a.sequence_run_id=r.id WHERE b.donor_identifier='{sample}' " + f"AND a.reference_build='chm13v2.0' LIMIT 1;") + if not aln: + return None + out = {} + for c in CONTIGS: + payload = sql(f"SELECT payload FROM analysis_artifact WHERE alignment_id={aln} " + f"AND kind='diploid_denovo:{c}';") + if not payload: + return None + for rec in json.loads(payload): + d = rec['dosage'] + if d is None or d < 0: + continue + al = set() + if d < 2: + al.add(rec['reference_allele']) + if d > 0: + al.add(rec['alternate_allele']) + out[(rec['contig'], rec['position'])] = al + return out + + +def load_panel(): + by_contig = collections.defaultdict(list) + with open(f'{SC}/panel.tsv') as f: + hdr = next(f).rstrip('\n').split('\t') + ni = [hdr.index(g) for g in NEANDERTHALS] + di = hdr.index(HELD_OUT) + for line in f: + p = line.rstrip('\n').split('\t') + by_contig[p[0]].append((int(p[1]), p[2], [p[i] for i in ni], p[di])) + return {c: sorted(v) for c, v in by_contig.items()} + + +def overlaps(regs, s, e): + for a, b in regs: + if min(e, b) > max(s, a): + return True + return False + + +def scores(panel, contig, s, e, carries): + """(neanderthal concordance, denisova concordance) over the segment; None when too few sites.""" + v = panel.get(contig, []) + keys = [x[0] for x in v] + lo, hi = bisect.bisect_left(keys, s), bisect.bisect_left(keys, e) + nh = [0] * len(NEANDERTHALS) + nd = [0] * len(NEANDERTHALS) + dh = dd = 0 + for pos, derived, ncalls, dcall in v[lo:hi]: + got = carries.get((contig, pos)) + has = bool(got) and derived in got + for i, c in enumerate(ncalls): + if c == 'D': + nd[i] += 1 + nh[i] += has + if dcall == 'D': + dd += 1 + dh += has + best = max((h / d for h, d in zip(nh, nd) if d >= 3), default=None) + den = (dh / dd) if dd >= 3 else None + return best, den + + +def main(): + panel = load_panel() + rows = [] + for grp, listfile, source, n in (('EUR', 'eur60.txt', f'{V2}/sweep', 20), + ('EAS', 'eas30.txt', f'{V2}/sweep_eas', 20)): + for s in [x.strip() for x in open(f'{SC}/{listfile}') if x.strip()][:n]: + f_ = f'{source}/{s}.json' + tp = f'{V2}/truth_{s}.bed' + if not (os.path.exists(f_) and os.path.exists(tp)): + continue + carries = carries_for(s) + if carries is None: + continue + truth = load_bed(tp, tol=1000) + doc = json.load(open(f_)) + doc = doc[RATIO] if RATIO in doc else doc + for seg in doc['segments']: + c = seg['contig'] + if c not in CONTIGS or seg['posterior'] < MIN_POST or seg['n_private'] < MIN_SITES: + continue + if seg['end'] - seg['start'] < MIN_BP: + continue + nea, den = scores(panel, c, seg['start'], seg['end'], carries) + rows.append({ + 'grp': grp, 'sample': s, 'contig': c, 's': seg['start'], 'e': seg['end'], + 'mb': (seg['end'] - seg['start']) / 1e6, + 'tp': overlaps(truth.get(c, []), seg['start'], seg['end']), + 'nea': nea, 'den': den, + }) + print(f'segments scored: {len(rows)} ' + f'(with a Neanderthal score: {sum(1 for r in rows if r["nea"] is not None)})\n') + + truth_mb = {} + for grp, listfile in (('EUR', 'eur60.txt'), ('EAS', 'eas30.txt')): + tot = n = 0 + for s in [x.strip() for x in open(f'{SC}/{listfile}') if x.strip()][:20]: + p = f'{V2}/truth_{s}.bed' + if os.path.exists(p): + t = load_bed(p, tol=1000) + tot += sum(e - a for v in t.values() for a, e in v) / 1e6 + n += 1 + truth_mb[grp] = tot / max(n, 1) + print(f'truth Mb/person: EUR {truth_mb["EUR"]:.3f} EAS {truth_mb["EAS"]:.3f} ' + f'-> ordering to reproduce {truth_mb["EAS"] / truth_mb["EUR"]:.3f}\n') + + npeople = {g: len({r['sample'] for r in rows if r['grp'] == g}) for g in ('EUR', 'EAS')} + print(f"{'threshold':>9s} {'kept':>6s} {'prec':>6s} {'EUR Mb':>7s} {'EAS Mb':>7s} " + f"{'ordering':>9s} {'Denisova kept/dropped':>22s}") + for thr in (0.0, 0.70, 0.80, 0.85, 0.90, 0.95, 1.0): + kept = [r for r in rows if r['nea'] is not None and r['nea'] >= thr] + dropped = [r for r in rows if r['nea'] is not None and r['nea'] < thr] + if not kept: + continue + prec = sum(1 for r in kept if r['tp']) / len(kept) * 100 + eur = sum(r['mb'] for r in kept if r['grp'] == 'EUR') / npeople['EUR'] + eas = sum(r['mb'] for r in kept if r['grp'] == 'EAS') / npeople['EAS'] + dk = [r['den'] for r in kept if r['den'] is not None] + dd = [r['den'] for r in dropped if r['den'] is not None] + dtxt = (f'{sum(dk) / len(dk) * 100:8.1f}% / ' + f'{sum(dd) / len(dd) * 100:6.1f}%') if dk and dd else f'{sum(dk) / len(dk) * 100:8.1f}% / -' + print(f'{thr:9.2f} {len(kept):6d} {prec:5.1f}% {eur:7.3f} {eas:7.3f} ' + f'{eas / eur if eur else 0:9.3f} {dtxt:>22s}') + print('\n ordering: 1.217 is the target. Denisova is the HELD-OUT check -- it never enters the') + print(' filter, so kept segments scoring higher on it than dropped ones means the filter is') + print(' selecting genuinely archaic sequence, not just Neanderthal-shaped noise.') + + +if __name__ == '__main__': + main() diff --git a/scripts/archaic-validation/observable/ordering_cause.py b/scripts/archaic-validation/observable/ordering_cause.py new file mode 100644 index 0000000..71f69e5 --- /dev/null +++ b/scripts/archaic-validation/observable/ordering_cause.py @@ -0,0 +1,169 @@ +"""Does the population difference in BACKGROUND archaic carriage explain the ordering inversion? + +The blocker: reported extent orders the populations backwards (we call 0.937x where the truth says +1.217x), because the false-positive load differs (precision 32.2% EUR vs 41.9% EAS) at identical +thresholds and identical in-tract contrast. + +The arbiter turned up a candidate: outside tracts, Europeans match archaic genomes at 59.0% against +East Asians' 45.5%. If Europeans carry archaic-derived alleles more often in non-introgressed +sequence, then a fixed threshold turns more European background into calls -- more false positives, +more spurious extent, and an inflated European total. + +The model's own background parameter should absorb this: p_background is estimated per individual +from the genome-wide carrying rate. So the question is whether that estimate actually tracks the +difference, or whether the two rates measure different things. + +Measured here, per population: + * the carrying rate the MODEL estimates (all diagnostic sites, what p_background sees) + * the carrying rate the ARBITER sees outside tracts (per-genome concordance) + * the resulting false-positive Mb per person +If the model's estimate does NOT differ while the arbiter's does, p_background is blind to exactly +the thing driving the inversion, and that is the defect to fix. +""" + +import bisect +import collections +import json +import os +import subprocess + +SC = os.path.dirname(os.path.abspath(__file__)) +V2 = f'{SC}/v2' +DB = os.path.expanduser('~/.decodingus/navigator-rs.db') +CONTIGS = ('chr21', 'chr22') +RATIO, MIN_POST, MIN_SITES, MIN_BP = '4.5', 0.98, 16, 5_000 +GENOMES = ['AltaiNeanderthal', 'Vindija33.19', 'Chagyrskaya8', 'Denisova3'] + + +def merge(iv, tol=0): + iv = sorted(iv) + out, cs, ce = [], *iv[0] + for s, e in iv[1:]: + if s > ce + tol: + out.append((cs, ce)) + cs, ce = s, e + else: + ce = max(ce, e) + out.append((cs, ce)) + return out + + +def load_bed(path, tol=0): + d = collections.defaultdict(list) + for line in open(path): + p = line.split() + if len(p) >= 3: + d[p[0]].append((int(p[1]), int(p[2]))) + return {c: merge(v, tol) for c, v in d.items()} + + +def sql(q): + return subprocess.run(['sqlite3', DB, q], capture_output=True, text=True).stdout.strip() + + +def carries_for(sample): + aln = sql(f"SELECT a.id FROM biosample b JOIN sequence_run r ON r.biosample_guid=b.guid " + f"JOIN alignment a ON a.sequence_run_id=r.id WHERE b.donor_identifier='{sample}' " + f"AND a.reference_build='chm13v2.0' LIMIT 1;") + if not aln: + return None + out = {} + for c in CONTIGS: + payload = sql(f"SELECT payload FROM analysis_artifact WHERE alignment_id={aln} " + f"AND kind='diploid_denovo:{c}';") + if not payload: + return None + for rec in json.loads(payload): + d = rec['dosage'] + if d is None or d < 0: + continue + al = set() + if d < 2: + al.add(rec['reference_allele']) + if d > 0: + al.add(rec['alternate_allele']) + out[(rec['contig'], rec['position'])] = al + return out + + +def overlaps(regs, s, e): + for a, b in regs: + if min(e, b) > max(s, a): + return True + return False + + +def in_reg(regs, p): + i = bisect.bisect_right([s for s, _ in regs], p) - 1 + return i >= 0 and regs[i][1] > p + + +def main(): + # every diagnostic site (what the MODEL's p_background is estimated over) + diag = collections.defaultdict(dict) + with open(f'{SC}/classify.tsv') as f: + next(f) + for line in f: + c, p, d, _ = line.rstrip('\n').split('\t') + diag[c][int(p)] = d + callable_r = load_bed(f'{SC}/callable.bed') + + print(f"{'pop':4s} {'model p_bg (all sites)':>24s} {'model p_bg OUTSIDE tracts':>27s} " + f"{'FP Mb/person':>13s}") + out = {} + for grp, listfile, source, n in (('EUR', 'eur60.txt', f'{V2}/sweep', 15), + ('EAS', 'eas30.txt', f'{V2}/sweep_eas', 15)): + all_hit = all_tot = out_hit = out_tot = 0 + fp_mb = 0.0 + people = 0 + for s in [x.strip() for x in open(f'{SC}/{listfile}') if x.strip()][:n]: + f_ = f'{source}/{s}.json' + tp = f'{V2}/truth_{s}.bed' + if not (os.path.exists(f_) and os.path.exists(tp)): + continue + carries = carries_for(s) + if carries is None: + continue + people += 1 + truth = load_bed(tp, tol=1000) + for c in CONTIGS: + cal = callable_r.get(c, []) + tr = truth.get(c, []) + for p, derived in diag[c].items(): + if not cal or not in_reg(cal, p): + continue + got = carries.get((c, p)) + hit = bool(got) and derived in got + all_tot += 1 + all_hit += hit + if not (tr and in_reg(tr, p)): + out_tot += 1 + out_hit += hit + doc = json.load(open(f_)) + doc = doc[RATIO] if RATIO in doc else doc + for seg in doc['segments']: + c = seg['contig'] + if c not in CONTIGS or seg['posterior'] < MIN_POST or seg['n_private'] < MIN_SITES: + continue + if seg['end'] - seg['start'] < MIN_BP: + continue + if not overlaps(truth.get(c, []), seg['start'], seg['end']): + fp_mb += (seg['end'] - seg['start']) / 1e6 + a = all_hit / all_tot * 100 + o = out_hit / out_tot * 100 + print(f'{grp:4s} {a:23.2f}% {o:26.2f}% {fp_mb / people:13.3f}') + out[grp] = (a, o, fp_mb / people) + + print() + e, a = out['EUR'], out['EAS'] + print(f' model p_background EUR {e[0]:.2f}% vs EAS {a[0]:.2f}% ratio {a[0] / e[0]:.3f}') + print(f' outside tracts EUR {e[1]:.2f}% vs EAS {a[1]:.2f}% ratio {a[1] / e[1]:.3f}') + print(f' false-positive Mb EUR {e[2]:.3f} vs EAS {a[2]:.3f} ratio {a[2] / e[2]:.3f}') + print() + print(' If the model estimate is flat across populations while the FP load is not, then') + print(' p_background is blind to whatever drives the extra European calls, and no amount of') + print(' per-individual estimation of it will fix the ordering.') + + +if __name__ == '__main__': + main() From b39db0b620b122dbd12c2c7823539986c02e7536 Mon Sep 17 00:00:00 2001 From: James Kane Date: Fri, 31 Jul 2026 13:31:22 -0500 Subject: [PATCH 8/9] Ship Tier B archaic segments as a WITHIN-POPULATION measure Re-enables ARCHAIC_SEGMENTS_ENABLED on the rebuilt caller, wires it through the app, and states the comparability limit where the number is read. WHAT EARNED IT BACK. Held out on 30 Europeans the fit never saw, per-individual extent correlates at r = +0.710 (p < 0.0001) where the withdrawn density caller managed -0.018 (p = 0.94). All 90 individuals across two populations score above their own random-placement null. The concordance filter takes precision from 54% to 90%, verified with Denisova held out of it entirely: kept segments score 74.9% on Denisova concordance against 21.5% for dropped ones. WHAT IT MUST NOT BE USED FOR, and why the caveat is in the UI rather than a doc: East Asian tracts match our four sequenced archaic genomes less well than European ones (83.4% vs 89.2%), so their extent is under-called and the reported figure orders the two populations backwards against the truth. That is a property of which archaic genomes have been sequenced -- Denisova is the best match for 32.2% of East Asian tracts against 11.2% of European ones -- not a threshold to tune. The caveat sits directly under the number, in amber, because a limit recorded anywhere else will not be read by the person reading the number. CHANGES: - archaic_match gains segment_concordance / carried_panel_sites / filter_by_concordance, with MIN_CONCORDANCE = 0.70 from the measured plateau. - The app method now runs matching + filter instead of the density caller, and loads the Tier A panel for the per-genome calls. The outgroup asset is no longer needed here: it existed to strip African-shared variants for a density model that no longer exists. - The reference is read one contig at a time -- whether a diagnostic site is informative depends on the reference base there, and holding all of CHM13 would cost 3.1 GB. - archaic_segment_sig now includes METHOD_VERSION. Without it a workspace holding results from the withdrawn caller would keep serving them: stale answers from a method removed for having no per-individual signal. Three of the eight tests pin failures that actually happened during development: concordance must condition on the GENOME not the subject (the other way scores ~100% everywhere including background, separating nothing); a segment with too few judgable sites is KEPT, since absence of evidence is not evidence of a bad call and dropping on it would quietly penalise sparse regions; and the summary is recomputed after filtering rather than carried over. VALIDATED ON chr21+chr22 ONLY. This ships genome-wide. The design records the same extrapolation trap from the previous caller -- projecting the chr21+22 target predicted 85.1 Mb against a measured 90.9, 6% low -- so this is a real gap, not a formality. A genome-wide run is the next thing to do, and the inputs for it (whole-genome calls on the ground-truth subject, hmmix's genome-wide European callset) are already local. Verified: clippy clean under -D warnings; 30 suites pass including i18n parity; `navigator archaic-segments` runs end to end through the app path. Co-Authored-By: Claude Opus 5 (1M context) --- .../navigator-analysis/src/archaic_match.rs | 274 +++++++++++++++++- crates/navigator-app/src/haplogroup.rs | 49 +++- crates/navigator-app/src/lib.rs | 84 ++++-- crates/navigator-domain/locales/en.txt | 1 + crates/navigator-domain/locales/es.txt | 1 + crates/navigator-ui/src/ui/central.rs | 10 + 6 files changed, 378 insertions(+), 41 deletions(-) diff --git a/crates/navigator-analysis/src/archaic_match.rs b/crates/navigator-analysis/src/archaic_match.rs index 6d70fc7..3ce9202 100644 --- a/crates/navigator-analysis/src/archaic_match.rs +++ b/crates/navigator-analysis/src/archaic_match.rs @@ -146,11 +146,26 @@ use std::collections::BTreeMap; use serde::{Deserialize, Serialize}; -use crate::archaic::{ArchaicCallable, ArchaicClassify, DiagnosticClass}; +use crate::archaic::{ArchaicCallable, ArchaicClassify, ArchaicMarkerPanel, DiagnosticClass, ARCHAIC_GENOMES}; use crate::archaic_segments::{ArchaicSegment, ArchaicSegmentResult, ArchaicSource, ArchaicSummary}; use crate::caller::SiteGenotype; use crate::ibd::GeneticMap; +/// Bumped whenever a change would alter the segments this module produces. +/// +/// Persisted results are keyed on it (`archaic_segment_sig`), so a workspace holding output from an +/// earlier method — notably the withdrawn private-variant density caller — re-derives instead of +/// serving answers the current code would never produce. +pub const METHOD_VERSION: u32 = 1; + +/// Concordance a segment must reach to be kept — measured, from the threshold sweep where precision +/// plateaus (54 % -> 90 % at 0.70, and no better above it while recall keeps falling). +pub const MIN_CONCORDANCE: f64 = 0.70; + +/// Sites a genome needs in a segment before its concordance is trusted. Without a floor a genome +/// called at a single site scores 1.0 and wins every segment. +pub const MIN_CONCORDANCE_SITES: usize = 3; + /// One diagnostic site, reduced to what the HMM consumes. #[derive(Debug, Clone, Copy, PartialEq)] pub struct SiteObs { @@ -454,10 +469,115 @@ pub fn call_from_observations( ArchaicSegmentResult { segments, summary } } +/// How well a segment matches each archaic genome: of the sites where a given archaic genome +/// carries the derived allele, what fraction does the subject also carry. The best genome wins. +/// +/// Conditioning on the **genome**, not on the subject, is what makes this discriminate. The +/// intuitive version — over the sites the subject carries, how many does an archaic genome share — +/// scores ~100 % everywhere including background, because at an informative site some archaic +/// carries the derived allele by construction. Read this way, background sits at the subject's +/// genome-wide carrying rate and an inherited haplotype sits far above it. +/// +/// Returns `None` when no genome has enough called sites in the span to judge. +pub fn segment_concordance( + panel: &ArchaicMarkerPanel, + contig: &str, + start: i64, + end: i64, + carried: &BTreeMap<(&str, i64), bool>, + min_sites: usize, +) -> Option { + let mut hits = [0usize; ARCHAIC_GENOMES.len()]; + let mut dens = [0usize; ARCHAIC_GENOMES.len()]; + for s in &panel.sites { + if s.contig != contig || s.position < start || s.position > end { + continue; + } + let subject_has = carried.get(&(contig, s.position)).copied().unwrap_or(false); + for (i, call) in s.calls.iter().enumerate() { + if call.carries_derived() { + dens[i] += 1; + hits[i] += usize::from(subject_has); + } + } + } + (0..ARCHAIC_GENOMES.len()) + .filter(|&i| dens[i] >= min_sites) + .map(|i| hits[i] as f64 / dens[i] as f64) + .fold(None, |best: Option, r| Some(best.map_or(r, |b| b.max(r)))) +} + +/// Which sites of the Tier A panel the subject carries the archaic-derived allele at. +/// +/// A site with no variant record is hom-reference and therefore **not** a carrier; sites where the +/// derived allele is the reference base are excluded upstream by the panel's own orientation. +pub fn carried_panel_sites<'a>( + panel: &'a ArchaicMarkerPanel, + calls: &'a [SiteGenotype], +) -> BTreeMap<(&'a str, i64), bool> { + let by_pos: BTreeMap<(&str, i64), &SiteGenotype> = + calls.iter().map(|c| ((c.contig.as_str(), c.position), c)).collect(); + let mut out = BTreeMap::new(); + for s in &panel.sites { + let Some((k, g)) = by_pos.get_key_value(&(s.contig.as_str(), s.position)) else { continue }; + let carries = g.dosage > 0 && g.alternate_allele.starts_with(s.archaic_derived_allele); + out.insert(*k, carries); + } + out +} + +/// Drop segments that do not look like an inherited archaic haplotype. +/// +/// This is the single largest quality lever measured: **precision 54 % -> 90 %** on real data. It +/// works because it consults evidence the segment caller never sees — which archaic genome carries +/// what — so it is genuinely new information rather than a re-reading of the same signal. +/// +/// A segment with too few judgable sites is **kept**: absence of evidence is not evidence of a bad +/// call, and dropping on it would silently penalise sparse regions. +/// +/// Note the population caveat in the module docs: East Asian tracts match the four sequenced archaic +/// genomes less well than European ones, so this filter removes proportionally more of them. It +/// improves precision everywhere and makes extent **less** comparable between populations. +pub fn filter_by_concordance( + result: ArchaicSegmentResult, + panel: &ArchaicMarkerPanel, + calls: &[SiteGenotype], + min_concordance: f64, + min_sites: usize, +) -> ArchaicSegmentResult { + let carried = carried_panel_sites(panel, calls); + let kept: Vec = result + .segments + .into_iter() + .filter(|seg| { + match segment_concordance(panel, &seg.contig, seg.start, seg.end, &carried, min_sites) { + Some(c) => c >= min_concordance, + None => true, + } + }) + .collect(); + let total_mb: f64 = kept.iter().map(|s| s.length_mb()).sum(); + let callable_mb = result.summary.callable_mb; + ArchaicSegmentResult { + summary: ArchaicSummary { + total_mb, + pct_callable: if callable_mb > 0.0 { total_mb * 100.0 / callable_mb } else { 0.0 }, + callable_mb, + neanderthal_mb: 0.0, + denisovan_mb: 0.0, + unknown_mb: total_mb, + n_segments: kept.len(), + }, + segments: kept, + } +} + #[cfg(test)] mod tests { use super::*; - use crate::archaic::{CallableContig, ClassifyContig, PositionStream}; + use crate::archaic::{ + ArchaicCall, ArchaicPanelThresholds, ArchaicSite, CallableContig, ClassifyContig, PositionStream, + }; fn gmap(contig: &str, len: i32) -> GeneticMap { GeneticMap::uniform(1.0, &[(contig, len)]) @@ -583,6 +703,156 @@ mod tests { assert!(out.iter().all(|o| !o.carries), "no calls means nothing carried"); } + fn panel_site(pos: i64, derived: char, calls: [ArchaicCall; 4]) -> ArchaicSite { + ArchaicSite { + contig: "chr21".into(), + position: pos, + reference_allele: 'T', + alternate_allele: derived, + archaic_derived_allele: derived, + calls, + diagnostic_class: DiagnosticClass::SharedArchaic, + afr_freq: 0.0, + grch37: None, + grch38: None, + } + } + + fn genotype(pos: i64, alt: char, dosage: i32) -> SiteGenotype { + SiteGenotype { + name: String::new(), + contig: "chr21".into(), + position: pos, + reference_allele: "T".into(), + alternate_allele: alt.to_string(), + ploidy: 2, + dosage, + gq: 60, + depth: 30, + ref_depth: 15, + alt_depth: 15, + pls: Vec::new(), + gt: None, + allele_depths: None, + } + } + + /// Concordance must be read per ARCHAIC GENOME, not per carried site. The intuitive version — + /// over the sites the subject carries, how many does some archaic share — scores ~100 % + /// everywhere including background, because at an informative site some archaic is derived by + /// construction. That version was written first and separated nothing. + #[test] + fn concordance_conditions_on_the_genome_not_the_subject() { + use ArchaicCall::{HomAncestral as A, HomDerived as D}; + // Altai is derived at all 4 sites; Denisova at only the first. + let panel = ArchaicMarkerPanel { + build: "chm13v2.0".into(), + thresholds: ArchaicPanelThresholds { max_afr_freq: 0.01, min_non_afr_freq: 0.0005 }, + sites: vec![ + panel_site(1_000, 'A', [D, A, A, D]), + panel_site(2_000, 'A', [D, A, A, A]), + panel_site(3_000, 'A', [D, A, A, A]), + panel_site(4_000, 'A', [D, A, A, A]), + ], + }; + // The subject carries 3 of Altai's 4 → 0.75 for Altai, 1.0 for Denisova but on ONE site. + let calls = vec![ + genotype(1_000, 'A', 1), + genotype(2_000, 'A', 1), + genotype(3_000, 'A', 1), + ]; + let carried = carried_panel_sites(&panel, &calls); + // min_sites = 3 excludes Denisova's single site, so Altai's 0.75 is the answer. Without + // that floor a 1/1 genome would win every segment. + let c = segment_concordance(&panel, "chr21", 0, 5_000, &carried, 3).expect("a score"); + assert!((c - 0.75).abs() < 1e-9, "expected Altai's 0.75, got {c}"); + } + + /// A segment with too few judgable sites must be KEPT. Absence of evidence is not evidence of a + /// bad call, and dropping on it would quietly penalise sparse regions — which are exactly the + /// regions where a caller most needs the benefit of the doubt. + #[test] + fn filter_keeps_segments_it_cannot_judge() { + let panel = ArchaicMarkerPanel { + build: "chm13v2.0".into(), + thresholds: ArchaicPanelThresholds { max_afr_freq: 0.01, min_non_afr_freq: 0.0005 }, + sites: vec![panel_site(1_000, 'A', [ArchaicCall::HomDerived; 4])], + }; + let seg = ArchaicSegment { + contig: "chr21".into(), + start: 500_000, + end: 600_000, // no panel sites here at all + posterior: 0.99, + n_private: 40, + source: ArchaicSource::Unknown, + neanderthal_matches: 0, + denisovan_matches: 0, + }; + let r = ArchaicSegmentResult { + summary: ArchaicSummary { + total_mb: 0.1, + pct_callable: 1.0, + callable_mb: 10.0, + neanderthal_mb: 0.0, + denisovan_mb: 0.0, + unknown_mb: 0.1, + n_segments: 1, + }, + segments: vec![seg], + }; + let out = filter_by_concordance(r, &panel, &[], 0.9, 3); + assert_eq!(out.segments.len(), 1, "an unjudgable segment must survive"); + } + + /// The filter's whole purpose: a segment that does not look like an inherited archaic haplotype + /// goes, one that does stays. This is the 54 % -> 90 % precision lever. + #[test] + fn filter_drops_poorly_matching_segments() { + use ArchaicCall::{HomAncestral as A, HomDerived as D}; + let mut sites = Vec::new(); + for i in 0..10 { + sites.push(panel_site(1_000 + i * 100, 'A', [D, A, A, A])); // good segment + } + for i in 0..10 { + sites.push(panel_site(50_000 + i * 100, 'A', [D, A, A, A])); // bad segment + } + let panel = ArchaicMarkerPanel { + build: "chm13v2.0".into(), + thresholds: ArchaicPanelThresholds { max_afr_freq: 0.01, min_non_afr_freq: 0.0005 }, + sites, + }; + // Carries 9/10 in the first span, 1/10 in the second. + let mut calls: Vec = (0..9).map(|i| genotype(1_000 + i * 100, 'A', 1)).collect(); + calls.push(genotype(50_000, 'A', 1)); + + let mk = |start: i64, end: i64| ArchaicSegment { + contig: "chr21".into(), + start, + end, + posterior: 0.99, + n_private: 20, + source: ArchaicSource::Unknown, + neanderthal_matches: 0, + denisovan_matches: 0, + }; + let r = ArchaicSegmentResult { + summary: ArchaicSummary { + total_mb: 0.002, + pct_callable: 1.0, + callable_mb: 10.0, + neanderthal_mb: 0.0, + denisovan_mb: 0.0, + unknown_mb: 0.002, + n_segments: 2, + }, + segments: vec![mk(900, 2_000), mk(49_900, 51_000)], + }; + let out = filter_by_concordance(r, &panel, &calls, 0.7, 3); + assert_eq!(out.segments.len(), 1, "the poorly-matching segment should go"); + assert_eq!(out.segments[0].start, 900); + assert_eq!(out.summary.n_segments, 1, "the summary must be recomputed, not carried over"); + } + /// A no-call is hom-reference, i.e. NOT carrying. Conditioning on "has a call" instead is what /// made an early version of this analysis report ~80 % carrying against a 4.3 % background. #[test] diff --git a/crates/navigator-app/src/haplogroup.rs b/crates/navigator-app/src/haplogroup.rs index c2d9c3a..03d0cae 100644 --- a/crates/navigator-app/src/haplogroup.rs +++ b/crates/navigator-app/src/haplogroup.rs @@ -3612,19 +3612,24 @@ impl App { } let rb = ReferenceBuild::Chm13v2; + // The outgroup track is no longer needed: the density caller used it to strip + // African-shared variants, and this method does not model density at all. for p in [ - crate::archaic_outgroup_path(rb), crate::archaic_classify_path(rb), crate::archaic_callable_path(rb), + crate::archaic_markers_path(rb), ] { self.ensure_ancestry_asset(rb, &p).await?; } let load = |p: PathBuf| -> Result, AppError> { crate::read_verified_asset(rb, &p)?.ok_or_else(|| AppError::AncestryPanelMissing(p.clone())) }; - let outgroup = ArchaicOutgroup::from_bytes(&load(crate::archaic_outgroup_path(rb))?)?; let classify = ArchaicClassify::from_bytes(&load(crate::archaic_classify_path(rb))?)?; let callable = ArchaicCallable::from_bytes(&load(crate::archaic_callable_path(rb))?)?; + // Tier A's panel carries the per-archaic-genome calls the concordance filter needs — the + // single largest quality lever measured (precision 54 % -> 90 %). + let panel = ArchaicMarkerPanel::from_bytes(&load(crate::archaic_markers_path(rb))?)?; + let (_, reference) = self.alignment_bam_reference(aln).await?; // Genetic map over the contigs actually present, so transitions are recombination-scaled. let mut lengths: std::collections::BTreeMap = std::collections::BTreeMap::new(); @@ -3636,17 +3641,39 @@ impl App { let gmap = crate::load_genetic_map(rb, &pairs); eprintln!("archaic segments: {} calls over {contigs_present} autosome(s)", calls.len()); - let result = tokio::task::spawn_blocking(move || { - navigator_analysis::archaic_segments::call_archaic_segments( + let result = tokio::task::spawn_blocking(move || -> Result<_, AppError> { + use navigator_analysis::archaic_match as am; + + let mut by_contig: std::collections::BTreeMap> = + Default::default(); + for c in &calls { + by_contig.entry(c.contig.clone()).or_default().insert(c.position, c); + } + let mut observations = std::collections::BTreeMap::new(); + for (contig, pos_map) in &by_contig { + // One contig's reference at a time: whether a diagnostic site is informative depends + // on the reference base there, and holding all of CHM13 would cost 3.1 GB. + let seq = navigator_analysis::reader::read_contig_sequence(&reference, contig)?; + let obs = am::observations_for_contig( + contig, + &classify, + pos_map, + |p| seq.get((p - 1).max(0) as usize).copied().map(|b| b.to_ascii_uppercase()), + &callable, + am::MatchConfig::default().min_callable_fraction, + ); + observations.insert(contig.clone(), obs); + } + let called = am::call_from_observations(&observations, &gmap, &callable, &am::MatchConfig::default()); + Ok(am::filter_by_concordance( + called, + &panel, &calls, - &outgroup, - &classify, - &callable, - &gmap, - &navigator_analysis::archaic_segments::ArchaicConfig::default(), - ) + am::MIN_CONCORDANCE, + am::MIN_CONCORDANCE_SITES, + )) }) - .await?; + .await??; consensus_archaic_segments::upsert( self.store.pool(), diff --git a/crates/navigator-app/src/lib.rs b/crates/navigator-app/src/lib.rs index 6364c84..7416831 100644 --- a/crates/navigator-app/src/lib.rs +++ b/crates/navigator-app/src/lib.rs @@ -1123,7 +1123,15 @@ fn archaic_marker_dist_path(build: ReferenceBuild) -> PathBuf { /// Cache signature for a Tier B segment result: the alignment it came from plus the caller's /// genotype version, so re-calling with a newer caller invalidates it. pub(crate) fn archaic_segment_sig(alignment_id: i64) -> String { - format!("aln{alignment_id}:gt{}", navigator_analysis::caller::GENOTYPE_VERSION) + // The METHOD version is part of the key, not just the genotype version. Tier B was rebuilt from + // a private-variant density model to archaic-genome matching, and without this a workspace + // carrying results from the old caller would keep serving them — stale answers from a method + // that was withdrawn for having no per-individual signal. + format!( + "aln{alignment_id}:gt{}:m{}", + navigator_analysis::caller::GENOTYPE_VERSION, + navigator_analysis::archaic_match::METHOD_VERSION + ) } fn archaic_outgroup_path(build: ReferenceBuild) -> PathBuf { @@ -2451,36 +2459,56 @@ pub struct AncientFitRow { pub const ANCIENT_ANCESTRY_ENABLED: bool = true; /// Whether Tier B **archaic segments** (the introgressed-tract caller and its chromosome browser) -/// are computed, read back, or shown. **Off**: validation against an external per-individual truth -/// set showed the caller carries no per-person signal. -/// -/// Tier B was shipped on the strength of one number — its total extent landed at 1.01x the hmmix -/// European mean. Validating it properly, against hmmix's own calls **for the same individuals** -/// (n=20 Europeans, chr21+22), showed that number is all there is: -/// -/// - **Locations disagree.** For HG00096, 2.1 % of hmmix's archaic bases are also called by us, -/// against a 5.0 % expectation (p95 9.4 %) for segments of our own lengths placed at *random* in -/// the same span. Below chance. Not a coordinate artefact: the overlap-vs-shift curve is flat -/// across +/-2 Mb with no peak, and 70.7 % of the truth lies inside our callable territory, so -/// the tracts were reachable. -/// - **Amounts do not track the individual.** Across the 20, Pearson r = -0.018 (p = 0.94) and -/// Spearman rho = -0.020 against a true range of 1.19-2.97 Mb. Our own spread is 0.63x the -/// truth's. The two individuals with the least archaic ancestry drew our two highest calls. -/// -/// The mean ratio really is ~0.92 — the caller reproduces the cohort average and nothing else, -/// which is what three fitted parameters were tuned to do. An honest report needs a measurement of -/// *this person*, so the feature is withheld rather than shown with a caveat. +/// are computed, read back, or shown. /// -/// The machinery stays, unit-tested, behind this flag — same discipline as `attribute_lineage` and -/// the ancient-ancestry precedent above. Re-enabling needs a method change (the design records -/// Skov-2020 haplotype matching as the path), not a threshold tweak, and a re-run of the validation -/// harness that produced these numbers. +/// **ON**, for the rebuilt caller ([`navigator_analysis::archaic_match`]) — and specifically as a +/// **within-population** measure. The history below is the first implementation, which was withdrawn. /// -/// Tier A — the marker **count** and percentile — is a different method on a different asset and is -/// **not** gated by this. +/// What changed: the caller no longer counts private-variant density, it matches the archaic genomes +/// directly. Held out on 30 Europeans the fit never saw, per-individual extent correlates at +/// **r = +0.710 (p < 0.0001)** where the old caller managed −0.018 (p = 0.94); every individual of +/// 90 scores above their own random-placement null; and a concordance filter takes precision from +/// 54 % to 90 %. /// -/// See `documents/design/ArchaicAncestry_Design.md`, "Tier B validation (2026-07-30)". -pub const ARCHAIC_SEGMENTS_ENABLED: bool = false; +/// **The one thing it must not be used for is comparing people of different ancestries.** East Asian +/// tracts match our four sequenced archaic genomes less well than European ones, so extent is +/// under-called for them — the reported figure orders the two populations backwards against the +/// truth. That is a property of which archaic genomes have been sequenced, not a threshold, so the +/// UI states the limit rather than implying a universal percentage. +pub const ARCHAIC_SEGMENTS_ENABLED: bool = true; + +// The prior gate's rationale, kept because it is the reason the current caller exists. +// +// **Off**: validation against an external per-individual truth +// set showed the caller carries no per-person signal. +// +// Tier B was shipped on the strength of one number — its total extent landed at 1.01x the hmmix +// European mean. Validating it properly, against hmmix's own calls **for the same individuals** +// (n=20 Europeans, chr21+22), showed that number is all there is: +// +// - **Locations disagree.** For HG00096, 2.1 % of hmmix's archaic bases are also called by us, +// against a 5.0 % expectation (p95 9.4 %) for segments of our own lengths placed at *random* in +// the same span. Below chance. Not a coordinate artefact: the overlap-vs-shift curve is flat +// across +/-2 Mb with no peak, and 70.7 % of the truth lies inside our callable territory, so +// the tracts were reachable. +// - **Amounts do not track the individual.** Across the 20, Pearson r = -0.018 (p = 0.94) and +// Spearman rho = -0.020 against a true range of 1.19-2.97 Mb. Our own spread is 0.63x the +// truth's. The two individuals with the least archaic ancestry drew our two highest calls. +// +// The mean ratio really is ~0.92 — the caller reproduces the cohort average and nothing else, +// which is what three fitted parameters were tuned to do. An honest report needs a measurement of +// *this person*, so the feature is withheld rather than shown with a caveat. +// +// The machinery stays, unit-tested, behind this flag — same discipline as `attribute_lineage` and +// the ancient-ancestry precedent above. Re-enabling needs a method change (the design records +// Skov-2020 haplotype matching as the path), not a threshold tweak, and a re-run of the validation +// harness that produced these numbers. +// +// Tier A — the marker **count** and percentile — is a different method on a different asset and is +// **not** gated by this. +// +// See `documents/design/ArchaicAncestry_Design.md`, "Tier B validation (2026-07-30)" and +// "Why it failed" for the full record. Historical only: the live gate is the `true` above. /// The persisted method name of the deep-ancestry breakdown — re-exported so the UI reads the /// rebuilt method by name and can never fall back to a retired one. diff --git a/crates/navigator-domain/locales/en.txt b/crates/navigator-domain/locales/en.txt index e447eef..47effd7 100644 --- a/crates/navigator-domain/locales/en.txt +++ b/crates/navigator-domain/locales/en.txt @@ -807,3 +807,4 @@ simple.relatives.estimate=Estimated from shared signals — connect to measure h # Tier B archaic segments — withheld pending a method that carries per-person signal. archaicSegments.withheld=Not reported: we could not show you where your archaic DNA sits accurately enough to be worth reporting. archaicSegments.withheldWhy=Checked against an independent published callset for the same 20 people, our tract locations matched no better than chance, and the total did not track the individual. The overall amount for a population came out right, which is not the same as being right about you. We would rather report nothing than a map that looks precise and is not. +archaicSegments.withinPopulation=Compare this only with people of similar ancestry. The amount is measured against four sequenced archaic genomes, and they match some ancestries better than others — so this figure is not a like-for-like number between, say, a European and an East Asian. diff --git a/crates/navigator-domain/locales/es.txt b/crates/navigator-domain/locales/es.txt index 696fd82..a9500bc 100644 --- a/crates/navigator-domain/locales/es.txt +++ b/crates/navigator-domain/locales/es.txt @@ -792,3 +792,4 @@ simple.relatives.estimate=Estimado a partir de señales compartidas: conecta par # Segmentos arcaicos de nivel B: retenidos hasta disponer de un método con señal individual. archaicSegments.withheld=No se informa: no hemos podido mostrarte dónde se sitúa tu ADN arcaico con la precisión necesaria para publicarlo. archaicSegments.withheldWhy=Al contrastarlo con un conjunto de datos publicado e independiente para las mismas 20 personas, la ubicación de nuestros tramos no acertó más que el azar, y el total no seguía a cada individuo. La cantidad global de una población sí salía bien, lo cual no equivale a acertar contigo. Preferimos no informar nada antes que ofrecer un mapa que parece preciso y no lo es. +archaicSegments.withinPopulation=Compara esta cifra solo con personas de ascendencia similar. La cantidad se mide frente a cuatro genomas arcaicos secuenciados, que se ajustan mejor a unas ascendencias que a otras: no es una cifra equiparable entre, por ejemplo, una persona europea y una del este de Asia. diff --git a/crates/navigator-ui/src/ui/central.rs b/crates/navigator-ui/src/ui/central.rs index 5aa4575..f515d98 100644 --- a/crates/navigator-ui/src/ui/central.rs +++ b/crates/navigator-ui/src/ui/central.rs @@ -663,6 +663,16 @@ impl NavigatorApp { "{} archaic tracts, {:.2}% of the {:.0} Mb we could read reliably.", r.summary.n_segments, r.summary.pct_callable, r.summary.callable_mb )); + ui.add_space(6.0); + // The comparability limit sits directly under the number, not in a + // footnote: this figure is measured against four sequenced archaic + // genomes that represent some ancestries better than others, so + // comparing it between people of different ancestry is the one use it + // cannot support. Stated where the number is read, or it will not be. + ui.label( + egui::RichText::new(self.tr("archaicSegments.withinPopulation")) + .color(egui::Color32::from_rgb(230, 180, 90)), + ); ui.add_space(4.0); // Lineage split is withheld, not merely absent — say so. ui.label( From 9fca4c1234e3fde65ca06e7f708b02fbc8582b5d Mon Sep 17 00:00:00 2001 From: James Kane Date: Sat, 1 Aug 2026 05:41:43 -0500 Subject: [PATCH 9/9] Genome-wide validation passes; fix a cache key that served partial results forever Ran the validation the previous commit called for, and it found a real bug before it produced a number. THE BUG. `archaic_segment_sig` keyed on alignment + genotype version + method version, but NOT on which contigs were actually called. So a result computed from a partial call set stayed "current" indefinitely. Caught in the act: HG00096 had all 23 contigs cached and the app cheerfully served 33 segments over 2 contigs (1.94 Mb) from an earlier chr21+22 run. Anyone who called one chromosome, then later ran the whole genome, would have kept the one-chromosome answer. Fixed by folding the called contigs into the key. Reading them needed a cheap query -- `list_for_alignment` returns payloads, which here is ~1 GB of JSON across 22 contigs -- so `artifact::list_kinds` returns kinds alone. With the fix HG00096 recomputes to 1,208 segments over 22 contigs (83.56 Mb). GENOME-WIDE VALIDATION, three Europeans against hmmix's genome-wide callset: ours hmmix ratio sens prec null (max of 400) HG00096 83.6 Mb 93.0 0.90 40.3% 44.9% 5.5% HG00102 83.9 Mb 89.3 0.94 42.4% 45.1% 4.9% HG00112 82.1 Mb 91.0 0.90 42.9% 47.5% 5.1% All three above the ENTIRE random-placement null. Both sensitivity and precision come out BETTER genome-wide than on chr21+22 (40-43% against 31.6%, ~46% against 34.9%), so the two-chromosome numbers were conservative rather than optimistic. That direction matters. The design records the previous caller being burned by the opposite -- a chr21+22 target that extrapolated 6% low -- which is exactly why shipping on a two-chromosome result was worth objecting to even though the result here turned out fine. The objection was right; the answer was just favourable. Extent lands at 0.90-0.94 of hmmix's, i.e. we under-call by about 10% genome-wide rather than over-call. Module docs now lead with these figures. Co-Authored-By: Claude Opus 5 (1M context) --- .../navigator-analysis/src/archaic_match.rs | 18 ++++++++++ crates/navigator-app/src/haplogroup.rs | 7 ++-- crates/navigator-app/src/lib.rs | 36 +++++++++++++++---- crates/navigator-store/src/artifact.rs | 14 ++++++++ 4 files changed, 66 insertions(+), 9 deletions(-) diff --git a/crates/navigator-analysis/src/archaic_match.rs b/crates/navigator-analysis/src/archaic_match.rs index 3ce9202..e3e636d 100644 --- a/crates/navigator-analysis/src/archaic_match.rs +++ b/crates/navigator-analysis/src/archaic_match.rs @@ -38,6 +38,24 @@ //! //! # Validation //! +//! ## Genome-wide (the shipping configuration) +//! +//! Three Europeans called across all 22 autosomes and scored against hmmix's genome-wide callset for +//! the same individuals: +//! +//! | | ours | hmmix | ratio | sensitivity | precision | null (max of 400 draws) | +//! |---|---|---|---|---|---|---| +//! | HG00096 | 83.6 Mb | 93.0 | 0.90 | 40.3 % | 44.9 % | 5.5 % | +//! | HG00102 | 83.9 Mb | 89.3 | 0.94 | 42.4 % | 45.1 % | 4.9 % | +//! | HG00112 | 82.1 Mb | 91.0 | 0.90 | 42.9 % | 47.5 % | 5.1 % | +//! +//! All three sit above the *entire* random-placement null. Both sensitivity and precision are +//! **better** genome-wide than on chr21+22 (40–43 % against 31.6 %, ~46 % against 34.9 %), so the +//! two-chromosome figures below are conservative rather than optimistic — worth stating because the +//! previous caller's design was burned by the opposite, extrapolating a chr21+22 target 6 % low. +//! +//! ## chr21+22, with a train/test split +//! //! Scored against hmmix's own calls for the same individuals, 60 Europeans on chr21+22, split 30 //! **train** / 30 **test** on a fixed seed. Thresholds were fitted on train only; every figure below //! is the held-out half. The split exists because the previous caller was tuned until a cohort diff --git a/crates/navigator-app/src/haplogroup.rs b/crates/navigator-app/src/haplogroup.rs index 03d0cae..d8821df 100644 --- a/crates/navigator-app/src/haplogroup.rs +++ b/crates/navigator-app/src/haplogroup.rs @@ -3532,7 +3532,8 @@ impl App { let Some(aln) = self.alignment_with_diploid_calls(biosample_guid).await? else { return Ok(None); }; - if row.source_sig == archaic_segment_sig(aln) { + let contigs = crate::called_diploid_contigs(&self.store, aln).await?; + if row.source_sig == archaic_segment_sig(aln, &contigs) { Ok(Some(serde_json::from_str(&row.segments)?)) } else { Ok(None) @@ -3583,7 +3584,9 @@ impl App { "archaic segments currently require a CHM13 alignment (the Tier B assets are CHM13-only)".into(), )); } - let sig = archaic_segment_sig(aln); + // Computed from the contigs actually cached, so a later genome-wide pass invalidates a + // partial result instead of inheriting it. + let sig = archaic_segment_sig(aln, &crate::called_diploid_contigs(&self.store, aln).await?); if let Some(row) = consensus_archaic_segments::get(self.store.pool(), biosample_guid).await? { if row.source_sig == sig { return Ok(serde_json::from_str(&row.segments)?); diff --git a/crates/navigator-app/src/lib.rs b/crates/navigator-app/src/lib.rs index 7416831..e06ab16 100644 --- a/crates/navigator-app/src/lib.rs +++ b/crates/navigator-app/src/lib.rs @@ -1122,18 +1122,40 @@ fn archaic_marker_dist_path(build: ReferenceBuild) -> PathBuf { /// Tier B: positions variable in the African outgroup, for stripping shared variants. /// Cache signature for a Tier B segment result: the alignment it came from plus the caller's /// genotype version, so re-calling with a newer caller invalidates it. -pub(crate) fn archaic_segment_sig(alignment_id: i64) -> String { - // The METHOD version is part of the key, not just the genotype version. Tier B was rebuilt from - // a private-variant density model to archaic-genome matching, and without this a workspace - // carrying results from the old caller would keep serving them — stale answers from a method - // that was withdrawn for having no per-individual signal. +pub(crate) fn archaic_segment_sig(alignment_id: i64, called_contigs: &[String]) -> String { + // Three things make a result stale, and all three are in the key. + // + // The METHOD version, because Tier B was rebuilt from a private-variant density model to + // archaic-genome matching; without it a workspace carrying output from the withdrawn caller + // would keep serving it. + // + // The CONTIGS ACTUALLY CALLED, because the result covers only those. A subject called on chr21 + // alone and later called genome-wide would otherwise keep the two-chromosome answer forever — + // observed doing exactly that during genome-wide validation, reporting 1.94 Mb over 2 contigs + // while 22 sat cached and ready. + let mut contigs: Vec<&str> = called_contigs.iter().map(String::as_str).collect(); + contigs.sort_unstable(); format!( - "aln{alignment_id}:gt{}:m{}", + "aln{alignment_id}:gt{}:m{}:c[{}]", navigator_analysis::caller::GENOTYPE_VERSION, - navigator_analysis::archaic_match::METHOD_VERSION + navigator_analysis::archaic_match::METHOD_VERSION, + contigs.join(",") ) } +/// The autosomes that have cached de-novo diploid calls for `alignment_id`, as bare contig names. +pub(crate) async fn called_diploid_contigs( + store: &navigator_store::Store, + alignment_id: i64, +) -> Result, AppError> { + const PREFIX: &str = "diploid_denovo:"; + Ok(navigator_store::artifact::list_kinds(store.pool(), alignment_id) + .await? + .into_iter() + .filter_map(|k| k.strip_prefix(PREFIX).map(str::to_string)) + .collect()) +} + fn archaic_outgroup_path(build: ReferenceBuild) -> PathBuf { ancestry_asset_path("NAVIGATOR_ARCHAIC_OUTGROUP", "archaic_outgroup_af", build, "bin") } diff --git a/crates/navigator-store/src/artifact.rs b/crates/navigator-store/src/artifact.rs index 30132a7..242acfa 100644 --- a/crates/navigator-store/src/artifact.rs +++ b/crates/navigator-store/src/artifact.rs @@ -136,6 +136,20 @@ pub async fn list_for_alignment(pool: &SqlitePool, alignment_id: i64) -> Result< rows.into_iter().map(Row::into_domain).collect() } +/// The `kind` of every artifact on `alignment_id`, without their payloads. +/// +/// Exists so a caller can key a cache on *which* analyses are present without reading them. The +/// genome-wide de-novo calls run to ~1 GB of JSON across 22 contigs, so [`list_for_alignment`] is +/// the wrong tool for a question about coverage. +pub async fn list_kinds(pool: &SqlitePool, alignment_id: i64) -> Result, StoreError> { + let rows: Vec<(String,)> = + sqlx::query_as("SELECT kind FROM analysis_artifact WHERE alignment_id = ? ORDER BY kind") + .bind(alignment_id) + .fetch_all(pool) + .await?; + Ok(rows.into_iter().map(|(k,)| k).collect()) +} + /// Every artifact belonging to any of `alignment_ids`, in one query. The caller indexes the result /// by `(alignment_id, kind)` itself — this replaces a `get` per (alignment, kind), which for a /// project report meant one round-trip per cell. An empty `alignment_ids` yields no query.