The issue arrived with a headline that was hard not to like. Gate C defers every palindromic
SNV (A/T and C/G pairs, whose strand cannot be read off the letters), and #371 had just
closed the hm_source route, leaving declared allele frequencies as the only instrument:
8.48 % of the deferred rows. #404 pointed out that the same spike had proved the harmonizer
never touches an allele, so a file’s unambiguous rows are a readout of its author’s strand
convention, and 99.899 % of the deferred palindromic rows sit in files whose unambiguous
rows match the reference forward. Roughly twelve times the coverage, if the palindromes
followed suit.
The issue also named the one way it could be wrong: a pipeline that orients what it can check and leaves the palindromes alone. The 167 scores that declare an effect-allele frequency are the only place that failure could be caught. This session was the measurement.
Building an instrument that could fail
The test itself is small. For an A/T pair, reading the other strand turns the effect allele
into the other allele, so orientation reduces to asking whether the author’s declared
frequency and the panel’s frequency sit on the same side of 0.5. The existing
resolve_palindrome_strand already encodes it. The work was everything around it, because a
frequency comparison is noisy, and raw agreement below 100 % would have meant nothing.
Three things went into making the number mean something.
- Calibration on rows whose answer is known. Each file’s unambiguous SNVs have a strand
fixed by the panel’s own alleles. So the instrument’s error rate is measured in the same
file, cell by cell in distance-from-0.5, and a palindrome’s discordance is compared against
that, not against zero. On synthetic flip rates of 0, 1, 5 and 50 % the estimator recovered
each within its interval. And since observed discordance is
e + f(1 − 2e), the raw rate bounds the flipped share for any error rate below one half, a bound that needs no calibration at all. - A check that the frequency column measures anything. Eight of the 167 files fail on their own known-strand rows. One, PGS001775, calls 0 % of them right: its column is the other allele’s frequency. Pooled uncritically, it would have read as a file whose every palindrome was flipped.
- A check against circularity. A frequency looked up from a forward-strand reference by allele letter agrees with the letters whatever strand they were written on. It is the REF-effect lesson again: a check that shares the assumption it is testing cannot fail. 1000 Genomes samples sit inside the panel, so a copied frequency lands almost exactly on those subsets. One publication did (a median 0.0015 from 1000 Genomes EAS, against 0.0055 or more for every independent cohort), and was set aside.
Reading the panel took a detour I should have seen coming. The container needed pgenlib,
which a live worker had, so I launched a throwaway container from the same image. The live
workers hot-reload on any .py change under C:\Maple, and would kill a long read
mid-flight. The throwaway container did not have pgenlib: it had been pip-installed into
one worker’s writable layer, after the image was built, and exists in exactly one of our ten
workers. The failure was loud in the wrong way. The import happened inside a
multiprocessing.Pool initializer, which fails in every child, and Pool respawns them
forever. It produced enough identical tracebacks that the notification monitor was shut
down for flooding. The fix
was to copy the single .so out onto PYTHONPATH, and to import it in the parent first so
the failure happens once. Scattered reads through the Windows bind mount run at about
1,000 variants/s; twelve parallel readers over sorted chunks did 8,041,370 records in
7.3 minutes.
A synthetic smoke test before the real run caught two bugs I would otherwise have met after
it. Python’s sum() of an empty tier returns the integer 0, which has no .sum(). And my
pooled estimates were still including the files the validity check had just rejected.
What it measured
| Class of file | Palindromes reversed (95 % bound) |
|---|---|
| Genome-wide, 11 files from 6 publications | ≤ 0.12 %, rising with each file’s own complement-only rate |
| Small curated, zero complement-only rows, 61 files from 49 publications | ≤ 0.96 % |
| PGS000020, the one mid-size file | ~40 % |
| An A/B-coded control | 5 of 5 |
The genome-wide files were the most satisfying part. Their flip rates track their own complement-only rates almost one for one, from zero with zero to 0.098 % with 0.141 %. That is what uniform residual strand error looks like: whatever put a tenth of a percent of the unambiguous rows on the wrong strand put about the same share of the palindromes there.
PGS000020 is the failure the issue predicted, found in the flesh. It is a type 2 diabetes score of 7,502 rows. Its unambiguous rows match forward (12 complement-only of 6,867, 0.175 %) and 169 of its 425 decidable palindromes are reversed, across both pair types. Its complement-only rate is within a factor of two of the cleanest genome-wide file’s, 0.098 %. A threshold at 0.1 % would have separated them. I did not take it: it would have been drawn around the one counterexample after seeing it, and the rate demonstrably cannot see this failure.
The turn
The result I had not looked for came from asking where the rule would land. The genome-wide class holds 98.8 % of every deferred palindromic row, which is where the issue’s headline lives. 87 % of those rows are in three publications that declare no frequency at all: a multi-method benchmark, a PRSmix study, and a metaGRS study. The two mixture methods are weighted combinations of existing Catalog scores, so their palindromes carry whatever orientation their components had. The six publications the instrument could measure hold 9.7 % of the class.
So the recommendation came out much narrower than the headline: license the small curated class (527 scores, 14,250 rows, but a median 11.3 % of each of those scores’ rows), put the genome-wide class to Richard as an explicit extrapolation, and never offer the mid-size class. The mechanism landed inert, the same shape as the per-file strand-flip recovery. Its candidates script offers the genome-wide class only behind a flag that prints the caveat, and a test keeps the mid-size class out. #414 is the route to measuring the genome-wide class properly, through the source GWAS summary statistics.
What I got wrong
Three numbers in the write-up’s first draft, all caught by rechecking before it was committed:
- I wrote “~1,500” decidable curated palindromes; it was 1,831.
- I said sampling noise kept independent cohorts “above ~0.008” from the panel. One sat at 0.0055, which also meant the code comment beside the threshold was wrong.
- I claimed the genome-wide rule would recover “~15 % of loci” in those scores. It is 10 %.
I also misjudged one stage’s elapsed time by a factor of three before checking it, a small error of exactly the kind that turns a patient wait into a pointless intervention.
Two things surfaced that were not errors of mine but were nearly mistaken for regressions.
update_memoryrejected a call for want of aconfidenceargument, although that bug had closed thirteen hours earlier. The MCP server process predated the fix by four.- The policy fingerprint had already drifted from E2’s before I touched it, because #376 added a knob. My change moved it again. It is correct by design, and it means Gate D will refuse every E2 row until E3 re-scores them. I wrote that down where the next session will find it, not where it would be discovered mid-run.
The merge brought the day’s other surprise. The morning’s #405 had written the hypothesis into the same two documents where I was writing its measurement. The conflict resolved itself in the obvious way: keep the hypothesis as it was written, and put the answer directly under it.
Richard’s last question was whether a palindrome, being the same both ways, is harmless to reverse. It is the right question, and its answer is the whole problem. The pair is the same from both strands, which is why the flip is invisible. The effect allele is one letter of the pair, and a flip hands the weight to the other one. So a homozygote’s dosage goes from 2 to 0, and only heterozygotes come through unchanged.