diff --git a/crates/navigator-analysis/examples/archaic_callable_dump.rs b/crates/navigator-analysis/examples/archaic_callable_dump.rs new file mode 100644 index 0000000..a2fbbc6 --- /dev/null +++ b/crates/navigator-analysis/examples/archaic_callable_dump.rs @@ -0,0 +1,51 @@ +//! Dump the Tier B callability mask as BED, so what the segment caller can and cannot see is +//! checkable against an external callset rather than assumed. +//! +//! Windows below `min_frac` of `window_bp` callable are excluded by the caller itself, so the same +//! threshold is applied here — the output is the territory a segment could actually be called in. +//! +//! ```sh +//! cargo run --release -p navigator-analysis --example archaic_callable_dump -- \ +//! ~/.decodingus/ancestry/archaic_callable_chm13v2.0.bin 0.5 chr21 chr22 > callable.bed +//! ``` + +use navigator_analysis::archaic::ArchaicCallable; + +fn main() -> Result<(), Box> { + let mut a = std::env::args().skip(1); + let path = a.next().expect("usage: archaic_callable_dump [min_frac] [contig ...]"); + let min_frac: f64 = a.next().and_then(|s| s.parse().ok()).unwrap_or(0.0); + let want: Vec = a.collect(); + + let cal = ArchaicCallable::from_bytes(&std::fs::read(&path)?).map_err(|e| e.to_string())?; + eprintln!( + "callable track: {} contigs, window {} bp, {:.1} Mb callable total", + cal.contigs.len(), + cal.window_bp, + cal.callable_mb() + ); + + let w = cal.window_bp; + for c in &cal.contigs { + if !want.is_empty() && !want.contains(&c.contig) { + continue; + } + let (mut run_start, mut in_run) = (0i64, false); + for (i, &bp) in c.callable_bp.iter().enumerate() { + let pos = c.start + (i as i64) * w; + let ok = (bp as f64) / (w as f64) >= min_frac; + if ok && !in_run { + run_start = pos; + in_run = true; + } else if !ok && in_run { + println!("{}\t{}\t{}", c.contig, run_start, pos); + in_run = false; + } + } + if in_run { + let end = c.start + (c.callable_bp.len() as i64) * w; + println!("{}\t{}\t{}", c.contig, run_start, end); + } + } + Ok(()) +} diff --git a/crates/navigator-analysis/examples/cram_query_probe.rs b/crates/navigator-analysis/examples/cram_query_probe.rs new file mode 100644 index 0000000..82ff181 --- /dev/null +++ b/crates/navigator-analysis/examples/cram_query_probe.rs @@ -0,0 +1,92 @@ +//! Time the phases of an indexed region read, so a slow BAM/CRAM path can be attributed to the +//! part actually responsible (open, header, first query, warm query, bulk region iteration) rather +//! than to whichever call the wall-clock happened to land in. +//! +//! Written to diagnose a CRAM that took ~500x longer than an equivalent BAM for one region query. +//! +//! ```sh +//! cargo run --release -p navigator-analysis --example cram_query_probe -- \ +//! [span_bp=1000000] +//! ``` + +use std::path::Path; +use std::time::Instant; + +use navigator_analysis::reader::open_indexed; +use noodles::core::Region; + +fn main() { + let a: Vec = std::env::args().collect(); + if a.len() < 5 { + eprintln!("usage: cram_query_probe [span_bp=1000000]"); + std::process::exit(2); + } + let (path, refp, contig) = (Path::new(&a[1]), Path::new(&a[2]), a[3].clone()); + let pos: usize = a[4].parse().expect("pos"); + let span: usize = a.get(5).and_then(|s| s.parse().ok()).unwrap_or(1_000_000); + + let t0 = Instant::now(); + let (header, mut reader) = open_indexed(path, Some(refp)).expect("open"); + println!("open + header : {:>8.2?}", t0.elapsed()); + + // A single-base region: cost here is per-query overhead, not per-record work. + let one = |p: usize| -> Region { format!("{contig}:{p}-{p}").parse().expect("region") }; + + let t = Instant::now(); + let n: usize = reader.query(&header, &one(pos)).expect("q1").count(); + println!("first 1bp query ({n:>4} rec) : {:>8.2?} <- includes any lazy setup", t.elapsed()); + + let t = Instant::now(); + let n: usize = reader.query(&header, &one(pos)).expect("q2").count(); + println!("same query again ({n:>4} rec): {:>8.2?} <- warm: is the cost per-query or one-off?", t.elapsed()); + + let t = Instant::now(); + let n: usize = reader.query(&header, &one(pos + 5_000_000)).expect("q3").count(); + println!("distant 1bp query ({n:>4} rec): {:>8.2?} <- new container: does it re-decode?", t.elapsed()); + + let region: Region = format!("{contig}:{pos}-{}", pos + span).parse().expect("region"); + let t = Instant::now(); + let n: usize = reader.query(&header, ®ion).expect("bulk").count(); + let el = t.elapsed(); + println!( + "{:.1} Mb region ({n} rec) : {:>8.2?} = {:.2?}/Mb", + span as f64 / 1e6, + el, + el / (span as u32 / 1_000_000).max(1) + ); + + // VERIFY=1 checks the container-skipping query against noodles' own (whole-contig) Query on + // REAL data. The checked-in fixture is a single container, so only a large multi-container CRAM + // can catch a container wrongly skipped -- which would present as a faster caller, not a broken + // one. Slow by construction: the oracle is the implementation we replaced. + if std::env::var("VERIFY").is_ok_and(|v| v == "1") { + use noodles::cram; + + let key = |r: &noodles::sam::alignment::RecordBuf| { + ( + r.name().map(|n| n.to_vec()), + r.alignment_start().map(usize::from), + r.flags().bits(), + r.sequence().as_ref().to_vec(), + ) + }; + let mine: Vec<_> = reader.query(&header, ®ion).expect("mine").map(|r| key(&r.expect("rec"))).collect(); + + let repo = navigator_analysis::reader::build_repository(refp).expect("repo"); + let mut oracle = cram::io::indexed_reader::Builder::default() + .set_reference_sequence_repository(repo) + .build_from_path(path) + .expect("noodles open"); + let oh = oracle.read_header().expect("noodles header"); + let t = Instant::now(); + let theirs: Vec<_> = oracle + .query(&oh, ®ion) + .expect("noodles query") + .map(|r| key(&r.expect("rec"))) + .collect(); + println!("\nVERIFY: noodles' own Query took {:?} for the same region", t.elapsed()); + println!(" ours {} records, noodles {} records", mine.len(), theirs.len()); + assert_eq!(mine, theirs, "container skipping changed the records returned"); + println!(" IDENTICAL — container skipping is lossless"); + } +} diff --git a/crates/navigator-analysis/src/reader.rs b/crates/navigator-analysis/src/reader.rs index 2919567..f8a8568 100644 --- a/crates/navigator-analysis/src/reader.rs +++ b/crates/navigator-analysis/src/reader.rs @@ -12,7 +12,8 @@ use std::fs::File; use std::num::NonZeroUsize; use std::path::{Path, PathBuf}; -use noodles::core::Region; +use noodles::core::region::Interval; +use noodles::core::{Position, Region}; use noodles::sam::alignment::RecordBuf; use noodles::{bam, bgzf, cram, fasta, sam}; @@ -216,6 +217,37 @@ pub enum IdxReader { }, } +/// File offsets of the `.crai` containers that can hold records overlapping `interval` on `ref_id`. +/// +/// **This is the whole reason CRAM region queries are usable.** A CRAM container is the unit of +/// decode — you cannot decode part of one — so restricting *which containers get decoded* is the +/// only place a region query can save work. noodles' own `Query` (and, before this, our `for_each`) +/// selected containers by reference sequence alone and then discarded non-overlapping records +/// *after* decoding them, which made every query cost a whole chromosome no matter how small the +/// region: measured at 20.9 s for a 1 bp query on chr21 and 116 s on chr1, against 4–6 ms for the +/// same query on a BAM. chr21 of a 30x WGS holds 1,140 containers and a point query needs exactly +/// one of them. +/// +/// A container whose `alignment_start` is absent is **kept**: that is a container this index cannot +/// place, and dropping it would silently lose records. Skipping is only ever done on positive +/// evidence that the container lies outside the interval. +fn cram_container_offsets(index: &cram::crai::Index, ref_id: usize, interval: Interval) -> Vec { + index + .iter() + .filter(|r| r.reference_sequence_id() == Some(ref_id)) + .filter(|r| match r.alignment_start() { + Some(start) => { + // Span 0 would make an empty range, which intersects nothing; treat it as one base. + let span = r.alignment_span().max(1); + let end = Position::new(usize::from(start).saturating_add(span - 1)).unwrap_or(start); + interval.intersects((start..=end).into()) + } + None => true, + }) + .map(|r| r.offset()) + .collect() +} + /// Open `path` for indexed region queries (autoloads the `.bai`/`.crai`). `reference` is /// required for CRAM. pub fn open_indexed(path: &Path, reference: Option<&Path>) -> Result<(sam::Header, IdxReader), AnalysisError> { @@ -268,10 +300,87 @@ impl IdxReader { RecordBuf::try_from_alignment_record(header, &rec).map_err(|e| AnalysisError::io(&path, e)) }))) } - IdxReader::Cram { inner, path, .. } => { + // Hand-rolled rather than `inner.query(...)`: noodles' `Query` decodes every container + // of the contig and filters records afterwards, so a 1 bp query costs a whole + // chromosome (see [`cram_container_offsets`]). This decodes only the containers that + // can overlap, lazily — one container at a time, so a caller that stops early (a + // `.take(n)` probe, a cancelled walk) does not pay for the rest. + IdxReader::Cram { inner, repo, path } => { + use std::io::{Seek, SeekFrom}; + + use noodles::sam::alignment::Record as _; // alignment_start/_end on cram::Record + let path = path.clone(); - let q = inner.query(header, region).map_err(|e| AnalysisError::io(&path, e))?; - Ok(Box::new(q.map(move |r| r.map_err(|e| AnalysisError::io(&path, e))))) + let repo = repo.clone(); + let ref_id = header.reference_sequences().get_index_of(region.name()).ok_or_else(|| { + AnalysisError::Message(format!( + "contig {} not in {} header", + String::from_utf8_lossy(region.name()), + path.display() + )) + })?; + let interval = region.interval(); + let mut offsets = cram_container_offsets(inner.index(), ref_id, interval).into_iter(); + + let mut pending: std::vec::IntoIter = Vec::new().into_iter(); + let mut container = cram::io::reader::Container::default(); + Ok(Box::new(std::iter::from_fn(move || { + loop { + if let Some(rec) = pending.next() { + return Some(Ok(rec)); + } + // Next container that can overlap; `None` ends the iterator. + let offset = offsets.next()?; + let io_err = |e| AnalysisError::io(&path, e); + if let Err(e) = inner.get_mut().seek(SeekFrom::Start(offset)).map_err(io_err) { + return Some(Err(e)); + } + match inner.read_container(&mut container).map_err(io_err) { + Ok(0) => continue, + Ok(_) => {} + Err(e) => return Some(Err(e)), + } + let compression_header = match container.compression_header().map_err(io_err) { + Ok(h) => h, + Err(e) => return Some(Err(e)), + }; + let mut buf = Vec::new(); + for slice in container.slices() { + let slice = match slice.map_err(io_err) { + Ok(s) => s, + Err(e) => return Some(Err(e)), + }; + let (core, external) = match slice.decode_blocks().map_err(io_err) { + Ok(b) => b, + Err(e) => return Some(Err(e)), + }; + let records = match slice + .records(repo.clone(), header, &compression_header, &core, &external) + .map_err(io_err) + { + Ok(r) => r, + Err(e) => return Some(Err(e)), + }; + for rec in &records { + // Same per-record overlap test noodles applies post-decode — the + // container filter is a coarse prefilter, not a replacement for it. + if let (Some(Ok(start)), Some(Ok(end))) = (rec.alignment_start(), rec.alignment_end()) + { + if !interval.intersects((start..=end).into()) { + continue; + } + } else { + continue; + } + match RecordBuf::try_from_alignment_record(header, rec).map_err(io_err) { + Ok(r) => buf.push(r), + Err(e) => return Some(Err(e)), + } + } + } + pending = buf.into_iter(); + } + }))) } } } @@ -372,14 +481,12 @@ impl IdxReader { })?; let interval = region.interval(); - // Collect the file offsets of this contig's containers before borrowing `inner` - // mutably to seek/read (the `.crai` index borrow can't overlap the read borrow). - let offsets: Vec = inner - .index() - .iter() - .filter(|r| r.reference_sequence_id() == Some(ref_id)) - .map(|r| r.offset()) - .collect(); + // Collect the file offsets of the containers that can overlap the query before + // borrowing `inner` mutably to seek/read (the `.crai` index borrow can't overlap the + // read borrow). Selecting on the interval — not just the contig — is what keeps this + // proportional to the region instead of the chromosome; see + // [`cram_container_offsets`]. + let offsets = cram_container_offsets(inner.index(), ref_id, interval); let mut container = cram::io::reader::Container::default(); for offset in offsets { @@ -578,4 +685,79 @@ mod tests { assert!(!sink.0.is_empty(), "fixture should have chrM records"); assert_eq!(sink.0, via_query, "cram::Record path must match RecordBuf path"); } + + /// [`cram_container_offsets`] decides which containers are decoded at all, so its boundary + /// behaviour *is* the correctness of every CRAM region query: a container wrongly skipped is + /// reads silently missing from a variant call, which no downstream test would attribute to the + /// reader. Checked at the edges, where an off-by-one actually lives. + #[test] + fn container_offsets_select_only_overlapping_containers() { + let p = |n: usize| Position::new(n).unwrap(); + // ref 0 containers spanning [1000,1099], [2000,2099], [3000,3099]; one on ref 1; and one + // the index cannot place. + let idx: cram::crai::Index = vec![ + cram::crai::Record::new(Some(0), Some(p(1000)), 100, 10, 0, 0), + cram::crai::Record::new(Some(0), Some(p(2000)), 100, 20, 0, 0), + cram::crai::Record::new(Some(0), Some(p(3000)), 100, 30, 0, 0), + cram::crai::Record::new(Some(1), Some(p(2000)), 100, 40, 0, 0), + cram::crai::Record::new(Some(0), None, 0, 50, 0, 0), + ]; + let sel = |a: usize, b: usize| cram_container_offsets(&idx, 0, (p(a)..=p(b)).into()); + + // A point inside one container decodes that container — not the contig. This single + // assertion is the difference between 8 ms and 21 s on a real chr21. + assert_eq!(sel(2050, 2050), vec![20, 50]); + // Boundaries: touching the first/last base of a container counts as overlap. + assert_eq!(sel(2099, 2099), vec![20, 50], "last base of a container overlaps"); + assert_eq!(sel(2000, 2000), vec![20, 50], "first base of a container overlaps"); + assert_eq!(sel(2100, 2100), vec![50], "one past the end does not"); + assert_eq!(sel(1999, 1999), vec![50], "one before the start does not"); + // A span crossing several containers takes exactly those it crosses. + assert_eq!(sel(1050, 2050), vec![10, 20, 50]); + // The other reference is never selected, even at identical coordinates. + assert_eq!(cram_container_offsets(&idx, 1, (p(2050)..=p(2050)).into()), vec![40]); + // An unbounded interval keeps every container on the reference. + assert_eq!( + cram_container_offsets(&idx, 0, Region::new(b"x".to_vec(), ..).interval()), + vec![10, 20, 30, 50] + ); + } + + /// Our hand-rolled `query` must return exactly what noodles' own `Query` returns. We replaced + /// it for speed, and a reimplementation that quietly drops records would look like a faster + /// caller rather than a broken one — the failure mode worth a test. + #[test] + fn cram_query_matches_noodles_query() { + let dir = PathBuf::from(env!("CARGO_MANIFEST_DIR")).join("tests/fixtures"); + let cram = dir.join("coverage.cram"); + let reference = dir.join("ref.fa"); + + for region in [ + Region::new(b"chrM".to_vec(), ..), + Region::new(b"chrM".to_vec(), Position::new(1).unwrap()..=Position::new(200).unwrap()), + Region::new(b"chrM".to_vec(), Position::new(50).unwrap()..=Position::new(60).unwrap()), + ] { + let (header, mut ours) = open_indexed(&cram, Some(&reference)).expect("open"); + let mine: Vec = ours + .query(&header, ®ion) + .expect("query") + .map(|r| capture(&r.expect("rec"))) + .collect(); + + // noodles' own indexed query, unmodified, as the oracle. + let repo = build_repository(&reference).expect("repo"); + let mut theirs = cram::io::indexed_reader::Builder::default() + .set_reference_sequence_repository(repo) + .build_from_path(&cram) + .expect("noodles open"); + let nheader = theirs.read_header().expect("noodles header"); + let reference_impl: Vec = theirs + .query(&nheader, ®ion) + .expect("noodles query") + .map(|r| capture(&r.expect("rec"))) + .collect(); + + assert_eq!(mine, reference_impl, "region {region:?}: must match noodles' Query exactly"); + } + } } diff --git a/crates/navigator-app/src/haplogroup.rs b/crates/navigator-app/src/haplogroup.rs index 9722e43..c2d9c3a 100644 --- a/crates/navigator-app/src/haplogroup.rs +++ b/crates/navigator-app/src/haplogroup.rs @@ -3515,10 +3515,17 @@ impl App { /// The cached Tier B archaic segment result for a subject, if current for the alignment and /// caller version it was produced from. + /// + /// Gated by [`crate::ARCHAIC_SEGMENTS_ENABLED`] on the **read** path as well as the compute + /// path: rows persisted before the gate went in are still in the workspace, and a read-only + /// gate is the difference between withholding a result and merely declining to recompute it. pub async fn cached_archaic_segments( &self, biosample_guid: SampleGuid, ) -> Result, AppError> { + if !crate::ARCHAIC_SEGMENTS_ENABLED { + return Ok(None); + } let Some(row) = consensus_archaic_segments::get(self.store.pool(), biosample_guid).await? else { return Ok(None); }; @@ -3545,6 +3552,18 @@ impl App { &self, biosample_guid: SampleGuid, ) -> Result { + // Withheld: the caller reproduces the cohort mean and nothing about the individual. An + // error rather than an empty result, because every caller of this asked for a computation — + // silently returning zero segments would read as "you have no archaic ancestry", which is + // a far worse claim than "we are not reporting this". + if !crate::ARCHAIC_SEGMENTS_ENABLED { + return Err(AppError::Import( + "archaic segment calling is disabled: validated against hmmix's per-individual \ + calls (n=20), it reproduced the cohort mean but showed no per-person signal — \ + locations below chance and extent r = -0.02. See ARCHAIC_SEGMENTS_ENABLED." + .into(), + )); + } // Prefer an alignment that already HAS genome-wide diploid calls over the // highest-coverage one. Those calls are an hours-long per-alignment pass, so a subject with // several CHM13 alignments (this one has four) would otherwise be told to re-run work they diff --git a/crates/navigator-app/src/lib.rs b/crates/navigator-app/src/lib.rs index 67a1396..6364c84 100644 --- a/crates/navigator-app/src/lib.rs +++ b/crates/navigator-app/src/lib.rs @@ -2450,6 +2450,38 @@ pub struct AncientFitRow { /// See `documents/design/ancient-ancestry-rebuild.md` (start at §7.14). 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. +/// +/// 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)". +pub const ARCHAIC_SEGMENTS_ENABLED: bool = false; + /// 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. pub use navigator_analysis::ancestry::ANCIENT_ADMIXTURE; diff --git a/crates/navigator-domain/locales/en.txt b/crates/navigator-domain/locales/en.txt index 0ef5868..e447eef 100644 --- a/crates/navigator-domain/locales/en.txt +++ b/crates/navigator-domain/locales/en.txt @@ -803,3 +803,7 @@ simple.relatives.cmUnit=cM shared simple.relatives.segUnit=segments simple.relatives.confirmed=Measured by a DNA exchange — simple.relatives.estimate=Estimated from shared signals — connect to measure how much DNA you actually share. + +# 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. diff --git a/crates/navigator-domain/locales/es.txt b/crates/navigator-domain/locales/es.txt index b4e70ab..696fd82 100644 --- a/crates/navigator-domain/locales/es.txt +++ b/crates/navigator-domain/locales/es.txt @@ -788,3 +788,7 @@ simple.relatives.cmUnit=cM compartidos simple.relatives.segUnit=segmentos simple.relatives.confirmed=Medido por un intercambio de ADN: simple.relatives.estimate=Estimado a partir de señales compartidas: conecta para medir cuánto ADN compartís realmente. + +# 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. diff --git a/crates/navigator-ui/src/ui/central.rs b/crates/navigator-ui/src/ui/central.rs index db2cbab..5aa4575 100644 --- a/crates/navigator-ui/src/ui/central.rs +++ b/crates/navigator-ui/src/ui/central.rs @@ -617,8 +617,21 @@ impl NavigatorApp { } // Tier B: archaic SEGMENTS (WGS only — needs genome-wide de-novo calls). + // Withheld pending a working method: the card states that rather than + // disappearing, because a section that silently vanishes between releases reads + // as a bug, where a stated withholding is a finding about the data. ui.add_space(10.0); card(ui, self.tr("card.archaicSegments"), |ui| { + if !navigator_app::ARCHAIC_SEGMENTS_ENABLED { + ui.label(self.tr("archaicSegments.withheld")); + ui.add_space(4.0); + ui.label( + egui::RichText::new(self.tr("archaicSegments.withheldWhy")) + .weak() + .small(), + ); + return; + } ui.horizontal(|ui| { let have = self.archaic_segments.is_some(); let label = if have { diff --git a/documents/design/ArchaicAncestry_Design.md b/documents/design/ArchaicAncestry_Design.md index d8e6906..3174393 100644 --- a/documents/design/ArchaicAncestry_Design.md +++ b/documents/design/ArchaicAncestry_Design.md @@ -1,10 +1,11 @@ # Archaic Ancestry Report (Neanderthal / Denisovan) — Design -**Status:** **SHIPPED** in `v0.1.0-alpha.14` (2026-07-30). Tier A (#34) and Tier B (#35) are both in -`main` and in users' hands; §8 Phases 1 and 2 are complete, Phase 3 (M4) is not started and remains -optional. Drafted 2026-07-23; plan added 2026-07-26 on branch `feat/archaic-ancestry`; all three §9 -questions resolved. **Two things shipped differently from the plan below — see *Deviations from the -plan* at the end of §10 before trusting §7's expected percentage or M3's feature-gate rule.** +**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.** **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. @@ -956,14 +957,88 @@ genuine open gap, not a deviation that resolved itself.** ### Open follow-ups - **Write the Denisovan-floor regression test** (item 3 above). -- **Validate extent across several individuals**, not one. Needs WGS for people whose hmmix result is - known. Until then the 1.01× agreement is a calibration check, not a validation. +- ~~**Validate extent across several individuals**~~ — **DONE, and Tier B failed it.** See below. - **Lineage attribution** needs the Skov-2020 approach — match a segment's own haplotype against each archaic genome relative to a background expectation, rather than against pre-classified site - categories. Threshold tuning will not fix it; see the attribution section above. -- **Asset staging policy was never decided.** M1 said to check Asset 1's size against - `ON_DEMAND_PREFIXES` in `packaging/stage-assets.sh` "before deciding"; that decision was never - made. As a result `PATTERNS` lists no `archaic_*`, so dev-mode staging omits them, while - release-mode fetches everything the manifest names except `ancestry_haps_` and bundles all five - (105.3 MB). The two modes disagree — exactly what the comment above `PATTERNS` warns against — and - alpha.14's installers are ~60 MB larger than alpha.13's as a result. Tracked separately. + categories. Threshold tuning will not fix it; see the attribution section above. The same paper is + now also the path back for the segments themselves. +- ~~**Asset staging policy was never decided**~~ — decided 2026-07-30 (PR #38): Tier A bundled, + Tier B on demand. See `packaging-and-release.md`. + +--- + +## Tier B validation (2026-07-30) — no per-individual signal; GATED OFF + +The follow-up above ("validate extent across several individuals") was run, and Tier B did not +survive it. `ARCHAIC_SEGMENTS_ENABLED = false`. + +**What made it possible.** The workspace holds CHM13 CRAMs for **2,307 of the 2,309** individuals in +hmmix's own published callset, including 632 of the 633 Europeans — so the comparison is against +*the same people*, not against a cohort distribution. (Getting there required fixing a CRAM +region-query defect that made every query cost a whole chromosome; chr21+22 per sample went from +~3.5 hours to ~90 seconds. See PR #39.) + +### Locations disagree — below chance + +HG00096, chr21+22, our segments against hmmix's for the same individual, lifted hg38→CHM13: + +| | | +|---|---| +| our extent | 2.294 Mb in 112 segments | +| hmmix extent | 2.333 Mb in 48 tracts | +| **base overlap** | **0.050 Mb** | +| sensitivity / precision | **2.1 % / 1.5 %** | +| null: our own segment lengths placed at random in the same span | **5.0 %** (p95 9.4 %) | + +**We score below chance.** Every alternative explanation was tested and ruled out: + +- **Not a coordinate error.** Overlap-vs-shift is flat (1–4 %) across ±2 Mb with no peak; zero-shift + sits mid-range. The lift is faithful — lifted fragment lengths sum to the hg38 input exactly. +- **Not callability.** 70.7 % of hmmix's tracts lie inside our callable territory, so they were + reachable. Restricted to reachable truth, sensitivity is 3.0 %. +- **Not haplotype handling.** hmmix is per-haplotype and is unioned, not summed — which reproduces + their published EUR mean of 2.09 Mb exactly. +- **Not a harness artefact.** Overlap was re-derived brute-force. (One real harness bug was found and + fixed: CrossMap splits tracts at median 2 bp gaps, inflating 48 tracts to 423.) + +### Amounts do not track the individual either + +n = 20 Europeans, randomly drawn, spanning the natural range: + +| statistic | value | +|---|---| +| Pearson r (our extent vs theirs) | **−0.018** (permutation p = 0.94) | +| Spearman ρ | **−0.020** (p = 0.94) | +| segment-count r | −0.222 | +| mean ratio ours/theirs | **0.923** | +| SD: truth 0.496 Mb, ours 0.312 Mb | ours = **0.63×** the truth's spread | + +Truth ranged 1.19–2.97 Mb; our output barely moved. The two individuals with the *least* archaic +ancestry drew our two *highest* calls (1.71×, 1.90×), and the two with the most drew among our +lowest (0.54×, 0.63×). + +### What this means about the earlier "validation" + +The M3 result — 1.01× the hmmix EUR mean — was **three fitted parameters hitting the statistic they +were fitted to**, on one individual. It was recorded honestly as a calibration check rather than a +validation, and that caution was right: the fit does not transfer. A caller emitting a calibrated +constant passes a cohort-mean test and fails every test that asks about a person. + +The general lesson, which is the same one the ancient-ancestry work produced: **agreement with an +aggregate is not evidence of a per-individual measurement.** Any future re-enable must clear the +per-individual bar — positive r against held-out individuals, and location overlap well above the +random-placement null — not a mean. + +### What is NOT affected + +**Tier A is a different method and is unaffected.** The marker count is direct dosage over a fixed +panel with no HMM, no fitted thresholds, and was checked differently (per-site archaic rate on the +intersection with real 23andMe v5 chip content). It remains enabled and is what the Simple-mode +"Neanderthal ancestry" card and the Advanced count + percentile report. + +### 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. diff --git a/scripts/archaic-validation/README.md b/scripts/archaic-validation/README.md new file mode 100644 index 0000000..6e2d51e --- /dev/null +++ b/scripts/archaic-validation/README.md @@ -0,0 +1,77 @@ +# Tier B archaic-segment validation + +The harness that produced the numbers in +[`documents/design/ArchaicAncestry_Design.md`](../../documents/design/ArchaicAncestry_Design.md) +§ *Tier B validation*, and the gate for ever setting `ARCHAIC_SEGMENTS_ENABLED` back to `true`. + +## Why this exists + +Tier B was shipped on one number: its total archaic extent landed at 1.01× the hmmix European mean. +That number was produced by three parameters fitted until it did. **A caller that emits a calibrated +constant passes a cohort-mean test**, so the mean says nothing about whether the caller measures the +person in front of it. + +These scripts ask the two questions that do: + +1. **Locations** — do we call archaic sequence *where* an independent callset does, better than + random placement would? (`compare_locations.py`) +2. **Amounts** — across individuals, does our extent rise and fall with theirs? (`correlate_extent.py`) + +Tier B scored *below* the random null on (1) and r = −0.018 on (2), which is why it is gated off. + +## What you need + +- **hmmix's published 1000G callset** (Zenodo, CC BY 4.0) — the same record checkpoint A used. + `hmmix_segments_chr21_22.tsv` covers 2,310 individuals across four regions on chr21+22; + `hmmix_eur_all.tsv` is genome-wide but Europeans only. Both are hg38. +- **CHM13 alignments in the workspace for those same individuals.** The comparison is only worth + anything against *the same people* — a distribution comparison is what got us here. +- `CrossMap` and `~/.decodingus/liftover/hg38ToHs1.over.chain` for the lift. + +## Procedure + +**1. Call the contigs the truth covers.** chr21+22 is the right axis: the truth covers every +individual there, and `pct_callable` restricts its denominator to the contigs actually called, so a +partial run is directly comparable rather than silently wrong. + +```sh +ls runs || mkdir runs +for s in HG00096 HG00133 ...; do ./run_one.sh "$s" "$PWD" /path/to/navigator; done # ~90 s each +``` + +Run 2 at a time (`xargs -P 2`): one sample saturates ~6 cores of a 16-core machine, so two fill it. + +**2. Build the truth, lifted to CHM13.** Emit one BED line per hmmix segment **with a unique id**, +CrossMap it, then reassemble. + +Two traps, both of which produced wrong answers before they were caught: + +- **Reassemble from the lifted fragments, not from min..max of them.** CrossMap splits a segment + into many pieces, and a couple land across a rearranged region — taking the span inflated + HG00096's 2.3 Mb to 26.6 Mb. The fragment-length sum should equal the hg38 input exactly; check it. +- **Merge fragments with a ~1 kb tolerance.** The lift leaves median 2 bp gaps, which a strict union + will not close, turning 48 real tracts into 423 shards and making per-tract recovery meaningless. + +**3. Union across haplotypes — never sum.** hmmix reports per haplotype; our caller is unphased. +Summing doubles their figure and would make a caller look correctly calibrated when it is not. +Unioning reproduces their published EUR mean of 2.09 Mb on chr21+22, which is the check that the +truth-side arithmetic is right. + +**4. Compare.** + +```sh +python3 compare_locations.py runs/HG00096.json truth_HG00096.chm13.bed chr21 chr22 +python3 correlate_extent.py # reads runs/*.json, prints r, rho, permutation p, spread +``` + +## Reading the output + +- **Sensitivity must beat the null**, which `compare_locations.py` does not compute for you — draw + segments of your own lengths at random within the truth's span and measure the same overlap. For + Tier B that null was 5.0 % (p95 9.4 %) against an observed 2.1 %. +- **`r` near zero with a wide truth spread means a calibrated constant**, not a measurement. Check + the spread ratio too: an estimator of a varying quantity should vary about as much as the quantity + does. Tier B's was 0.63×. +- **Before concluding the caller is wrong, rule out the harness.** Coordinate offset (cross-correlate + overlap against ±shift; a real bug peaks off zero), callability (what fraction of the truth is even + reachable — use `cargo run --example archaic_callable_dump`), and the two lift traps above. diff --git a/scripts/archaic-validation/compare_locations.py b/scripts/archaic-validation/compare_locations.py new file mode 100644 index 0000000..98e08c1 --- /dev/null +++ b/scripts/archaic-validation/compare_locations.py @@ -0,0 +1,117 @@ +#!/usr/bin/env python3 +"""Compare Navigator's archaic segment calls against the lifted hmmix truth for one sample. + +Usage: compare_archaic.py [contig ...] + +`ours.json` is `navigator archaic-segments --json` output; `truth.bed` is the hmmix callset for the +same individual, lifted hg38 -> CHM13 and unioned across haplotypes (hmmix reports per haplotype, +our caller is unphased, so union -- summing would double their figure). + +Reports three things, in increasing order of how hard they are to fake: + 1. total extent -- a caller that only matches the cohort mean passes this + 2. base-level overlap -- did we find the SAME sequence, not just the same amount + 3. per-segment recovery -- how many of their tracts we hit at all +""" +import json +import sys +from collections import defaultdict + + +def load_bed(path, keep=None): + iv = defaultdict(list) + for line in open(path): + p = line.rstrip("\n").split("\t") + if len(p) < 3 or p[0].startswith("#"): + continue + 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()} + + +def union(iv): + iv = sorted(iv) + out, cs, ce = [], *iv[0] + for s, e in iv[1:]: + if s > ce: + out.append((cs, ce)) + cs, ce = s, e + else: + ce = max(ce, e) + out.append((cs, ce)) + return out + + +def total(d): + return sum(e - s for v in d.values() for s, e in v) + + +def intersect(a, b): + """Base-level intersection of two contig->intervals dicts.""" + out = 0 + for c, av in a.items(): + bv = b.get(c) + if not bv: + continue + i = j = 0 + while i < len(av) and j < len(bv): + lo, hi = max(av[i][0], bv[j][0]), min(av[i][1], bv[j][1]) + if hi > lo: + out += hi - lo + if av[i][1] < bv[j][1]: + i += 1 + else: + j += 1 + return out + + +def main(): + ours_path, truth_path = sys.argv[1], sys.argv[2] + contigs = set(sys.argv[3:]) or None + + doc = json.load(open(ours_path)) + segs = doc.get("segments", doc if isinstance(doc, list) else []) + ours_iv = defaultdict(list) + for s in segs: + c = s.get("contig") or s.get("chrom") + if contigs and c not in contigs: + 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) + + o_mb, t_mb = total(ours) / 1e6, total(truth) / 1e6 + inter = intersect(ours, truth) / 1e6 + union_mb = o_mb + t_mb - inter + + print(f"contigs compared : {sorted(set(ours) | set(truth))}") + print() + print("1. EXTENT") + print(f" ours : {o_mb:7.3f} Mb in {sum(len(v) for v in ours.values()):5d} segments") + print(f" hmmix (truth) : {t_mb:7.3f} Mb in {sum(len(v) for v in truth.values()):5d} merged tracts") + print(f" ratio ours/theirs : {o_mb / t_mb:7.3f}" if t_mb else " ratio: n/a") + print() + print("2. BASE-LEVEL AGREEMENT (the test extent alone cannot pass)") + print(f" overlap : {inter:7.3f} Mb") + print(f" sensitivity : {inter / t_mb * 100:6.1f}% of their archaic bases we also call") + print(f" precision : {inter / o_mb * 100:6.1f}% of our archaic bases they also call" if o_mb else "") + print(f" Jaccard : {inter / union_mb:7.3f}" if union_mb else "") + print() + print("3. PER-TRACT RECOVERY") + hit = miss = 0 + for c, tv in truth.items(): + ov = ours.get(c, []) + for s, e in tv: + if any(min(e, oe) > max(s, os_) for os_, oe in ov): + hit += 1 + else: + 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 __name__ == "__main__": + main() diff --git a/scripts/archaic-validation/correlate_extent.py b/scripts/archaic-validation/correlate_extent.py new file mode 100644 index 0000000..2ba6e2c --- /dev/null +++ b/scripts/archaic-validation/correlate_extent.py @@ -0,0 +1,139 @@ +#!/usr/bin/env python3 +"""Does Tier B's archaic extent track the truth PER PERSON, or only on average? + +The shipped validation showed one genome landing at 1.01x the cohort MEAN. A caller whose output +is pure noise, scaled to the right average, passes that test. This asks the question that separates +the two: across individuals, does our extent rise and fall with hmmix's? + +Reports Pearson r (linear agreement) and Spearman rho (rank agreement, robust to a scale error), +each with a permutation p-value, plus the same for segment COUNT. A wide spread in the truth is what +gives the test power, so the observed spread is printed too. +""" +import collections +import glob +import json +import math +import os +import random + +SC = os.environ.get("ARCHAIC_RUNS", os.path.dirname(os.path.abspath(__file__))) + + +def union_bp(iv): + iv = sorted(iv) + tot, cs, ce = 0, *iv[0] + for s, e in iv[1:]: + if s > ce: + tot += ce - cs + cs, ce = s, e + else: + ce = max(ce, e) + return tot + ce - cs + + +def truth_extent(): + """hmmix extent per individual on chr21+22, haplotypes UNIONED (they report per haplotype).""" + src = os.path.expanduser('~/.decodingus/ancestry-build/tmp/hmmix_segments_chr21_22.tsv') + segs = collections.defaultdict(list) + nseg = collections.Counter() + with open(src) as f: + next(f) + for line in f: + p = line.rstrip('\n').split('\t') + segs[(p[0], p[4])].append((int(p[5]), int(p[6]))) + per, cnt = collections.defaultdict(int), collections.Counter() + for (n, _c), iv in segs.items(): + per[n] += union_bp(iv) + cnt[n] += len(iv) + return per, cnt + + +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) + i = 0 + while i < len(order): + j = i + while j + 1 < len(order) and v[order[j + 1]] == v[order[i]]: + j += 1 + avg = (i + j) / 2 + 1 + for k in range(i, j + 1): + r[order[k]] = avg + i = j + 1 + return r + + +def perm_p(x, y, stat, draws=20000): + obs = abs(stat(x, y)) + rnd = random.Random(0) + yy = list(y) + hits = 0 + for _ in range(draws): + rnd.shuffle(yy) + if abs(stat(x, yy)) >= obs: + hits += 1 + return (hits + 1) / (draws + 1) + + +def main(): + per, cnt = truth_extent() + rows = [] + for f in sorted(glob.glob(f'{SC}/runs/*.json')): + name = os.path.basename(f)[:-5] + try: + d = json.load(open(f)) + except Exception: + continue + s = d.get('summary') or {} + if name not in per or not s: + continue + rows.append((name, per[name] / 1e6, s['total_mb'], cnt[name], s['n_segments'])) + + if len(rows) < 5: + print(f'only {len(rows)} samples — not enough') + return + rows.sort(key=lambda r: r[1]) + print(f'n = {len(rows)} European individuals, chr21+22\n') + print(f"{'sample':10s} {'hmmix Mb':>9s} {'ours Mb':>8s} {'ratio':>6s} {'hmmix segs':>11s} {'our segs':>9s}") + for n, t, o, tc, oc in rows: + print(f'{n:10s} {t:9.3f} {o:8.3f} {o / t:6.2f} {tc:11d} {oc:9d}') + + t = [r[1] for r in rows] + o = [r[2] for r in rows] + print(f'\ntruth spread: {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} <- the statistic the shipped validation used') + + r = pearson(t, o) + rho = pearson(rank(t), rank(o)) + print(f'\nPER-INDIVIDUAL AGREEMENT (the statistic that was never measured)') + print(f' Pearson r = {r:+.3f} permutation p = {perm_p(t, o, pearson):.4f}') + print(f' Spearman rho = {rho:+.3f} permutation p = {perm_p(rank(t), rank(o), pearson):.4f}') + tc = [float(r_[3]) for r_ in rows] + oc = [float(r_[4]) for r_ in rows] + print(f' segment count: Pearson r = {pearson(tc, oc):+.3f}') + + # Spread carries the same message without needing a correlation: a measurement of a varying + # quantity should vary about as much as the quantity does. A calibrated constant does not. + def sd(v): + m = sum(v) / len(v) + return math.sqrt(sum((x - m) ** 2 for x in v) / (len(v) - 1)) + + st, so = sd(t), sd(o) + print(f'\nSPREAD truth SD = {st:.3f} Mb ours SD = {so:.3f} Mb (ours/truth = {so / st:.2f})') + print(f' truth CV = {st / (sum(t) / len(t)):.3f} ours CV = {so / (sum(o) / len(o)):.3f}') + print('\n r ~ 0 with this spread => the total is a calibrated constant, not a measurement of') + print(' the individual. r > 0 => the headline number carries real per-person signal even') + print(' though the LOCATIONS do not.') + + +if __name__ == '__main__': + main() diff --git a/scripts/archaic-validation/run_one.sh b/scripts/archaic-validation/run_one.sh new file mode 100755 index 0000000..1fbda0c --- /dev/null +++ b/scripts/archaic-validation/run_one.sh @@ -0,0 +1,10 @@ +#!/usr/bin/env bash +set -uo pipefail +S="$1"; SC="$2"; BIN="$3" +[ -s "$SC/runs/$S.json" ] && { echo "SKIP $S (done)"; exit 0; } +t0=$SECONDS +for c in chr21 chr22; do + "$BIN" call --subject "$S" --contig $c --out /dev/null >/dev/null 2>>"$SC/runs/$S.log" || { echo "FAIL-call $S $c"; exit 1; } +done +"$BIN" archaic-segments --subject "$S" --json > "$SC/runs/$S.json" 2>>"$SC/runs/$S.log" || { echo "FAIL-seg $S"; exit 1; } +echo "OK $S $((SECONDS-t0))s"