The decision that wasn’t one
2026-09-24, later the same day. The plan was to take E1’s measurements back to Gate D. Two of the three things E1 measured said something other than what the docs predicted.
What we set out to do
This morning’s fleet run computed 2,321 polygenic scores and wrote them all
quarantined. The documented next step was Gate D — normalize those raw sums
against a reference population — then re-run Gate E under the final policy.
The primer was explicit that E1 was a measurement run whose numbers would settle the open policy questions. That framing turned out to be exactly right, and slightly uncomfortable, because two of the settled questions resolved against what the document itself claimed.
The blocker that was 0.005 % of a file
docs/prs_primer.md called the unplaceable_rows policy “the highest-leverage
open decision by an order of magnitude” — 1,437 scores deferred because some
rows lacked a harmonized genomic position. It was flagged as a judgement call
about partial coverage, the kind CLAUDE.md says to bring to Richard rather than
decide alone.
So I measured it before bringing it, which is how it stopped being a judgement call. Those 1,437 scores are blocked by a median of 0.005 % of their rows. p95 is 0.138 %. The worst is 5 %. Between them they carry 1.27 billion variant rows — more than the entire eligible set they were being excluded from.
The threshold sweep is the part worth keeping:
| threshold | scores admitted |
|---|---|
| 0.1 % | 1,338 |
| 1 % | 1,431 |
| 5 % | 1,437 |
The curve is flat. Any defensible number admits essentially all of them, which means the threshold is not buying coverage — it exists only to refuse a score with a genuinely meaningful hole. We set it at 1 % (a score must be ≥ 99 % placeable) and the eligible set went 2,322 → 3,753 scores, 894.5M → 2.166B rows. Nothing lost eligibility.
What made it safe rather than merely attractive is a property I had to go
looking for: unplaceability is a defect of the scoring file, not of
Maureen’s genome. It deflates her and the reference panel identically, and a
percentile is a monotone transform over that deflation, so Gate D absorbs it
exactly. The cost is a noisier estimate of the same liability — variance, not
bias. That is the tolerable direction for a model feature, and it is a different
argument from the one the primer used to make, which leaned on the model’s
learned weight absorbing it. That argument was retracted in an earlier revision
because losses.py L2-normalizes the weights and can’t rescale a feature. The
normalization absorbs it. The model does not.
Reclassifying without re-reading 59 GB
v4 changes the promotion policy and no feature extraction, so every structural count already in the manifest is still correct and eligibility is a pure function of columns the table holds. A full re-classify would re-read 59 GB of scoring files to arrive at numbers identical to the ones already stored.
That is an argument, not a proof, and the manifest was silently stale for two
months the last time someone reasoned this way. So --verify N re-parses N
scoring files end to end and checks that a real classification agrees with the
recomputed row on class, categories and eligibility — biased toward the scores
whose answer actually changes, since verifying only unchanged rows would pass
without testing the new policy at all. 30 of 30 agreed, zero mismatches. Then
apply: 5,377 rows, ~90 seconds, instead of ~103 minutes.
The number this document had wrong
The primer says palindromic SNVs cost “1.56 % of rows overall, median 2.47 % per affected score.” The aggregate is right — E1 deferred 13,929,705 rows against 13,930,267 predicted, which is a satisfying confirmation of the manifest.
The per-score median is 10.24 %. Four times what was written. p95 is 21.49 %, and 1,088 of 2,321 scores lose more than 5 % of their rows.
Splitting by score size shows why, and inverts the reading:
| score size | scores | mean % palindromic |
|---|---|---|
| < 1k rows | 1,193 | 9.21 % |
| 1k–100k | 452 | 4.60 % |
| 100k–1M | 293 | 0.00 % |
| > 1M | 383 | 0.87 % |
A few hundred genome-wide scores hold most of the rows and carry almost no palindromes. So “1.56 % of rows” is a statement about those scores, not about the fleet. The burden lands on the small curated hit-list scores — the interpretable ones — and every one of them counts as exactly one feature, the same as a nine-million-row score does.
The document had already been wrong here once in the other direction, from a 60-score sample, and the correction over-corrected. The durable lesson is not “check your numbers.” It is that a row-weighted aggregate and a per-item median answer different questions, and you have to weight the statistic by whatever actually consumes the thing. A feature vector consumes one value per score. I have filed this next to the sharding bug from this morning, which was the same error wearing a different hat: slicing by file count when rows were the cost.
The practical consequence is that Gate D’s palindrome recovery is worth considerably more than the docs implied — and only 192 of 2,321 scores (8.3 %) carry an effect allele frequency of their own, so the scoring files cannot resolve this themselves. The reference panel is the only route, which turns Gate D from “normalization, eventually” into the thing standing between us and about a tenth of the small scores’ content.
Gate D is not Gate E times two thousand
The primer describes Gate D as “itself a fleet run over thousands of panel genomes, so it consumes Gate E’s machinery.” Read literally that is 2,504 × the E1 run: ~4.55 GB of dense callset each, ~4.3 hours of scoring each. Eleven terabytes and ten thousand CPU-hours.
It is also unnecessary, because a reference panel is not shaped like Maureen’s data. Her gVCF is one sample with reference blocks, and the hard problem is telling hom-ref from not-sequenced — which is what Gate A’s FASTA-verified callset exists to solve. A panel arrives as a multi-sample matrix: one row per site, every sample’s genotype already explicit. Nothing needs reconstructing.
So the panel pass is one stream over the union with an accumulator of
n_samples × n_scores — about 46 MB for 2,504 × 2,322 — rather than one pass
per genome. Better still, it factors as a sparse matrix product: dosages
(samples × loci) against weights (loci × scores), chunked by locus block so
neither side has to be resident. The per-row Python overhead that set E1’s pace
disappears into a C-level matmul.
What must not change is who does the scoring. It has to be MODS’s scorer over
that panel, not pgsc_calc’s, because a percentile is only meaningful if the
patient and the panel were deflated the same way — and pgsc_calc mean-imputes
missing dosages to 2p and resolves palindromes by frequency where Gate C drops
and defers. Scoring the panel under a different policy than the patient
reintroduces precisely the bias Gate D exists to remove.
Measure the 153 MB file before downloading the 14.9 GB one
The panel we chose is pgsc_HGDP+1kGP_v1 — 4,150 samples, 80 populations, and
the same panel the dual-scorer witness uses, so a MODS-vs-pgsc_calc divergence
stays attributable to arithmetic rather than to panel choice.
It is 14.86 GB, and there was an open risk under it: the packaged panel is a variant subset of unknown extent, and if it does not carry our union loci it is the wrong substrate no matter how convenient. The PGS Catalog publishes the panel’s variant list separately as a 153 MB sorted BED. Both it and our union shards are sorted by contig and position and both use Ensembl naming, so the overlap is a linear merge — minutes and 153 MB, against hours and 14.9 GB.
Measure first. The habit is the same one that wrote check_worker_env.py this
morning after a missing pyfaidx killed a three-hour union build: find out what
the expensive step needs before paying for it.
One thing the small download taught us for free: ftp.ebi.ac.uk served it at
~180 KB/s on a single connection, the same per-connection throttling as
ftp.ensembl.org. At that rate the panel itself would take thirty hours
single-stream, so the fetcher reuses the ranged parallel downloader already
written and tested for the reference FASTA.
Two numbers that were not both right
The panel covered 94.77 % of our locus union and 100.0000 % of every one of 44
sampled scores. I spent hours reconciling that, and the reconciliation was
elegant: the union counts each distinct locus once while a per-score figure is
row-weighted, so rare variants — negligible per score, plentiful in a
distinct-locus count — would make the union over-weight exactly what a panel of
common variants lacks. I falsified three competing hypotheses to get there
(positional clustering, concentration in large scores, %.2f rounding), then
ran an attribution scan that found nine scores in one accession band explaining
96 % of the gap, with PGS005353 at 77 % coverage. I wrote up a design
requirement off the back of it: Gate D needs a per-score coverage flag, because
one score’s distribution would rest on three quarters of its loci while looking
normal in the manifest.
Almost none of that was real. I had measured against
pgsc_all+hm3_variants_hg38_latest.txt.gz — a 153 MB file in the panel’s own FTP
directory. It is not the panel’s variant list. It holds 33.8M variants; the
panel’s .pvar holds 84,325,213. Real union coverage is 98.93 %, uniform
across every contig, and PGS005353 is 99.65 % covered. There is no outlier
and nothing needs flagging.
I picked that file because its name looked right and because a size calculation seemed to confirm it — 16 GB over 4,150 samples at a quarter-byte per genotype is about 15M variants, near enough to 33.8M to feel corroborated. Two weak signals agreeing is not verification. It is one signal.
What catches this kind of thing is reading the artifact. inspect_reference_panel.py
did it on its first run, in one line: 84,325,213 variants over 24 contigs. I had
written that script an hour earlier to check something else entirely — whether the
panel’s REF alleles agree with our pinned FASTA, which they do, 20,000 for 20,000.
The part that stings is the shape of the error. I never tested the cheapest hypothesis available: that the thing I was measuring against was wrong. Four increasingly sophisticated explanations for a discrepancy, and the discrepancy was an instrument artifact. An unexplained gap is at least as likely to be a broken denominator as a finding, and “check what you’re comparing to” belongs before “explain the comparison.”
The attribution technique itself was sound and I’d use it again — a 44-of-2,322 sample has about a 98 % chance of missing a 9-item cluster, so enumerating the missing items and walking their sources genuinely was the right move where sampling harder was not. It was simply aimed at a phantom.
Where this leaves things
The locus union is rebuilding now at 3,753 scores — roughly eight hours, running
on local-prod5, heartbeating. Because the worker resolves eligibility from the
database at runtime rather than from code, it picked up the expanded set without
a deploy, which also means no worker restart, which is the third way I killed a
job this morning.
Gate D proper — panel staged, panel scored, PCs fitted — is the next build, and then E2 re-runs Maureen under the final policy with palindromes resolved. (As it turned out, E2 ran before Gate D with palindromes still deferred; the packaged panel holds 3,942 samples, 3,330 after the relatedness exclusion, not 4,150; and union coverage at the v4 scope this rebuild was for measured 90.7 %, where 98.93 % below is the v3 figure. The PRS primer’s §5.7 has the corrected numbers.) The sequencing matters more than it looks: the panel extraction is keyed to the union, so expanding eligibility after staging the panel would mean staging it twice. That is the whole reason today went to a policy decision instead of straight to Gate D.
Still zero released scores, which remains correct. Gate C writes quarantined
and only Gate F promotes.