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_match_probe.rs b/crates/navigator-analysis/examples/archaic_match_probe.rs new file mode 100644 index 0000000..45c1213 --- /dev/null +++ b/crates/navigator-analysis/examples/archaic_match_probe.rs @@ -0,0 +1,109 @@ +//! 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())? + }; + + // `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!( + "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/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_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/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/crates/navigator-analysis/src/archaic_match.rs b/crates/navigator-analysis/src/archaic_match.rs new file mode 100644 index 0000000..e3e636d --- /dev/null +++ b/crates/navigator-analysis/src/archaic_match.rs @@ -0,0 +1,909 @@ +//! 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. +//! +//! # 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 +//! statistic matched and the statistic was then reported as evidence. +//! +//! | | density caller | this, uncalibrated | this, calibrated | +//! |---|---|---|---| +//! | 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. +//! +//! ## 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). +//! +//! ## 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. +//! +//! ## 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; + +use serde::{Deserialize, Serialize}; + +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 { + 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, + // 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, + // 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, + } + } +} + +/// 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 } +} + +/// 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::{ + ArchaicCall, ArchaicPanelThresholds, ArchaicSite, 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)); + // 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() + }; + 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"); + } + + 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] + 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/crates/navigator-app/src/haplogroup.rs b/crates/navigator-app/src/haplogroup.rs index c2d9c3a..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)?); @@ -3612,19 +3615,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 +3644,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..e06ab16 100644 --- a/crates/navigator-app/src/lib.rs +++ b/crates/navigator-app/src/lib.rs @@ -1122,8 +1122,38 @@ 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 { - format!("aln{alignment_id}:gt{}", navigator_analysis::caller::GENOTYPE_VERSION) +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{}:c[{}]", + navigator_analysis::caller::GENOTYPE_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 { @@ -2451,36 +2481,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-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. 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( 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/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/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/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() 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() 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__": 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/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/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/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() 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/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() 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/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() 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()