Gate Tier B archaic segments off: no per-individual signal - #40
Merged
Conversation
…g containers
Every CRAM region query cost the same as reading the entire contig, no matter
how small the region. Measured on a 30x 1kGP CRAM, a ONE-BASE query:
chr21 (45 Mb) 20.92 s chr1 (248 Mb) 116.4 s
same query on a BAM: 4-6 ms
Cost tracked contig length exactly (248/45 = 5.5x length, 116/21 = 5.5x time)
and was completely independent of region size -- a 1 bp and a 1 Mb query on
chr21 both took 20.91 s.
CAUSE. noodles' `Query::read_next_container` skips a container only when its
`reference_sequence_id` does not match; it never consults the requested
interval, so it decodes every container of the chromosome and filters records
afterwards. Our own `for_each` had replicated the same
reference-sequence-only filter. The `.crai` already carries `alignment_start`
and `alignment_span` per container -- the information needed to skip was
present and unused. On HG00096, chr21 holds 1,140 containers and a point query
needs exactly 1.
FIX. `cram_container_offsets` selects containers by interval overlap, shared by
both read paths. `query` no longer delegates to noodles: it decodes only the
containers that can overlap, lazily, one at a time so an early-stopping caller
does not pay for the rest. A container the index cannot place
(`alignment_start` = None) is KEPT -- skipping is only ever done on positive
evidence a container lies outside the interval, because a wrongly skipped
container is reads silently missing from a variant call.
RESULT, same file, same machine:
chr21 1 bp 20.92 s -> 8.3 ms (~2,500x)
chr1 1 bp 116.4 s -> 113 ms (~1,030x)
chr21 1 Mb 20.91 s -> 438 ms (~48x)
`navigator call --contig chr21` >96 min -> 44 s (~130x)
The last line is the one that matters: the whole-chromosome cost was paid once
per caller chunk AND once per realignment candidate, which is why a single
chromosome ran over an hour and a half without finishing.
CORRECTNESS. Faster-but-lossy would present as a faster caller rather than a
broken one, so equivalence is tested, not assumed:
- `cram_query_matches_noodles_query` runs our query and noodles' own Query over
the checked-in fixture and compares every field of every record.
- `container_offsets_select_only_overlapping_containers` pins the selection
boundaries (first/last base of a container overlap, one past either end does
not, unbounded intervals keep everything, other references are never
selected, unplaceable containers are kept) -- the fixture is a single
container and cannot catch an off-by-one here.
- On the real 11 GB CRAM, `VERIFY=1 cram_query_probe` compares against noodles
record for record: 305 records over 1 kb and 47,898 over 200 kb, byte
identical on name, position, flags and sequence.
Why this went unnoticed: the ground-truth subject's alignment is a BAM, and its
whole genome calls in ~5 minutes. Every CRAM in the workspace -- including all
3,216 1kGP alignments -- has been paying the whole-contig cost. Found while
trying to validate archaic segment calls against the hmmix 1000G callset, which
needs those CRAMs.
Adds `cram_query_probe`, the harness these numbers come from, as the regression
tool: it attributes cost per phase (open, first query, warm query, bulk region)
so a future slowdown names the part responsible.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Tier B shipped in alpha.14 on the strength of one number -- its total archaic extent landed at 1.01x the hmmix European mean. Validating it against hmmix's own calls FOR THE SAME INDIVIDUALS shows that number is all there is. The workspace turned out to hold CHM13 CRAMs for 2,307 of the 2,309 people in hmmix's published callset, so the comparison is per-person rather than against a cohort distribution. (Reaching them needed the CRAM region-query fix in #39: chr21+22 per sample went from ~3.5 hours to ~90 seconds.) LOCATIONS -- BELOW CHANCE. HG00096, chr21+22, hmmix lifted hg38->CHM13: ours 2.294 Mb / 112 segments hmmix 2.333 Mb / 48 tracts base overlap 0.050 Mb -> sensitivity 2.1%, precision 1.5% null (our own segment lengths placed at random in the same span): 5.0%, p95 9.4% Every alternative explanation was tested and rejected: overlap-vs-shift is flat across +/-2 Mb with no peak (not a coordinate error, and lifted fragment lengths sum to the hg38 input exactly); 70.7% of the truth lies inside our callable territory, and sensitivity restricted to reachable truth is still 3.0% (not callability); unioning haplotypes reproduces hmmix's published 2.09 Mb EUR mean exactly (not haplotype handling); overlap re-derived brute force (not the harness -- though one real harness bug WAS found: CrossMap splits tracts at median 2 bp gaps, inflating 48 tracts to 423). AMOUNTS -- NO CORRELATION. n=20 Europeans, randomly drawn, truth spanning 1.19-2.97 Mb: Pearson r = -0.018 (permutation p = 0.94) Spearman = -0.020 (p = 0.94) mean ratio ours/theirs = 0.923 SD: truth 0.496 Mb, ours 0.312 Mb -> our spread is 0.63x the truth's The two individuals with the LEAST archaic ancestry drew our two highest calls (1.71x, 1.90x); the two with the most drew among our lowest (0.54x, 0.63x). So the caller reproduces the cohort average and nothing about the person -- which is exactly what three fitted parameters were tuned to do. The design recorded that result honestly as "a calibration check, not a validation"; this is what the validation found. THE GATE covers compute, read-back and display, per M3's rule (which the original ship did not follow -- see the design's Deviations section): - `call_archaic_segments_for_subject` returns an ERROR, not an empty result. Every caller asked for a computation, and silently returning zero segments would read as "you have no archaic ancestry" -- a much worse claim than "we are not reporting this". - `cached_archaic_segments` is gated too: rows written before the gate are still in workspaces, and a read gate is the difference between withholding a result and merely declining to recompute one. - The UI card STATES that it is withheld and why, rather than disappearing. A section that silently vanishes between releases reads as a bug; a stated withholding is a finding about the data. - No publish/export path exists for Tier B (M4 never started), so there is nothing further to cover. TIER A IS NOT AFFECTED and is not gated by this: the marker count is direct dosage over a fixed panel with no HMM and no fitted thresholds, checked against the per-site archaic rate on the intersection with real 23andMe v5 chip content. It remains what the Simple-mode card and the Advanced count report. Adds `scripts/archaic-validation/` -- the harness that produced these numbers, since the design now makes re-running it the condition for re-enabling, and it cannot be a gate if it only exists in a scratch directory. Its README records the three traps that produced wrong answers along the way (span-vs-fragment reassembly, the 2 bp lift gaps, and summing haplotypes instead of unioning). Also adds `archaic_callable_dump`, which answered "could we even have found them". 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. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
JamesKane
added a commit
that referenced
this pull request
Jul 31, 2026
…k the instruments (#41) Tier B was gated off (#40) for carrying no per-individual signal. This records WHY, because the caller was not mis-tuned -- it was built on an observable that leaves about one bit of evidence per tract, and no amount of fitting could have rescued it. THE NUMBER THAT SETTLES IT. Detectability of one 36 kb tract (the measured median): evidence sens @5% FP @1% FP private-variant density (built) ~1 variant 14.3% 4.3% archaic-allele matching ~30 sites 95.1% 80.9% Reaching 80% sensitivity at 5% false positives takes 36 kb by allele matching and is not reached by 500 kb with density. Both observables carry the SAME ~3x contrast (2.89x vs 3.04x), so contrast was never the problem: evidence per tract was. THE WRONG TURN. §3 chose Skov 2018 (hmmix), whose premise is detecting introgression WITHOUT archaic reference genomes -- it infers them from private-mutation density precisely because it assumes you do not have them. We do: all four archaic genomes, and archaic_classify already ships 2,031,406 diagnostic sites derived from them. §3 dismissed the reference-based alternative in one clause ("IBDmix is reference-based and heavier"), and that clause is the root of every failure below. An introgressed tract is a haplotype inherited intact; the question is whether a stretch MATCHES an archaic genome, not whether it is slightly more mutated. Runs, not counts. WHAT WAS RULED OUT FIRST, each of which looked right at the time: - the 6.0x hard-coded multiple IS wrong (real enrichment 2.89x), but correcting it changes little; - per-individual Baum-Welch EM -- the thing hmmix does and we do not -- makes it WORSE: unconstrained it learns 22.4x and ~9 kb tracts, calling 7x the truth; - oracle parameters taken from the truth still sit at their random-placement null, precision pinned at 4-5%, so the model CLASS is wrong, not its parameters; - the background is 14.6x overdispersed against the Poisson it is modelled with, with a 5.3x p10-p90 spread -- noise larger than the signal; - a mutation-rate map would not have sufficed: the best proxy available explains 38% of that variance, leaving 7.4x; - quality filtering does not help (it lowers dispersion only by discarding variants); - it is NOT our variant calling: 1000G's own calls for the same person give the same contrast (1.98x vs 2.08x), though ours are 6x noisier -- 2.12x the carried SNVs on chr21 and 97% of private calls unique to us, which is worth its own investigation; - it is NOT the truth set or my lift: in native hg38 with no lifting, hmmix's tracts are enriched 1.84x for their own archaic SNPs against a 1.04x null. WHAT THE REFRAME ALSO FIXES: diagnostic sites become the denominator, so their uneven density cancels and the mutation-rate map is unnecessary; the background's dispersion stops applying because nothing is modelled as a density; the assets already exist; and per-site matching split by lineage class IS Skov's post-hoc annotation, so the gated attribute_lineage may be recoverable with it. STILL OPEN: the background carrying rate measures 13.0% here against the 4.3% recorded earlier in the design -- unreconciled, and it sets the contrast. ~30 sites per tract assumes adequate coverage at diagnostic sites. And the reframe is n=1 so far. Banked: three diagnostic examples (archaic_private_dump, archaic_outgroup_density, archaic_classify_dump) and the eight analysis scripts, with a README mapping each to the question it answered. Two harness bugs that produced confident wrong answers are documented there too -- conditioning on "has a call" (reported ~80% carrying against a known 4.3% background) and measuring tract contrast in 100 kb bins (reported 1.14x where the correct figure is 1.98x). Documentation and dev-only examples; no shipped behaviour changes. Tier B stays gated. Co-authored-by: Claude Opus 5 (1M context) <noreply@anthropic.com>
JamesKane
added a commit
that referenced
this pull request
Aug 1, 2026
* Archaic Tier B: diagnose the failure as the wrong OBSERVABLE, and bank the instruments Tier B was gated off (#40) for carrying no per-individual signal. This records WHY, because the caller was not mis-tuned -- it was built on an observable that leaves about one bit of evidence per tract, and no amount of fitting could have rescued it. THE NUMBER THAT SETTLES IT. Detectability of one 36 kb tract (the measured median): evidence sens @5% FP @1% FP private-variant density (built) ~1 variant 14.3% 4.3% archaic-allele matching ~30 sites 95.1% 80.9% Reaching 80% sensitivity at 5% false positives takes 36 kb by allele matching and is not reached by 500 kb with density. Both observables carry the SAME ~3x contrast (2.89x vs 3.04x), so contrast was never the problem: evidence per tract was. THE WRONG TURN. §3 chose Skov 2018 (hmmix), whose premise is detecting introgression WITHOUT archaic reference genomes -- it infers them from private-mutation density precisely because it assumes you do not have them. We do: all four archaic genomes, and archaic_classify already ships 2,031,406 diagnostic sites derived from them. §3 dismissed the reference-based alternative in one clause ("IBDmix is reference-based and heavier"), and that clause is the root of every failure below. An introgressed tract is a haplotype inherited intact; the question is whether a stretch MATCHES an archaic genome, not whether it is slightly more mutated. Runs, not counts. WHAT WAS RULED OUT FIRST, each of which looked right at the time: - the 6.0x hard-coded multiple IS wrong (real enrichment 2.89x), but correcting it changes little; - per-individual Baum-Welch EM -- the thing hmmix does and we do not -- makes it WORSE: unconstrained it learns 22.4x and ~9 kb tracts, calling 7x the truth; - oracle parameters taken from the truth still sit at their random-placement null, precision pinned at 4-5%, so the model CLASS is wrong, not its parameters; - the background is 14.6x overdispersed against the Poisson it is modelled with, with a 5.3x p10-p90 spread -- noise larger than the signal; - a mutation-rate map would not have sufficed: the best proxy available explains 38% of that variance, leaving 7.4x; - quality filtering does not help (it lowers dispersion only by discarding variants); - it is NOT our variant calling: 1000G's own calls for the same person give the same contrast (1.98x vs 2.08x), though ours are 6x noisier -- 2.12x the carried SNVs on chr21 and 97% of private calls unique to us, which is worth its own investigation; - it is NOT the truth set or my lift: in native hg38 with no lifting, hmmix's tracts are enriched 1.84x for their own archaic SNPs against a 1.04x null. WHAT THE REFRAME ALSO FIXES: diagnostic sites become the denominator, so their uneven density cancels and the mutation-rate map is unnecessary; the background's dispersion stops applying because nothing is modelled as a density; the assets already exist; and per-site matching split by lineage class IS Skov's post-hoc annotation, so the gated attribute_lineage may be recoverable with it. STILL OPEN: the background carrying rate measures 13.0% here against the 4.3% recorded earlier in the design -- unreconciled, and it sets the contrast. ~30 sites per tract assumes adequate coverage at diagnostic sites. And the reframe is n=1 so far. Banked: three diagnostic examples (archaic_private_dump, archaic_outgroup_density, archaic_classify_dump) and the eight analysis scripts, with a README mapping each to the question it answered. Two harness bugs that produced confident wrong answers are documented there too -- conditioning on "has a call" (reported ~80% carrying against a known 4.3% background) and measuring tract contrast in 100 kb bins (reported 1.14x where the correct figure is 1.98x). Documentation and dev-only examples; no shipped behaviour changes. Tier B stays gated. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * Archaic Tier B v2: call tracts by matching the archaic genomes, not counting mutations First increment of the rework diagnosed in the previous commit. New module `archaic_match`, its probe, and the harness fix needed to score it honestly. Not yet wired to the app; Tier B stays gated. THE MODEL. An introgressed tract is a haplotype inherited intact from an archaic ancestor, so it carries the archaic allele at a large share of the diagnostic sites it spans. That is a two-state HMM whose observation is ONE BIT PER DIAGNOSTIC SITE -- carried or not -- with Bernoulli emissions, indexed BY SITE rather than by base pair, and recombination-scaled transitions between consecutive sites. Indexing by site is what makes this robust where density was not: diagnostic sites become the denominator, so their uneven density cancels. The mutation-rate map the density model needed -- and which no available proxy supplied, the best explaining 38% of a 14.6x overdispersion -- is simply not required. MEASURED ON REAL DATA, HG00096 chr21+22, against hmmix's own calls: density (shipped) matching (this) sensitivity 2.1% 38.7% precision 1.5% 19.8% Jaccard 0.009 0.151 tracts hit 4/47 16/47 random-placement null 5.0% 6.5% (p95 12.3, max 20.2) verdict AT OR BELOW CHANCE ABOVE THE NULL'S FULL RANGE 38.7% is 5.9x the null mean and higher than any of 400 random draws. The same data, the same truth, the same harness -- only the observable changed. TWO CHOICES WORTH REVIEWING: - Parameters are MEASURED, not fitted. p_background is estimated directly from the subject's genome-wide carrying rate (dominated by non-archaic sequence, so it estimates the background); p_archaic derives from it by a measured 3.04x ratio. Unconstrained Baum-Welch on the previous caller diverged to a degenerate fit -- 22x emission ratio, 9 kb tracts, 7x the truth -- so EM is deliberately not used here. - Sites where the archaic-derived allele IS the reference base are DROPPED. Every reference-matching genome trivially carries the derived allele there, so the site separates nothing; worse, because the caller emits only variant records, a no-call at such a site means the subject DOES carry it, the opposite of what a no-call means everywhere else. Getting this wrong is not hypothetical: an early version of the analysis conditioned on "has a call" and reported an 80% carrying rate against a known 4.3% background. Five unit tests cover the model without assets or I/O: a real run is found, a background-only contig calls nothing, scattered background-rate carriers do NOT accumulate into a tract (the density caller's defining failure, restated in this observable), reference-derived sites are dropped, and a missing call is a hom-reference non-carrier. HARNESS FIX, without which the numbers above are unreadable: compare_locations.py now computes the random-placement null itself and merges the truth with a 1 kb tolerance. Sensitivity rises with how much sequence you call, so a bare figure means nothing -- the density caller's 2.1% looked like weak performance when it was below chance. The tolerance closes the median 2 bp gaps the lift leaves, which otherwise report 423 shards where there are 47 real tracts. STILL OPEN, and none of it is hidden by the numbers above: we over-call by 1.95x (4.56 Mb against 2.34), precision is 19.8%, and this is n=1. Calibration is now legitimate because the locations are real, but it must be fitted on some individuals and validated on held-out ones -- fitting to a cohort statistic is exactly what produced the last caller. Attribution stays off. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * Archaic Tier B v2: calibrate on 30 Europeans, report on 30 held-out ones Thresholds fitted by grid search on a TRAIN split and reported on a TEST split that never touched the fit. The split exists because the density caller was tuned until a cohort statistic matched and the statistic was then reported as evidence; this makes that failure visible instead of absorbing it. HELD-OUT TEST, 30 Europeans, chr21+22, against hmmix's own calls: density defaults calibrated base-level F1 -- 27.9% 32.5% precision 1.5% 20.2% 28.8% extent ratio ours/theirs 1.45 2.23 1.30 per-individual extent r -0.018 +0.520 +0.642 (p=0.94) (p=0.0001) Locations were already decisive on the earlier 20-sample cohort: 20/20 individuals above their OWN random-placement null, computed at their own called extent, mean sensitivity 45.3% against a 7.1% null. The density caller scored 2.1% against 5.0% -- below chance. Objective was base-level F1, not sensitivity: sensitivity is bought by calling more sequence, and the uncalibrated caller over-called 2.2x while still scoring 45%. F1 makes over-calling cost something. THE ARGMAX WAS NOT TAKEN, deliberately. The grid's best TRAIN F1 wanted min_segment_bp = 40 kb, which discards 61% of real tracts by construction -- hmmix's median tract is 31-36 kb and its p10 is 7 kb. It bought 0.1 F1 points over the 5 kb floor and was WORSE on both things the number exists for: sensitivity (37.4% vs 38.2%) and per-individual extent correlation (+0.642 vs +0.658). The design doc already records this trap once, at 50 kb. Chosen: min_posterior 0.95, min_sites 24, min_segment_bp 5 kb. One reporting bug fixed while writing this: the harness flagged ANY train-test gap over 5 points as overfitting, including the case where TEST scores higher. Test above train is split-to-split variation at n=30, not the grid finding the split. Overfitting is train >> test, and only that now trips the warning. STILL NOT GOOD ENOUGH TO RE-ENABLE, and the module says so: precision is 28.8% and the extent ratio 1.30, so the caller over-calls by about a third. ARCHAIC_SEGMENTS_ENABLED stays false until that closes and the result reproduces outside Europe -- the cohort here is European only, on two chromosomes, and `archaic_ratio` was held at its measured value rather than swept because it changes the emissions and needs a re-run per value. Banks calibrate.py (the split + grid) and cohort_score.py (per-individual scoring against each individual's own null) alongside the existing harness. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * Archaic Tier B v2: sweep the emission ratio; over-calling is gone The previous calibration swept only the three post-hoc thresholds, because archaic_ratio changes the emissions and cannot be re-filtered from a finished run. Sweeping it is what removed the systematic over-calling. The probe now sweeps it IN PROCESS (`ARCHAIC_RATIOS=2.0,2.5,...`): the expensive work -- reading the reference, walking the diagnostic sites -- is per sample, not per ratio, so eight values cost one pass rather than eight. 60 samples x 8 ratios runs in well under a second each. HELD-OUT TEST, 30 Europeans that never touched the fit: density uncalibrated thresholds + ratio base-level F1 -- 27.9% 33.1% 34.5% precision 1.5% 20.2% 31.4% 34.9% extent ratio ours/theirs 1.45 2.23 1.11 0.98 per-individual extent r -0.018 +0.520 +0.628 +0.710 (p=0.94) (p<0.0001) The extent ratio is the line that matters: 0.98 means the caller is no longer systematically over-calling, which no threshold could fix because the emissions were wrong. archaic_ratio is now FITTED, not measured, and the docs say so. The observed enrichment inside real tracts is 3.04x but the model separates best at 4.5x -- not a contradiction: 3.04x is the average over an external tract set that is itself weakly supported, while the emission ratio is what makes the HMM selective enough to place boundaries. p_background remains estimated per individual from the subject's own data, unfitted. THE ARGMAX WAS AGAIN NOT TAKEN. The grid preferred min_segment_bp = 10 kb; 5 kb is chosen. Within the plateau the two differ by 0.1 F1, 5 kb is slightly BETTER on per-individual correlation (+0.710 vs +0.706), and it discards half as many real tracts (8% of the truth under 5 kb against 16% under 10 kb). An earlier sweep wanted 40 kb, which discards 61%; the design records the same trap at 50 kb. Structural exclusion of real tracts is not worth a tenth of a point. One unit test needed pinning rather than fixing: `finds_a_run_of_carried_sites` broke when min_posterior rose to 0.98 and trimmed the run's edges. That is correct behaviour and the wrong thing for a model test to be sensitive to, so it now sets its own thresholds instead of inheriting the calibrated ones. STILL GATED, with specific limits rather than general unease: precision 34.9% means two thirds of called sequence is not in the reference callset; the cohort is European only and chr21+22 only; and the reference callset is itself weakly supported (hmmix's tracts are enriched just 1.84x for their own archaic SNPs), so agreement with it caps well below 100% even for a correct caller. East Asians are the sharp next test -- the truth predicts ~1.18x more archaic sequence there, which a caller merely tracking European structure would miss. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * Archaic Tier B v2: detection transfers to East Asians; the reported number does not Ran 30 East Asians with the parameters FROZEN at the European fit -- nothing refitted -- because a caller tuned on one population reproducing that population's statistics proves nothing. DETECTION TRANSFERS: Europe (fitted) East Asia (new) above own null 60/60 30/30 sensitivity 31.6% 31.6% precision 32.2% 41.9% per-individual r +0.620 +0.545 Identical sensitivity and BETTER precision on a population the thresholds never saw. The calibration learned archaic structure, not European structure. THE REPORTED EXTENT ORDERS THE POPULATIONS BACKWARDS, and this is now the single reason the feature stays gated. Truth puts East Asian archaic extent at 1.217x Europe's; we call 0.937x. A user would be told an East Asian carries LESS archaic ancestry than a European. Cause: reported extent is true positives plus false positives, and the false-positive load is population-dependent (precision 32.2% vs 41.9%), so Europeans accumulate more spurious extent. A RESULT I ALMOST REPORTED AS A SUCCESS, recorded so it is not quoted later: "detected sequence reproduces 1.22x" is CIRCULAR. detected = sensitivity x truth, and sensitivity is equal across the two populations, so that ratio matches by construction. It restates the invariance; it does not test the ordering. The honest read is that per-individual detection transfers and the population-level number does not. THREE CAUSES RULED OUT, each measured rather than argued: - background contamination of p_background -- carrying rates 11.9% (EUR) vs 12.2% (EAS), and both emission states scale together so the contrast is preserved; - tract length -- median 29 kb in both populations; East Asians simply have MORE tracts (54 vs 46 per person), consistent with the 1.22x truth; - panel ascertainment -- in-tract contrast 2.99x (EUR) vs 3.04x (EAS), ratio 1.014, so the diagnostic panel is equally informative in both. This was the most plausible candidate and it is not the cause. Banks cross_population.py (the frozen-parameter transfer test, with a permutation test on the population difference) and observable/ascertainment.py (the per-population contrast measurement). Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * Archaic Tier B v2: an arbiter that does not ask another caller's opinion Precision has been scored against hmmix, but a call they did not make is not automatically wrong -- their callset is incomplete by an unknown amount (their own tracts are enriched just 1.84x for their own archaic SNPs). Settling that needed a referee independent of both callers. The Tier A panel records, per site, which of the four archaic genomes carries the derived allele. The Tier B caller NEVER SEES THIS -- it reads only a derived base and a lineage class from ArchaicClassify -- so per-genome concordance is evidence it cannot have been fitted to. Of the sites where a given archaic genome is derived, what fraction does the subject carry: true positive FALSE positive background Europe 93.6% 81.3% 59.0% East Asia 93.5% 72.9% 45.5% Our "false positives" sit 64% (EUR) and 57% (EAS) of the way from background to true positive. They are a MIXTURE: real tracts hmmix missed, genuine noise, and calls that are correctly placed but over-extended. So precision against hmmix understates this caller -- though not enough to dismiss it, and F1 stays a usable objective, which means the calibration was not invalidated. THREE ATTEMPTS, TWO OF THEM NON-ANSWERS, all kept because each looked convincing: - carrying rate inside segments: CIRCULAR. The caller selects on carrying rate, so every segment it emits has a high one. It reported false positives at 98% of the way to true positive, which is meaningless. - overlap with other individuals' tracts: SATURATED. Across 90 individuals the hmmix tracts cover 67% of callable territory, so a random region hits them ~60% of the time and the observed excess was 2-7 points. - the first version of the arbiter itself: VACUOUS. It conditioned on the subject already carrying the derived allele, but at a discordant site at least one genome is derived by construction, so best-genome concordance was ~100% EVERYWHERE INCLUDING BACKGROUND. Fixed by conditioning on the GENOME instead -- of the sites where genome g is derived, how many does the subject carry -- at which point background correctly lands near the genome-wide carrying rate. A LEAD ON THE ORDERING INVERSION: background concordance differs by population, 59.0% in Europeans against 45.5% in East Asians. Europeans carry archaic-derived alleles more often OUTSIDE tracts, which is a plausible mechanism for the population-varying false-positive load and hence for the reported extent ordering the populations backwards. Adds archaic_panel_dump (the per-genome calls as TSV) and banks all three attempts under scripts/archaic-validation/, since the two failures are the more instructive half. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * Archaic Tier B v2: concordance filter takes precision to 90%, and finds a harder limit Scoring each called segment against the archaic genomes and dropping poor matches raises PRECISION FROM 54% TO 90%. The filter is sound on a held-out check: with Denisova excluded from it entirely, kept segments score 74.9% on Denisova concordance against 21.5% for dropped ones -- a 3.5x separation on a genome the filter never saw. IT DOES NOT FIX THE POPULATION ORDERING, and tightening it makes the ordering worse (0.927 unfiltered -> 0.895 at 90% precision -> 0.831 at the strictest threshold, against a 1.098 target on this subset). At 90% precision the reported extent is mostly true positives, so false positives are no longer the explanation. What remains is recovery: ~46% of European truth against ~38% of East Asian. THE CAUSE, and it is not a threshold: EUR true positives Neanderthal-best 89.2% Denisova 69.3% EAS true positives Neanderthal-best 83.4% Denisova 77.4% Denisova beats every Neanderthal: EUR 11.2% EAS 32.2% That 2.9x is the Denisovan ancestry East Asians carry and Europeans essentially lack -- the data reproduces known biology, which is reassuring about the observable. But it also means our four sequenced archaic genomes UNDER-REPRESENT EAST ASIAN ARCHAIC DIVERSITY, so any reference-based filter under-calls East Asians. Fixing it would take archaic genomes closer to the populations that introgressed into East Asia, which do not exist. So a cross-population comparable number is not achievable this way. The caller is defensible WITHIN a population and not BETWEEN them, and the module says so. HOLDING DENISOVA OUT WAS THE WRONG CHOICE, for a reason worth recording: it is not a neutral held-out genome. Excluding it biases the filter against exactly the population that carries Denisovan ancestry, and the best-of-three ordering (0.824) is worse than best-of-four (0.895). A held-out validator has to be neutral with respect to the thing being measured, and this one was not. Also records why p_background cannot fix the ordering: the model's estimate is flat across populations (14.17% EUR vs 14.04% EAS, ratio 0.991) while the false-positive load is not (ratio 0.877), so the parameter is blind to the driver. Banks concordance_filter.py (the filter with the held-out design) and observable/ordering_cause.py (the p_background diagnosis). Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * Ship Tier B archaic segments as a WITHIN-POPULATION measure Re-enables ARCHAIC_SEGMENTS_ENABLED on the rebuilt caller, wires it through the app, and states the comparability limit where the number is read. WHAT EARNED IT BACK. Held out on 30 Europeans the fit never saw, per-individual extent correlates at r = +0.710 (p < 0.0001) where the withdrawn density caller managed -0.018 (p = 0.94). All 90 individuals across two populations score above their own random-placement null. The concordance filter takes precision from 54% to 90%, verified with Denisova held out of it entirely: kept segments score 74.9% on Denisova concordance against 21.5% for dropped ones. WHAT IT MUST NOT BE USED FOR, and why the caveat is in the UI rather than a doc: East Asian tracts match our four sequenced archaic genomes less well than European ones (83.4% vs 89.2%), so their extent is under-called and the reported figure orders the two populations backwards against the truth. That is a property of which archaic genomes have been sequenced -- Denisova is the best match for 32.2% of East Asian tracts against 11.2% of European ones -- not a threshold to tune. The caveat sits directly under the number, in amber, because a limit recorded anywhere else will not be read by the person reading the number. CHANGES: - archaic_match gains segment_concordance / carried_panel_sites / filter_by_concordance, with MIN_CONCORDANCE = 0.70 from the measured plateau. - The app method now runs matching + filter instead of the density caller, and loads the Tier A panel for the per-genome calls. The outgroup asset is no longer needed here: it existed to strip African-shared variants for a density model that no longer exists. - The reference is read one contig at a time -- whether a diagnostic site is informative depends on the reference base there, and holding all of CHM13 would cost 3.1 GB. - archaic_segment_sig now includes METHOD_VERSION. Without it a workspace holding results from the withdrawn caller would keep serving them: stale answers from a method removed for having no per-individual signal. Three of the eight tests pin failures that actually happened during development: concordance must condition on the GENOME not the subject (the other way scores ~100% everywhere including background, separating nothing); a segment with too few judgable sites is KEPT, since absence of evidence is not evidence of a bad call and dropping on it would quietly penalise sparse regions; and the summary is recomputed after filtering rather than carried over. VALIDATED ON chr21+chr22 ONLY. This ships genome-wide. The design records the same extrapolation trap from the previous caller -- projecting the chr21+22 target predicted 85.1 Mb against a measured 90.9, 6% low -- so this is a real gap, not a formality. A genome-wide run is the next thing to do, and the inputs for it (whole-genome calls on the ground-truth subject, hmmix's genome-wide European callset) are already local. Verified: clippy clean under -D warnings; 30 suites pass including i18n parity; `navigator archaic-segments` runs end to end through the app path. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * Genome-wide validation passes; fix a cache key that served partial results forever Ran the validation the previous commit called for, and it found a real bug before it produced a number. THE BUG. `archaic_segment_sig` keyed on alignment + genotype version + method version, but NOT on which contigs were actually called. So a result computed from a partial call set stayed "current" indefinitely. Caught in the act: HG00096 had all 23 contigs cached and the app cheerfully served 33 segments over 2 contigs (1.94 Mb) from an earlier chr21+22 run. Anyone who called one chromosome, then later ran the whole genome, would have kept the one-chromosome answer. Fixed by folding the called contigs into the key. Reading them needed a cheap query -- `list_for_alignment` returns payloads, which here is ~1 GB of JSON across 22 contigs -- so `artifact::list_kinds` returns kinds alone. With the fix HG00096 recomputes to 1,208 segments over 22 contigs (83.56 Mb). GENOME-WIDE VALIDATION, three Europeans against hmmix's genome-wide callset: ours hmmix ratio sens prec null (max of 400) HG00096 83.6 Mb 93.0 0.90 40.3% 44.9% 5.5% HG00102 83.9 Mb 89.3 0.94 42.4% 45.1% 4.9% HG00112 82.1 Mb 91.0 0.90 42.9% 47.5% 5.1% All three above the ENTIRE random-placement null. Both sensitivity and precision come out BETTER genome-wide than on chr21+22 (40-43% against 31.6%, ~46% against 34.9%), so the two-chromosome numbers were conservative rather than optimistic. That direction matters. The design records the previous caller being burned by the opposite -- a chr21+22 target that extrapolated 6% low -- which is exactly why shipping on a two-chromosome result was worth objecting to even though the result here turned out fine. The objection was right; the answer was just favourable. Extent lands at 0.90-0.94 of hmmix's, i.e. we under-call by about 10% genome-wide rather than over-call. Module docs now lead with these figures. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> --------- Co-authored-by: Claude Opus 5 (1M context) <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Tier B shipped in alpha.14 on the strength of one number — its total archaic extent landed at 1.01× the hmmix European mean. Validating it against hmmix's own calls for the same individuals shows that number is all there is.
Locations — below chance
HG00096, chr21+22, hmmix lifted hg38→CHM13:
Every alternative explanation was tested and rejected:
Amounts — no correlation
n = 20 Europeans, randomly drawn, truth spanning 1.19–2.97 Mb:
The two individuals with the least archaic ancestry drew our two highest calls (1.71×, 1.90×); the two with the most drew among our lowest (0.54×, 0.63×).
The caller reproduces the cohort average and nothing about the person — which is exactly what three fitted parameters were tuned to do. The design recorded that result honestly as "a calibration check, not a validation". This is what the validation found.
The gate
Covers compute, read-back and display, per M3's rule — which the original ship did not follow (see the design's Deviations section):
call_archaic_segments_for_subjectreturns an error, not an empty result. Every caller asked for a computation; silently returning zero segments would read as "you have no archaic ancestry", a much worse claim than "we are not reporting this".cached_archaic_segmentsis gated too — rows written before the gate are still in workspaces, and a read gate is the difference between withholding a result and merely declining to recompute one.Tier A is not affected
Not gated by this. The marker count is direct dosage over a fixed panel — no HMM, no fitted thresholds — and was checked differently, against the per-site archaic rate on the intersection with real 23andMe v5 chip content. It remains what the Simple-mode "Neanderthal ancestry" card and the Advanced count + percentile report.
Also
scripts/archaic-validation/— the harness that produced these numbers. The design now makes re-running it the condition for re-enabling, and it can't be a gate if it only exists in a scratch directory. Its README records the three traps that produced wrong answers along the way: span-vs-fragment reassembly, the 2 bp lift gaps, and summing haplotypes instead of unioning. Plusarchaic_callable_dump, which answered "could we even have found them".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.
Clippy clean under
-D warnings; 30 suites pass including i18n parity.Depends on #39 for the CRAM fix (not required to merge, but the validation isn't rerunnable without it).
🤖 Generated with Claude Code