Four wrong explanations for one bad file
2026-09-25, continuing straight out of the previous day. Gate D’s reference panel staged and verified, the locus union rebuilt at 3,753 scores, and one number in it that took four attempts to explain.
Where this picked up
Classifier v4 had taken the eligible set from 2,322 to 3,753 scores, so the locus union and dense callset had to be rebuilt at that scope before Gate D’s panel work could key to them. That rebuild is the spine of this entry, and almost everything interesting happened around it rather than in it.
Two things landed cleanly and are worth stating before the messy part, because they are the actual deliverables:
The panel is staged and verified three ways. 14.86 GB of pgsc_HGDP+1kGP_v1
over 8 ranged connections (~4.6 h; EBI throttles per connection at ~180 KB/s, so a
single stream would have been ~30 h), md5 matching the catalog’s published
md5s.txt, then every extracted GRCh38 file checked against the checksums.txt
shipped inside the archive. 3,942 samples with SuperPop and Population
labels, plus a 613-sample relatedness exclusion list — exactly what the PCA needs.
I reclaimed 14 GB by deleting the GRCh37 half the archive also ships, which a
GRCh38-pinned pipeline will never touch.
The orientation precondition holds. The panel’s REF alleles agree with our pinned FASTA on 20,000 of 20,000 sampled loci. That is the one check Gate D cannot proceed without: Gate C resolves orientation from the reference allele alone, so if the panel disagreed with the FASTA at a locus, the patient and the panel would be scored under different orientations there and the percentile would be comparing two different things.
The union rebuilt and reconciles. 19,215,298 loci from 3,753 contributing scores over 2,166,261,885 rows, and the four callability outcomes sum to the union size exactly: 11,822,237 + 4,192,781 + 949,225 + 2,251,055. 7.2 hours of wall time, on the third attempt.
The thing I got most wrong
Before any of the science: I told Richard that merging a pull request was safe and
only a git pull into C:\Maple was dangerous, because a pull changes the .py
files the workers watch and WORKER_AUTO_RELOAD=1 turns that into a restart.
That was wrong, and it cost about four hours twice over.
(Correction, 2026-09-26, GH #361: this section’s conclusion is itself wrong. Nothing updates the local checkout when a pull request merges; the fast-forwards below were other sessions syncing it by hand, and my first statement was right. The reloader has been off since 2026-09-28.)
The reflog on C:\Maple’s main:
a500dff main@{2026-09-24 19:47:36 -0700}: merge a500dff82...: Fast-forward
d672467 main@{2026-09-24 16:14:19 -0700}: merge d6724677...: Fast-forward
574c5de main@{2026-09-24 10:38:47 -0700}: pull --ff-only: Fast-forward
And the worker’s own log, which is the part I should have read hours earlier:
19:20:27 claims task 1041504 <- attempt 1
23:16:36 [reloader] 1 file(s) changed <- kills it at 50.76%
23:21:52 claims task 1041504 again <- attempt 2
02:50:21 [reloader] 10 file(s) changed <- kills it at 49.28%
02:55:15 worker10 claims it <- attempt 3, the one that finished
Something fast-forwards C:\Maple to origin/main within about three minutes of a
merge. So a merge is a worker restart, transitively. Worse, d672467 is not one
of my commits — another session merging to main is enough to do it. A 3.5-hour job
that any unrelated merge can reset, with no checkpoint, is structurally fragile,
and the spill phase is the part that gets discarded.
There is a smaller error nested inside this one. When attempt 1 died I had just run
a test suite that loads a 4.55 GB structure, so I told Richard I had starved the
worker. The log says [reloader] 1 file(s) changed — a memory-starved process does
not print that. I reasoned from coincidence in time and announced a cause I had not
checked, in a session where the evidence was one docker logs away. The pytest run
was genuinely reckless on a loaded box; it just was not what killed the job.
The 12% anomaly, and four explanations for it
structural_fail went from 2.10 % of the union to 11.71 % — 2,251,055 loci that
cannot be scored. And the reason strings were self-contradicting:
1:833068 FASTA REF (G) not among scoring alleles ('G', 'G')
G is plainly among ('G','G'). That message is false as written, and it is the
thread that led me wrong twice.
Explanation 1: memory pressure. Already dealt with above. Wrong; it was the reloader.
Explanation 2: the union builder is corrupting allele pairs. This one felt
strong. The source file says 1:833068 effect='A' other='G', the classifier
returns ('A','G'), and the shard on disk says ('G','G'). Something between them
had to be mangling it — the spill round-trip, the intern table, the shard writer.
The test that settled it is one I like: build a union from a single score. One
score cannot collide with itself, so any cross-row artifact has to be the builder’s.
Result: n_targets=8,598,112, multiallelic_sites=0, 0.00 % degenerate, every
target position emitted exactly as declared. The builder round-trips correctly.
Hypothesis dead, and dead cheaply — minutes, not a re-read of 59 GB.
Explanation 3: it is concentrated in a few scores, findable by sampling. I screened 25 eligible scores through the real classifier. Zero degenerate pairs. Then I did it again on 1,470 genome-wide scores and gave up on wall-clock.
This was the same mistake I had written a memory about the previous morning — that sampling for a rate and attributing items to sources are different questions needing different scans. I wrote it down, and then sampled anyway.
Explanation 4, the right one: one published file is entirely malformed.
Screening every contributing score on its first few thousand rows rather than its whole content turned an overnight job into minutes, and a 100 %-density contributor cannot hide in that window:
PGS005228: rows=2,379,429 placed=2,379,397 degenerate=2,379,397 (100.0000%)
PGS004691: rows=950,353 placed=950,342 degenerate=16,657 (1.7527%)
PGS005228 is “FI_Frailty”, a Frailty Index score from Foote IF et al., Nature
Genetics 2025. Every single placeable row declares the same base in both allele
columns:
rs7909677 10 111955 G G 2e-06 hm_inferOtherAllele=(empty)
rs4881491 10 117636 G G -2.8e-05 hm_inferOtherAllele=(empty)
2,396,054 degenerate rows contributed against 2,251,055 structural_fail — the
accounting closes, with the ~145k gap being degenerate loci deduped against each
other plus those at variant sites, which behave differently (below).
Why a degenerate pair is fatal, and why it is invisible
I coined “degenerate allele pair” during this and should have flagged it as my
term rather than a standard one. It means the two declared alleles are the same
base — ('G','G') — degenerate in the mathematical sense: a structure meant to
distinguish two things has collapsed to one value.
Dosage asks how many copies of the effect allele does this person carry, 0/1/2.
That only means something against an alternative. With both alleles G there is no
contrast, no variant described, nothing to count — and orientation is unanswerable
in a peculiar way, because the FASTA’s REF=G matches both declared alleles, so
effect cannot be told from baseline. That is the real structural failure, and it is
exactly what the error message garbles into “not among”.
The contrast that makes it click:
| alleles | recoverable? | |
|---|---|---|
palindromic (A/T, C/G) | two different alleles that are complements | yes, from a reference panel’s allele frequencies — that is a large part of what Gate D is for |
degenerate (G/G) | one allele, twice | never, by any amount of reference data |
And Gate B cannot see it. 'G' is a perfectly valid allele sequence, so the
n_ambiguous_other check (which catches A/T candidate sets, N, -, star
alleles) does not fire. Nothing tests whether the two alleles are the same. So a
2.38M-row file that is 100 % unscoreable classified as production-eligible, and
we found out as 2.25M structural failures at the end of a seven-hour union build.
That is the finding worth more than the bad file: the classifier’s allele checks
ask “is each allele a resolvable sequence?” and never “do these two alleles
describe a contrast?”. Filed as GH #334, with the fix being an n_degenerate_pair
feature that is non-promoted, and Class D for a file that is entirely
degenerate — which is precisely what “no usable variant rows” already means.
The downstream effect, stated more carefully than I first stated it
My first instinct was that 2.25M loci of the callset were wrong and E2 was blocked. That was an overstatement and I corrected it within the session, but it is worth recording what the artifact actually looks like, because the two site types differ:
- hom-ref sites — the degenerate row structural-fails, because there is no non-reference allele to pair with the FASTA REF. One row emitted. Honest.
- variant sites — it is accepted, because the builder takes REF/ALT from the VCF record rather than from the scoring pair, and emits a duplicate identical row:
variant site A <- emitted twice
variant site B <- emitted twice
hom-ref site C <- emitted once; the degenerate path failed here
DenseCallset.get() returns a list of candidates and _select_call picks among
them, so two identical candidates should yield the same dosage. Scoring is very
likely correct on an inflated artifact rather than wrong — which is a materially
weaker claim than the one I made first, and it needs confirming rather than
assuming before E2.
What is definitely true: the union is ~12 % padding, multiallelic_sites (39,807 →
2,432,297) is largely fiction because a degenerate pair forms its own allele set
({'G'} versus {'A','G'}) and so reads as multiallelic wherever another score
covers the same locus, and one eligible score will score exactly nothing.
Two other things this session cost, both mine
A retraction. Everything I reported about panel coverage on the 24th was
measured against the wrong file. I used pgsc_all+hm3_variants_hg38_latest.txt.gz
on the strength of its name sitting in the panel’s own FTP directory, plus a size
estimate that seemed to corroborate it — 16 GB over 4,150 samples at a quarter-byte
per genotype is about 15M variants, near enough to 33.8M to feel confirmed. Two weak
signals agreeing is not verification; it is one signal. The panel’s actual .pvar
holds 84,325,213 variants. Real union coverage is 98.93 %, not 94.77 %, and
PGS005353 is 99.65 % covered, not 77 %. The whole nine-score “accession band”
finding and the per-score-coverage-flag design requirement I derived from it were
artifacts of my file choice, and are retracted in place in prs_primer.md §5.7
rather than quietly edited out. What caught it was inspect_reference_panel.py,
which reads the actual artifact, on its first run, in one line of output.
A wrong bug report. I filed GH #315 claiming the union handler reports
progress_pct by file count. It does not. The formula is a 3 %→9 % band over rows
against a hardcoded 894_528_461 — the v3 row count, now stale at 2,166,261,885 —
which reproduces the observed 3.88 % exactly and means progress freezes at 9 %
once rows pass 41 % of the real union. An indicator that stops moving on a healthy
job is worse than a wrong ETA. Corrected the issue and retitled it, but I had
filed a plausible inference without checking it against the code, in an issue whose
entire subject was a misleading metric.
What actually shipped
Two PRs merged: #329 (classifier v4 — unplaceable_rows as a ≥99 %-placeable
threshold, eligible 2,322 → 3,753, 894.5M → 2.166B rows, nothing lost) and #331
(artifact identity — the union reuse guard now compares contributing_pgs against
the requested scores instead of testing that a manifest merely exists, which was a
live path to silently scoring against a callset missing 1,431 scores’ loci).
Three issues filed: #330 (artifact paths encode who and how, never over what scope, so a 6-locus sentinel and an 11.4M-locus production build claim the same filename — with an acceptance test that rebuilds the destroyed sentinel and compares checksums, and explicit instructions for what a mismatch means), #334 (this one), and #315 corrected.
The union is rebuilt, the panel is staged and REF-verified, and Gate D’s actual
build — panel scoring, PCs, the mean/variance regression — is next. It needs a
rebuilt worker image for pgenlib, which currently exists only as an ad-hoc
install inside a running container.
The lesson I would keep if I could keep one
Four explanations for one anomaly, and the order I tried them in was almost exactly backwards. Memory pressure, then our builder, then sampling harder — and only last the cheapest hypothesis available: that the input was bad. The same inversion showed up in the coverage retraction, where I spent hours explaining a discrepancy from a reference I had never verified.
An unexplained gap is at least as likely to be an instrument artifact or bad input as it is to be a finding. “Check what you are measuring, and what you are measuring against, before you explain the measurement” would have saved most of two days.