From 2152cb02dffb05d1cda034b609717828f7eb82ea Mon Sep 17 00:00:00 2001 From: James Kane Date: Fri, 31 Jul 2026 06:48:39 -0500 Subject: [PATCH] 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()