The manifest was lying, and so was I
2026-09-22 → 09-23. Staging a reference genome, finding two bugs on the way, and discovering that a number we’d been steering by had been wrong for two months.
Background: the PRS primer §2.5 on the pinned reference genome, and the genomics glossary §4.1 on where it comes from.
The job was supposed to be a download
Gate A of the PRS rebuild needs the human reference genome on local disk — the pinned Ensembl release-109 GRCh38 primary assembly, 3.15 GB uncompressed. There is a script. It downloads, verifies a checksum, decompresses. One command.
It was going to take two and a half hours.
ftp.ensembl.org was serving about 102 KB/s. That is not a broken link and it is not a
big file problem — 840 MB compressed should not take an afternoon. The tell is that the
number is suspiciously round and stable. Servers that are genuinely saturated give you a
jagged rate. Servers that are throttling give you a flat line.
The 30-second experiment: does it honour HTTP Range, and do parallel connections scale?
HTTP/1.1 206 Partial Content
Accept-Ranges: bytes
Content-Range: bytes 0-1023/881211416
Six parallel range requests, while the slow single stream was still running: 1,585 KB/s aggregate. A 15× speedup, and running both at once is what proves the cap is per-connection rather than shared.
The real fetch took nine minutes.
The lesson isn’t “use more connections.” It’s that a multi-hour ETA deserves thirty seconds of scepticism before you accept it, and that the fix belongs in the committed script — not in a shell command that dies with the terminal. A fresh machine, a worker image, or a future session will run that script too.
One design point worth stating because it’s what makes the parallel path safe: correctness does not rest on the reassembly. The uncompressed sha256 is still checked against the pin. A mis-ordered or short assembly fails loudly rather than quietly becoming the genome you score a patient against.
Then the reference genome vindicated itself
With the FASTA staged (194 contigs, checksum matching, reference bases correct at three known loci), the real test is whether it works in situ. So: a Gate A dense-callset build against Maureen’s actual 30× germline gVCF, using a union of seven colorectal-relevant scoring files.
3,338 PGS loci + 6 sentinels = 3,344 targets
1,464 homref + 1,791 variant + 78 missing + 11 structural = 3,344 reconciles
All six sentinel loci reproduced July’s result exactly.
1,464 loci — 43.8% — are callable hom-ref. Every one of those scored zero under the old variant-only store. That is the REF-effect bug, the one that invalidated the entire legacy score fleet, demonstrated across thousands of real loci instead of six sentinels.
But the interesting part was the 11 structural failures. 0.33% is small. At catalog scale it’s ~40,000 loci, and “structural failure” folds together two causes with opposite implications: a REF-vs-FASTA disagreement would implicate the reference I had just staged, whereas unresolvable indel orientation is a known and harmless case.
So I itemised them. Two were genuine strand flips — the complement of a scoring allele equalled the FASTA base, an author-strand submission that harmonisation missed. At 2 of 3,338, or 0.06%, that’s the forward-strand premise the whole scorer rests on, confirmed empirically. The primer had predicted exactly this: “near zero confirms the premise; a large count refutes it.”
The FASTA was fine. The other nine were something else entirely.
Nine loci that were not what they claimed
PGS003436 11:74704207 effect_allele='G' other_allele='' hm_inferOtherAllele='A/T'
A/T is not a sequence. It is the harmonizer saying “the other allele is A or T — I
can’t tell which.” An ambiguity marker.
The classifier tested alleles by length:
if (ea and len(ea) != 1) or (oa and len(oa) != 1):
self.indel += 1len('A/T') == 3. So it counted as an indel.
That would be a cosmetic mislabel except for one thing: indel is a promoted Class-B
category — we’ve signed off that we can compute scores containing indels — and ambiguity
deliberately is not. So rows we cannot even orient were helping scores qualify for
production.
Measured over 2,962 scores: 223 carried these tokens, 132 were production-eligible, and
for all 132 n_indel == n_ambiguous exactly. Not one of them contained a single real
indel. Their entire indel signal — and with it part of their eligibility — was an artifact.
Length-1 tokens were worse. N (unknown base), -, P, I, D all have length 1, so
they passed every test and classified as clean Class A. Invisible.
Before writing the fix I scanned the catalog for what allele tokens actually occur, rather than assuming. Good thing:
| Column | Distinct non-ACGT tokens | Examples |
|---|---|---|
hm_inferOtherAllele | 11 (36,476 rows, 220 scores) | A/C, G/T, and three-way C/G/T, A/C/G |
other_allele | 23 (726 rows, 352 scores) | N, -, P, I, D, GTT,AT, APOE_E3, HLA-B_67_S |
A slash-only predicate — the obvious fix from the bug report — would have missed the comma lists, the unknown bases, and the star-allele names.
The fix is one category, ambiguous_other, un-promoted: an allele is resolvable only if it
matches ^[ACGT]+$. Deciding which of those four situations deserves which treatment is a
vocabulary-mapping call, and the house rule is not to make those unilaterally — so it went
to Richard as four options, and was filed as a follow-up once he picked one.
Then the same raw-length test turned up in a second place, where it was not cosmetic at
all: the locus-union builder sets is_indel, and the dense-callset builder reads that flag
to select stricter GQ/DP thresholds. A/T was being gated as though it were a noisy indel
call.
The number we’d been steering by
Re-classifying the catalog under the new rules took 103 minutes across six shards, zero failures. Then the result:
| before | after | |
|---|---|---|
| Class A | 1,997 | 1,077 |
| Class B | 2,950 | 3,870 |
| production-eligible | 3,863 | 2,322 |
| eligible variant rows | 1.89 B | 0.89 B |
Eligibility fell by 1,541 scores. I was about to write that up as the impact of the ambiguous-allele fix.
It mostly isn’t.
| Deferred category | before | after |
|---|---|---|
| Unplaceable rows | 18 | 1,437 |
ambiguous_other | — | 144 |
| Multiallelic | 588 | 132 |
The sole-blocker ranking inverted. And the cause is dated: a separate fix — stopping the parser falling back from the harmonized coordinate to the author’s original-build coordinate — merged on 2026-07-26, five days after the catalog was last classified on 07-21.
The manifest had been silently stale for two months. Every decision made from it since July was made on numbers that no longer described the code. Including the standing strategic advice, written in the primer, that “multiallelic handling is by far the highest-leverage open decision.” It isn’t, and hasn’t been since July.
There’s a process lesson underneath: a classifier version is stamped on every row precisely so a re-run re-classifies rather than skipping. That mechanism worked. What was missing was anything that noticed the stamp had gone stale relative to the code that produces it.
Eligibility isn’t random
With a fresh manifest, some structure shows up. Every eligible score carries an ontology term, so grouping is free:
| Namespace | Eligible | Deferred | Elig % | Meaning |
|---|---|---|---|---|
| MONDO | 1,165 | 985 | 54% | disease |
| EFO | 847 | 1,468 | 37% | trait / measurement |
| OBA | 279 | 521 | 35% | biological attribute |
| HP | 29 | 69 | 30% | phenotype |
Disease scores survive better than measurement scores. And by area, cancer does best at 62% eligible, against blood/haematology at 27%.
The driver turns out to be size, and it is not monotonic:
| Score size | Eligible % |
|---|---|
| < 100 variants | 71% |
| 100 – 1k | 73% |
| 1k – 10k | 44% |
| 10k – 100k | 21% ← trough |
| 100k – 1M | 28% |
| > 1M | 30% |
Small curated hit-lists — the pre-2015 style, a few dozen genome-wide-significant SNPs with rsIDs — are clean. The modern LD-aware genome-wide posteriors are messier. Cancer PRS skew old and small, which is why they survive, and that is a fact about publication history rather than about cancer.
What the scores are actually for
I closed the analysis by noting that 44 colorectal cancer scores are production-eligible, as though that were the payoff. It isn’t, and Richard’s correction is worth more than the number:
We already know she has colorectal cancer, and we know quite a lot about it. Calculating those scores can only tell us, at best, if she was lucky or unlucky in getting it. We aren’t interested in any one particular PRS score, or one particular disease. We are using these as inputs into an ML routine, so maximizing the number of data points we can use, without regard to what they detect or predict, is the main goal.
This was already written down — the primer’s own section on why Maple computes these says “breadth beats depth” and “a wrong number is worse than a missing number.” Having it documented did not stop me reaching for the clinically legible number.
It also changes a pending decision. unplaceable_rows now blocks 1,437 scores, and
promoting it means computing a score whose loci are only partly placeable. Clinically that
is unacceptable — you would be reporting a percentile you cannot stand behind.
But as an ML feature it’s different, and the difference is mechanical. The score enters the PDR2 loss multiplied by a learned per-disease weight. A roughly constant deflation from missing loci is largely absorbable by that weight — the model can rescale a feature. Fewer loci means more variance, not more bias, and a feature tolerates variance far better than bias. (⚠ Wrong as stated, corrected the next day: the released loss L2-normalises the weight vector, so the model cannot rescale a feature. What absorbs a constant deflation is the normalisation against the reference panel, Gate D. The conclusion below stands; the mechanism is the PRS primer’s §1.4.)
Which means the right form of the decision is not promote-or-defer. It’s a coverage threshold, set empirically. And the instruments to set it already exist: a coverage payload that records effective variance-weighted information, and a dual-scorer reconciliation that isolates exactly how much a partially-covered score moves when the missing dosages are filled in.
The threshold stops being a judgement call and starts being a measurement. Which is the direction everything here wants to move.
Ledger
Shipped: parallel ranged FASTA fetch with resume (2.5 h → 9 min); the reference genome
staged, checksum-verified and indexed; ambiguous_other as an un-promoted category, in both
places the raw-length test lived; classified_at refreshed on re-classify; the dashboard
taught to render the new category; the catalog re-classified, 5,377 scores, zero failures.
Found along the way: a 2-month-stale manifest; a strategic recommendation that had quietly expired; 591 scores whose entire indel signal was an artifact; 920 scores that had been classified “clean” while containing alleles we cannot resolve.
Left open, deliberately: splitting ambiguous_other into its four real sub-cases
(filed, not critical path); and the partial-coverage threshold, which is now a measurement
problem rather than an argument.
Still true: zero released polygenic scores. Gate C computes to quarantine; D, E and F remain. The feature vector that all of this exists to build has not been built yet.