Docs/archaic observable diagnosis - #42
Merged
Merged
Conversation
…k 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>
…ounting 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>
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>
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>
…umber 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>
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>
…ds 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>
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>
…sults 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>
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.
No description provided.