Docs/archaic observable diagnosis

crates/navigator-analysis/examples/archaic_classify_dump.rs

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

crates/navigator-analysis/examples/archaic_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<dyn std::error::Error>> {
+    let a: Vec<String> = std::env::args().skip(1).collect();
+    if a.len() < 5 {
+        eprintln!(
+            "usage: archaic_match_probe <classify.bin> <callable.bin> <reference.fa> \
+             <genetic_map.bin|-> <calls.json> [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<SiteGenotype> = Vec::new();
+    for p in &a[4..] {
+        let mut v: Vec<SiteGenotype> = serde_json::from_str(&std::fs::read_to_string(p)?)?;
+        calls.append(&mut v);
+    }
+    let mut by_contig: BTreeMap<String, BTreeMap<i64, &SiteGenotype>> = 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<String, Vec<SiteObs>> = 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(())
+}

crates/navigator-analysis/examples/archaic_outgroup_density.rs

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

crates/navigator-analysis/examples/archaic_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<dyn std::error::Error>> {
+    let mut a = std::env::args().skip(1);
+    let path = a.next().expect("usage: archaic_panel_dump <archaic_markers.bin> [contig ...]");
+    let want: Vec<String> = 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(())
+}

crates/navigator-analysis/examples/archaic_private_dump.rs

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