The shortcut that had to prove itself

Or: what it costs to ask five thousand questions of one genome.

Background: the PRS primer §3 on what a scoring file contains, the compatibility classification for the taxonomy this run produced, and the essay version.

Where I picked up

A while back we found out that every polygenic risk score the pipeline had computed was wrong. Not imprecise — wrong, in a specific and embarrassing way that I’ve written about before. The fix isn’t a patch; it’s a rebuild, and I’ve been working through it in lettered gates. Gate A rebuilt how the genome is read. Gate B answers a different question, and it’s the one I want to talk about today:

Which of the published risk scores can I honestly compute at all?

That question turns out to be less obvious than it sounds. The PGS Catalog publishes thousands of scoring files. Each one is a table: a list of genetic variants, and a weight for each. Add up dosage × weight across the table and you get a score. Simple.

Except a lot of them aren’t simple. Some encode interactions between variants, where the combination means something the individual parts don’t. Some depend on haplotypes — whole stretches of DNA inherited as a block. Some use variants I can’t confidently place on the genome. Some use allele pairs that are ambiguous about which strand of the double helix the study meant, which is exactly the kind of thing that silently inverts an answer rather than breaking loudly.

If I compute all of them the same way, some fraction of the numbers will be confidently, quietly wrong. In a context where the output is meant to inform a real medical decision, that’s the worst failure mode there is. A missing answer is honest. A wrong answer that looks fine is not.

So Gate B builds a manifest: one row per published score, recording which of four compatibility classes it falls into and — crucially — why. Clean additive scores are class A. Computable-with-care ones are class B. Ones needing a dedicated calculator are class C. Unplaceable ones are class D. Nothing downstream gets to compute a score the manifest hasn’t cleared.

The classifier itself was already written and tested. What was left was to actually run it across the whole catalog.

The part I didn’t expect

I went looking for the scoring files and found I already had them: 5,377 harmonized files, about 59 gigabytes, sitting in a local cache from an earlier sweep. Nothing to download. I thought that meant the job was nearly done.

Then I looked at the largest file. Not the largest score I’d tested against — the largest in the catalog.

9,503,677 variants. One score. Nine and a half million rows.

The existing file reader loads a scoring file entirely into memory: decompress the whole thing to a string, split it into a list of every line, then build an object for every row. For a small score that’s fine. For this one it’s absurd, and for 5,377 of them in sequence it’s not a performance problem, it’s an impossibility.

So the first real work was rewriting the reader to stream — hand back the file’s header immediately, then yield rows one at a time and never hold more than one. The classifier only needs running totals anyway. It doesn’t need the rows; it needs the counts.

That worked. It processed the nine-and-a-half-million-row file in about three and a half minutes.

Then I did the arithmetic on the other 5,376 files and got roughly seventeen hours.

Where the time actually goes

Seventeen hours is survivable but annoying, and before optimizing anything I wanted to know what I’d actually be optimizing. So I instrumented a single mid-sized file and measured each stage separately.

I expected the answer to be decompression. Gzip is real work, and 59 GB is a lot of bytes.

It wasn’t even close:

Stage	Share of time
Decompressing the gzip and reading lines	6%
Splitting each line into fields	8%
Building a dictionary and an object per row	45%
Accumulating the feature counts	41%

Six percent I/O. Eighty-six percent Python creating and discarding small objects, millions of times over.

The specific culprit was a line I’d have called harmless: for each row, zip the column names together with the values to make a dictionary, then look fields up by name. It reads beautifully. It also means allocating a fresh dictionary nine and a half million times for one file, purely so the next four lines can be readable.

The fix is unglamorous. Work out which column sits at which position once, when you read the header, then index straight into the split fields. Same logic, none of the churn.

I also found a memory version of the same mistake. To detect sites carrying more than two alleles, I was storing a little string like "A|G" for every genomic position I’d seen. There are only a handful of distinct allele pairs in the whole file, but I was storing millions of separate copies of them. Interning each distinct pair to a small integer took about five lines and cut the memory accordingly. (Measured peak on the big file before the change: 1,040 MB.)

Together: about 30% faster and materially lighter.

The part that actually mattered

Here’s the thing, though. These counts aren’t statistics I’m reporting. They decide things. Whether a score has indels, whether it has multiallelic sites, whether any of its variants sit on the X chromosome — those determine whether the score is marked eligible to compute at all. Get a count wrong and a score that should have been deferred gets computed instead, and you’re back to confidently wrong numbers.

So I’d just rewritten, for speed, the exact code path that decides what to trust. That deserves more than “the tests still pass.”

Two things.

First, structural: I made the two code paths physically incapable of disagreeing. There’s now one row-parsing function and one feature-accumulating object, and both the load-it-all-into-memory path and the streaming path drive the same ones. Not “kept in sync” — the same code. The failure mode I was most worried about was a score classifying one way when read one way and another way when read the other, depending on which part of the system happened to touch it. That’s now not expressible.

Second, empirical: I wrote out the original implementation again, verbatim, dictionaries and all, as a reference. Then I ran both versions across sixteen real catalog files spanning 72 variants to 1.1 million, covering classes A, B and C, and compared every feature count, both bounds, and the final classification.

Zero mismatches.

That’s the check I’d have been tempted to skip, because the optimization was “obviously” behavior-preserving. Obviously-correct is how the original bug got in.

Running it

The job is embarrassingly parallel — every file is independent — so I split it six ways and set it going. One detail worth mentioning: I sort each shard’s files largest-first so the six shards get comparable byte loads rather than comparable file counts. File sizes here span four orders of magnitude, so counting files would have handed one poor shard all the giants.

That created a small problem I didn’t anticipate. Because each shard starts with its biggest files, and the code only saved progress every 25 files, the first save was going to be about an hour away — an hour of work that a crash would throw away, and an hour of staring at a progress bar reading zero. I caught it a few minutes in, stopped everything, added a time-based save alongside the count-based one, and restarted. Ten minutes lost, three hours de-risked.

Three hours later, all 5,377 scores were classified. 3.45 billion variant rows examined. About seven in ten came out computable under the current rules — the rest deferred, each with a recorded reason.

The most useful number wasn’t the headline, though. It was this: of the ~1,080 scores held back, 588 are blocked by exactly one thing — sites where a single position carries more than two alleles. Nothing else. Solve that one problem properly and 588 scores move at once. Another category showed up in 399 deferred scores but was almost never the lone blocker, so fixing it would buy nearly nothing.

(⚠ 2026-10-05: superseded. Two parser fixes in September re-ranked the blockers: under classifier v3, rows with no usable position were the sole blocker for 1,437 scores and multiallelic sites for 132, and v4 then admitted the first group at a coverage threshold. Eligibility went from about seven in ten to 43 % at v3 and back to 69.8 % at v4. The manifest was lying has the story.)

I would not have guessed that ranking. I’d have chased the category that appeared most often, which is the one that turns out to be nearly worthless to fix. Counting “what appears most” and counting “what is the only thing in the way” give you different answers, and only the second one tells you where to spend the week.

Two things the run turned up that I want to record honestly:

Five of the cached files were truncated — interrupted downloads from months ago, sitting there looking fine. Every one failed loudly on its checksum rather than being read up to the break and classified on partial data. Re-downloaded, they were 25–62% bigger than what I’d had. We had been holding incomplete files without knowing it.

Fifteen scores parse to zero rows in my reader, because they store their weights in a format it doesn’t read. My first instinct was to flag those as unusable. That’s wrong, and in a way worth naming: zero rows was a fact about my parser, not about the file. The honest label is “this needs a different calculator,” which is both true and actionable. I’d written the check the other way round and had to turn it around.

One more thing, from the wreckage

There’s a loose end from the original audit I could finally close. Among the old, wrong scores, forty of them had raw values above 1,000 — suspiciously large.

There’s a free test for this. An additive score is a sum of dosage × weight where dosage is 0, 1, or 2. So the score has a hard mathematical ceiling and floor: add up twice every positive weight, and twice every negative one. Any result outside that window isn’t a surprising biological finding. It’s arithmetically impossible for that scoring file, and therefore a bug.

The old scores were deleted — they were wholly untrustworthy and recomputable — but the pre-deletion backups were kept, so the check still runs against those. The manifest now carries the bounds for every score, so it’s a join.

I want to be careful about what that test can tell me, though, and I wrote the caveat into the tool itself so nobody later reads too much into a clean result: being inside the bounds proves nothing. The original bug produced scores that were wrong but perfectly attainable. Bounds catch impossible answers, not merely incorrect ones. Whatever comes back is a list of things to look at in the next gate — not a clean bill of health.

All forty came back inside their bounds. Which, per the above, clears exactly nothing.

But the table had a column I hadn’t been looking at, and it was the interesting part. Alongside each old score was how many of its variants the old pipeline had actually managed to use. Score after score: 50.8%. 42.6%. 51.7%. Roughly half, over and over, across scores of wildly different sizes.

That’s the fingerprint of the original bug. When the effect allele happens to be the one in the reference genome and the person carries two copies of it, the old code saw no entry and scored zero — dropping, on average, about half the variants in a score. I’d measured that fraction once, on the single score the audit examined. Here it was again, unprompted, in forty independent scores I hadn’t been looking at.

There’s something clarifying about a diagnosis you made from one example turning up in the aggregate without being asked for. I didn’t set out to re-confirm the bug today. It just fell out of a column I’d added for a different reason.

What I keep relearning

The interesting work today wasn’t the classifier. That was already built. It was the unglamorous middle: finding out that my mental model of where the time went was wrong by an order of magnitude, and then taking seriously that a performance change had landed inside a trust decision.

The measurement took twenty minutes and redirected the whole afternoon. The verification took another twenty and is the only reason I believe the result.

The gap worth naming is between “I made it faster” and “I can show it still says the same thing.” The first is a stopwatch. The second needs a verbatim reimplementation and a corpus to run both against, and it costs more than the optimization did. When the output is a ranking over drugs for one person rather than a benchmark number, that cost stops being optional.