Polygenic Risk Scores — a primer, and how Maple computes them

What a polygenic score is, where the idea came from, what a scoring file actually contains, and the machinery MODS uses to turn one genome plus five thousand published scores into numbers we are willing to stand behind.

Status as of 2026-09-26: Gates A, B, C, E1 and E2 are done; D is part-staged and F is not started. (2026-09-29: stale on D and F. Both have now run for the v0 set, and the Gate F v0 release writer is built: backend/mods_prs_release.py, driven by backend/scripts/prs/release_gate_f.py (GH #296); see §5.5.) (2026-09-27: all three Gate D stages are now built — D3 in PR #386, D1 and D2 in GH #396 — and none has run against the real panel; §5.7.) The scorer has now been run over the whole v4 eligible set: 3,752 of 3,753 scores computed against a 16,964,243-record callset (19,215,298-locus union), every one reconciling and every one inside its theoretical bounds, with the fleet’s outcome counts summing to the manifest’s eligible row total exactly (§5.8). No released patient scores exist — Gate C writes at validity='quarantined' and only Gate F promotes, so every reader still correctly returns nothing (§5.5). (2026-09-28: superseded by the first release of her PRS: a v0 smoke release of 590 scores over 109 disease anchors (tag smoke_v0, release run c3be9df0), normalized against the 667-sample EUR reference, after D1–D3 ran on the 594 v0 candidates. The final Gate F is still ahead; §5.5.)

E1 was a measurement run by design, and it paid for itself the same day. Its histograms closed the largest open policy question (unplaceable_rows, §3.5) and corrected this document’s palindrome figures by 4× (§5.4), so classifier v4 promoted unplaceable_rows at a coverage threshold and the eligible set went 2,322 → 3,753 scores, 894.5M → 2.166B variant rows. Every “2,322” below is a statement about the E1 run, which is history and correct as written; the forward-looking number is 3,753.

Read “E2” carefully: the running order changed, and the name did not. The plan below was E1 → settle policy → D → E2, where E2 meant “the re-run under the final policy — the one that can be released.” That is not what E2 was. Gate D has never run (mods_prs_run holds zero gate_d rows), so E2 is better understood as E1 repeated at the corrected v4 scope — a second measurement run, not the post-D final one.

This was not drift for its own sake. E1’s histograms settled unplaceable_rows from E1’s own data, needing nothing from the panel, and that alone moved eligibility 2,322 → 3,753 — which invalidated the union and the callset. Since D1 extracts frequencies over the union, the union had to be rebuilt at final scope before D or it would be built twice. Rebuilding it was therefore correct and necessary; re-scoring the patient against it was the optional part, and it bought the measurements that close out the remaining policy questions (§5.8) plus a clean baseline to diff against later.

So E2 is final with respect to eligibility, and not final with respect to palindromic SNVs, which Gate C still defers — 64,810,321 rows, 2.99 % of the fleet.

Decided (Richard, 2026-09-26): no release candidate is declared until the palindrome question is settled, and an E3 follows D1. So E2 is explicitly not the release candidate. Note what “settled” can mean, because it is narrower than it sounds: only 167 of the 1,687 affected scores carry the effect-allele frequency that orientation needs, so D1 can recover at most 9.5 % of the deferred rows and the rest become a Gate F coverage gate rather than a pending decision. The staged plan, with what each stage is for, is §5.7b. Nothing else is outstanding except the NR weight-type call, which is Gate F’s and does not require re-scoring.

Companion documents. `genomics_glossary.md` defines the vocabulary for readers who want the terms first — allele, dosage, strand flip, gVCF, plus what “calling” means (§2.1, with GQ/DP/MIN_DP/FILTER), contig, and FASTA vs FASTQ including the provenance of our pinned reference (§4.1). `prs_validation/compatibility_classification.md` is the standalone explainer for the classifier in §4.5–4.6. This primer is the wider arc that connects them.

1. What a polygenic score is

1.1 The arithmetic

A polygenic score is a weighted sum, and that is genuinely all it is:

PRS = Σ  dosageᵢ × weightᵢ

For each of a set of positions in the genome, dosage is how many copies of the trait-associated allele this person carries — 0, 1, or 2, because we have two copies of each autosome — and weight is how much that allele shifted the trait in a large association study. Multiply, add, done. A thousand multiplications and one sum.

The entire difficulty of the field lives in three places, none of which is the arithmetic:

Getting the weights. That is the GWAS literature’s job, and it is upstream of us.
Lining the two sides up. Deciding that this row in a scoring file and this position in this person’s genome refer to the same thing, on the same strand, with the same allele playing the same role. This is where essentially every bug we have hit actually lived.
Interpreting the total. A raw sum is a number with no units and no meaning until it is placed against a reference distribution of people (§1.4).

1.2 A short history

When	What happened	Why it matters here
1918	R. A. Fisher shows that many small Mendelian effects, added up, produce the smooth bell-curve inheritance seen in traits like height — the infinitesimal model.	This is the theoretical licence for adding thousands of tiny effects together. The additive sum in §1.1 is Fisher’s model, unchanged.
1965	Falconer’s liability-threshold model: a yes/no disease is modelled as a continuous underlying liability with a cutoff.	Why it is coherent to give someone a continuous “risk score” for a binary outcome like colorectal cancer.
2002–2005	The HapMap project maps linkage disequilibrium — which variants travel together — making cheap tag-SNP genotyping arrays viable.	Made large-scale association studies economically possible. Also the reason scores can use proxy variants rather than causal ones.
2005	Klein et al. link complement factor H to age-related macular degeneration: the first widely accepted GWAS hit.	Proof the approach finds real biology.
2007	The Wellcome Trust Case Control Consortium publishes 7 diseases × ~14,000 cases.	The GWAS era starts in earnest; effect sizes turn out to be tiny, which kills the one-gene-one-disease hope and makes aggregation mandatory.
2009	The International Schizophrenia Consortium (Purcell et al.) takes thousands of SNPs that individually failed significance, weights them by effect size, sums them, and shows the aggregate predicts schizophrenia in an independent cohort.	This is the birth of the modern polygenic score. The key insight: signal too weak to be individually credible is still real in aggregate.
2015–2019	Method generation: LDpred, LDpred2, lassosum, PRS-CS, SBayesR, MegaPRS. Instead of picking significant SNPs, these re-weight every variant using an LD reference panel.	This is why the modern catalog contains scores with millions of variants and no rsIDs (§3.3) — they are genome-wide posterior weights, not curated hit lists.
2018	Khera et al. publish genome-wide scores for five common diseases and argue that the top few percent carry risk comparable to a monogenic mutation.	The paper that pushed PRS toward clinical translation — and toward controversy.
2017–2019	Martin et al. quantify the portability problem: scores trained in European-ancestry cohorts lose much of their accuracy in other ancestries.	A first-order limitation, not a footnote. It determines how we normalize (§1.4) and is why ancestry is a pipeline input, not a nicety.
2021	The PGS Catalog launches (Lambert et al.) with a standardized scoring-file format, harmonized coordinates, and a reference calculator (`pgsc_calc`).	The entire supply side of our pipeline. Without a common file format, §3 does not exist.
2022	T2T-CHM13: the first gapless human genome.	Not yet our reference (§2.2) — but it is where the field is going.
now	Prospective clinical trials (eMERGE, GenoVA, WISDOM), regulatory caution, and a large direct-to-consumer market of varying rigour.	The reason a pipeline like ours has to be able to say “not calculated, and here is why” out loud.

The through-line: the idea (Fisher) is a century old, the data (GWAS) is twenty years old, the practical method (ISC) is fifteen, and the standardized plumbing (PGS Catalog) is five. We are working entirely in that last layer.

1.3 What a polygenic score is not

Worth stating plainly, because the failure modes of this pipeline are mostly failures of interpretation rather than of code:

Not a diagnosis, and not a probability. A raw PRS is an unnormalized sum. On its own it is meaningless — it only acquires meaning relative to a population.
Not portable across ancestries without care. See 2017–2019 above.
Not the whole genetic story. A PRS captures common, additive liability. It says nothing about a rare high-penetrance variant. A rare high-penetrance germline variant is exactly the kind of thing a PRS structurally cannot see — which is why the variant-level pipeline and the PRS pipeline are separate inputs to the same model.
Not causal. Weights sit on tag variants in LD with whatever is actually doing the work.
Not comparable between scores. Two scores for the same trait, from different studies with different weight conventions, produce numbers on different scales.

1.4 Why Maple computes them

This matters for design, because our use is not the usual one. MODS does not compute PRS in order to hand somebody a clinical risk report. It computes them as a broad, patient-specific genetic feature vector that feeds the PDR2 personalized loss:

L = L_pred − λ₁·L_prs − λ₂·L_expr

where L_prs aggregates, over diseases i, a learnable weight wᵢ times the patient’s PRSᵢ times the model’s per-triple score. The polygenic score for each trait becomes a signed, learnable nudge on drug–disease edges incident to that disease node — the patient enters the objective, not just the inputs.

Take the shape, not the algebra. The paper’s literal formula (L_prs = Σᵢ wᵢ·PRSᵢ·s_pred) is unbounded as written, and the released code adds tanh(·) and L2 weight normalization that the text omits. We verified those nonlinearities are load-bearing. losses.py is the operational specification — see `pdr2_loss_text_vs_code_note.md` and `pdr_personalization_writeup.md`.

Three consequences follow, and they explain choices that would otherwise look strange:

Breadth beats depth. We want as many scores as the genome supports, not a curated clinical shortlist. A wide feature vector is the point. This is why the target is the ~5,377-score catalog and not thirteen diseases.

The trap this is guarding against, stated bluntly (Richard, 2026-09-23, after a session report led with “44 colorectal scores are eligible”): counting disease-relevant scores is not a result. We already know she has colorectal cancer, and we know a great deal about it. A colorectal PRS can, at best, tell us whether she was unlucky in getting it — which changes no decision. No individual score and no individual disease is the object of interest. These are ML inputs. The goal is to maximise the number of usable data points without regard to what they predict; anything that refines that data is a direction to move toward, and anything that degrades it a direction to move away from. When tempted to report “N scores for disease X”, report the size and quality of the feature vector instead.
The relative ordering across a cohort matters more than any absolute number — but we have a cohort of one, so a reference distribution is required to make the value mean anything at all. That is Gate D.
A wrong number is worse than a missing number. A missing feature is a zero the loss ignores. A wrong feature is a confident, learnable lie. Every “we defer this rather than approximate it” decision below traces back to this line.

A corollary that changes how the partial-coverage decision should be made — and a correction to how this document used to justify it. A score computed on only part of its loci is not simply “wrong” the way it would be clinically, so the coverage bar for “usable as an ML feature” sits genuinely lower than for “usable as a clinical risk estimate.” That conclusion stands. The reason previously given for it did not.

Superseded (2026-09-23). This section used to argue that a roughly constant deflation is “largely absorbable by wᵢ — the model can rescale a feature.” The released loss cannot do that. losses.py L2-normalizes the weight vector before the inner product and then squashes the result:
disease_weights_norm = F.normalize(self.disease_weights, p=2, dim=0)
prs_impact = torch.sum(disease_weights_norm * prs_values)
prs_term   = torch.tanh(prs_impact)
‖ŵ‖₂ = 1 is a hard constraint, so the model cannot rescale a feature — it can only re-allocate a fixed unit budget across diseases. Feature scale therefore is not absorbed; it decides the answer.

What actually absorbs a constant deflation is Gate D, not the model. A percentile (or a z-score) is computed against a reference distribution, and a multiplicative deflation applied equally to the patient and the panel is a monotone transform — it preserves rank exactly, so the normalised feature is unchanged. The absorption is in the normalisation, which is another reason Gate D is not optional (§5.5).

Three things keep it from being free, and the first is now the mechanism rather than a caveat:

The deflation must be consistent between the patient and the reference distribution (Gate D). Coverage that differs between Maureen and the panel biases the percentile, and nothing downstream can undo a bias introduced in normalisation. Score the panel through the same Gate C policy or this guarantee evaporates.
It is not a pure scale: which loci are missing matters, and a few large-weight loci can deflate a score disproportionately, which breaks the exact rank preservation above. Over many loci it approximates one.
Fewer loci means a noisier estimate of the same liability. That is variance, not bias — and variance is far more tolerable in a feature than bias is.

The corresponding trap on the other side: raw sums are not a usable feature vector. Across the 2,322 production-eligible scores the theoretical magnitude bound max(|lo|, |hi|) spans roughly twelve orders of magnitude — median 15.9, p95 706, max 2.59 × 10⁷, with 1,818 of 2,322 (78 %) already above 3. Under tanh(⟨ŵ, PRS⟩) with ‖ŵ‖₂ = 1, a vector on those scales saturates the squash, 1 − tanh² ≈ 0, and the gradient dies. Track A (all-zero PRS) produced no gradient because tanh(0) = 0; feeding raw sums would produce no gradient from the opposite end, and the “interpretable” disease weights would converge to a readout of which published score happens to carry the largest raw magnitude rather than which disease matters. Normalised features are what put ⟨ŵ, PRS⟩ back inside tanh’s live band.

The operational consequence: unplaceable_rows should not be a binary promote/defer flag. It should be a coverage threshold, chosen empirically by what it does to the feature vector. The instrument for choosing it already exists — the coverage_v1 payload’s third quantity (effective variance-weighted information, §5.4) is exactly the right statistic, and the dual-scorer reconciliation (§5.4) separates the imputation term, which measures precisely how much a partially-covered score changes when the missing dosages are filled in. Run that on the affected scores and the threshold stops being a judgement call.

2. Reference genomes

Every coordinate in this pipeline — in the patient’s gVCF, in a scoring file, in dbSNP — is meaningless except relative to a specific reference genome build. Mixing builds is the single easiest way to produce a plausible, entirely wrong score.

2.1 What the reference actually is

The human reference genome is not any one person’s genome. It is a curated haploid consensus assembled from a handful of anonymous donors (with one donor, RP11, contributing a majority), maintained by the Genome Reference Consortium. It carries the reference allele at every position by convention, not by biological privilege — at plenty of loci the “reference” base is the minor allele in most populations.

This matters directly for us: §5.1’s bug is precisely a case where being reference-matching is not the same as being effect-free.

2.2 The builds

Build	Released	Notes
NCBI36 / hg18	2006	Still the original build of 11 catalog scoring files.
GRCh37 / hg19	Feb 2009	The workhorse of the GWAS era. 3,951 of our 5,377 scoring files were authored on it (plus 415 more labelled `hg19`).
GRCh38 / hg38	Dec 2013	Current production standard. Fixed thousands of assembly errors, modelled centromeres, added alt loci for highly polymorphic regions (HLA).
GRCh38.p1 … p14	2014–2022	Patch releases. Patches add fix/novel scaffolds but never move an existing coordinate — so p13 and p14 are coordinate-compatible on the primary assembly.
T2T-CHM13	2022	The first truly complete, gapless genome — centromeres, acrocentric short arms, ~200 Mb of new sequence. Not yet our reference: dbSNP, the GWAS corpus, and the PGS Catalog are all still keyed to GRCh38.

GRCh37 and hg19 are the same assembly under UCSC vs GRC naming — with the irritating exception of the mitochondrial sequence, which differs. Likewise GRCh38 and hg38. The catalog uses all four spellings, which is why our build handling normalizes rather than string-matches.

2.3 Naming and assembly flavors

Two orthogonal traps:

Contig naming. UCSC writes chr1, chrX, chrM. Ensembl and the GRC write 1, X, MT. NCBI writes accessions (NC_000001.11). These name identical sequence. Every entry point in our code normalizes through connectors/reference_fasta.py::normalize_chrom (strip a leading chr), and we standardize on the Ensembl convention because that is what the patient’s gVCF and the PGS Catalog’s hm_chr column both use.

Assembly flavor. “GRCh38” can mean the primary assembly (chromosomes + unplaced scaffolds), the full analysis set (primary + alt contigs + decoys), or a no-alt subset. Alt contigs are alternate versions of polymorphic regions; if reads can align to both an alt and its primary counterpart, mapping quality collapses and variant calls quietly degrade. We use the primary assembly — 194 contigs — which matches how the gVCF was called.

2.4 Liftover, and why we do not do it at score time

Liftover converts coordinates between builds using a chain file of mapped blocks. It is lossy: some regions have no equivalent, some invert, some map many-to-one. The GRCh37→GRCh38 shift is large and irregular — on PGS000006, only 3 of 313 loci kept the same coordinate, with shifts up to 2.87 Mb.

Our position is that liftover is the PGS Catalog’s job, not ours. We download only the _hmPOS_GRCh38 harmonized files, so every scoring file arrives already on our build. The patient’s gVCF was called on GRCh38. No conversion happens anywhere in the scoring path, which removes an entire class of error.

2.5 Our pinned reference

Because the FASTA is the arbiter of what the reference allele is at every scored locus, it is pinned in code, by checksum, and a mismatch hard-fails rather than falling back.


Build	GRCh38.p13
Source	Ensembl release-109, DNA primary assembly
Path	`backend/data/reference/GRCh38/Homo_sapiens.GRCh38.dna.primary_assembly.fa` (host-only; `backend/data/` is gitignored)
sha256	`78777b0886e8dfa5e14e4957fbbaa53736fcbaa5668d59e09b6b7945fca93d8c`
Contigs	194
Naming	Ensembl (`1`, not `chr1`)
Pinned at	`connectors/reference_fasta.py::PINNED_GRCH38`
Access	`ReferenceFasta.ref_base` / `.ref_bases`, via pyfaidx

Validated base-identical to the patient’s gVCF REF field at three ancestry/clinically-informative loci: one clinically-informative locus, 5:33951588 = C (SLC45A2), 2:135851076 = G (LCT).

Two operational notes. Ensembl ships plain gzip, not bgzip, so the file must be decompressed to .fa before pyfaidx can index it — scripts/prs/prepare_reference_fasta.py downloads, verifies, and decompresses in one step. And the checksum is not verified on every load (hashing 3.15 GB is slow); it is verified once at pin time and recorded on the mods_prs_run row so any run is reproducible against the exact bytes used.

2.6 The two build fields inside a scoring file — a live gotcha

A harmonized scoring file carries two different build declarations, and they mean different things:

#genome_build= — the author’s original build. This is what pgscatalog.org displays as “Original Genome Build.”
#HmPOS_build= — the build the harmonized positions are on.

Catalog-wide across our 5,377 cached files (scanned from the headers, 2026-07-26):

`genome_build` — the author’s original	Files
GRCh37	3,951
GRCh38	457
hg19	415
NR (not reported)	357
hg38	182
NCBI36 (11) / hg18 (3) / NCBI35 (1)	15

HmPOS_build, by contrast, is GRCh38 for all 5,377 — because we only ever download the GRCh38 harmonized file. Per-row, the same split repeats: chr_position is on the original build, hm_pos is on GRCh38.

A practical consequence worth knowing: open_scoring_file stores HmPOS_build or genome_build, so mods_prs_score_manifest.genome_build is uniformly GRCh38 and the dashboard’s “Genome build” panel is a degenerate one-row table. That column is not the number pgscatalog.org shows; the table above is.

The gotcha. connectors/pgs_catalog_connector.py::_row_to_variant coalesces hm_pos → chr_position and hm_chr → chr_name as a fallback when the harmonized field is blank. For hm_rsID → rsID that fallback is safe, because rsIDs are build-independent. For position it silently substitutes an original-build coordinate into a pipeline that assumes GRCh38. Measured 2026-07-26: it fires in 85 of 200 randomly sampled files, at ~0.0–0.1 % of rows each; the worst case is PGS000016 (atrial fibrillation, authored on hg19) with 1,316 of 6,730,541 rows. Those rows are counted as placed_loci, not as unplaceable, so n_no_position under-reports.

Closed (GH #250), before Gate C computed anything. _ColumnIndex now decides once per file whether the original coordinates share the harmonized frame (original_coords_share_harmonized_frame: true only when the file is not harmonized at all, or when both builds resolve to the same assembly — “NR” does not count). A blank hm_pos on a score authored against another assembly is now an unplaceable row, which is what it always meant. Gate C counts those rows as unplaceable_row and they never enter a sum.

3. Scoring files

3.1 Where they come from

The PGS Catalog (EMBL-EBI + University of Cambridge) is the standardized public repository of published polygenic scores. Each score gets a stable PGS###### accession, structured metadata (trait, EFO/MONDO mapping, citation, training cohort, ancestry), and a downloadable scoring file.

Three file flavors exist per score:

Flavor	What it is	Do we use it?
Author-submitted	Exactly what the authors deposited. Build varies, may be rsID-only, columns inconsistent.	No.
`hmPOS` (harmonized positions)	Every variant mapped to a common build with `hm_chr`/`hm_pos`/`hm_rsID` columns added.	Yes — exclusively.
`hmVCF`	Harmonized and emitted as VCF.	No.

Our fetch path is PGSCatalogConnector.download_scoring_file(pgs_id, build="GRCh38"), which caches to <cache>/PGS######/PGS######_hmPOS_GRCh38.txt.gz alongside the REST metadata JSON. The full catalog is already on disk at C:\Maple\worker\data\pgs_catalog — 5,377 GRCh38-harmonized files, ~59 GB, bind-mounted into the workers. Do not re-download it. (DEFAULT_CACHE_DIR is a relative path, so it resolves differently per working directory — always pass cache_dir explicitly.)

3.2 File structure

A scoring file is a gzipped TSV with a #key=value header. Here is a real one in full — PGS000004, the Mavaddat 313-variant breast cancer score, annotated:

#format_version=2.0
##POLYGENIC SCORE (PGS) INFORMATION
#pgs_id=PGS000004
#pgs_name=PRS313_BC
#trait_reported=Breast Cancer          ← free text, as the authors described it
#trait_mapped=breast carcinoma         ← ontology term
#trait_efo=MONDO_0004989               ← the machine-readable trait ID
#genome_build=GRCh37                   ← the AUTHOR's build  (§2.6)
#variants_number=313
#weight_type=NR                        ← "not reported"      (§3.4)
##SOURCE INFORMATION
#pgp_id=PGP000002
#citation=Mavaddat N et al. Am J Hum Genet (2018). doi:10.1016/j.ajhg.2018.11.002
##HARMONIZATION DETAILS
#HmPOS_build=GRCh38                    ← the build hm_pos is on  (§2.6)
#HmPOS_date=2022-07-29

then a column header and rows:

chr_name  chr_position  effect_allele  other_allele  effect_weight  allelefrequency_effect  hm_source  hm_rsID  hm_chr  hm_pos     hm_inferOtherAllele
1         100880328     T              A             0.0373         0.4097                  liftover   (blank)  1       100414772  (blank)
1         10566215      G              A            -0.0586         0.329                   liftover   (blank)  1       10506158   (blank)
1         110198129     C              CAAA          0.0458         0.7755                  liftover   (blank)  1       109655507  (blank)

Four things are visible in three rows, and each has cost us something:

chr_position ≠ hm_pos. 100880328 vs 100414772 — a 466 kb shift from the GRCh37→GRCh38 harmonization. Reading the wrong column is a silent catastrophe.
hm_rsID is blank — and not just in these rows: all 313 of them. This score identifies its variants by position only. A pipeline that joins on rsID matches literally zero rows here and reports a confident zero.
Row 3 is an indel — C → CAAA, an insertion. Note the effect allele is the longer one; the shared anchor base C is what makes indel orientation unresolvable from the reference alone (§5.2).
allelefrequency_effect is present. This is the field that lets a scorer resolve a strand-ambiguous palindromic SNV — if the effect allele frequency is far from 0.5, the strand can be inferred.

Optional columns that change everything. Beyond the core schema, a scoring file may carry is_haplotype, is_diplotype, is_interaction, is_dominant, is_recessive, dosage_0_weight/dosage_1_weight/dosage_2_weight, variant_description, locus_name, OR, HR. Any of these means the file is not a plain additive sum, and quietly summing effect_weight over such a file produces a number that looks fine and means nothing. Detecting them is the classifier’s first job (§4).

Here is what one of those looks like — PGS004262, a GenoBoost psoriasis score:

chr_name  chr_position  effect_allele  other_allele  dosage_0_weight  dosage_1_weight  dosage_2_weight  ...
6         31274555      T              C            -0.0324954        0.1051552        0.2502743

There is no effect_weight column at all. The weights are per-genotype, and they are not linear in dosage (−0.032 → +0.105 → +0.250). An effect_weight-based parser reads zero rows from this file. Fifteen catalog scores are exactly this case, and the right conclusion is “needs a dedicated calculator,” not “the file is broken” — a distinction §4.4 encodes deliberately.

3.3 Scale and complexity

The catalog is not uniform. It is two populations of file wearing one format.

Bucket	Files	Bytes
> 50 MB	219	22.3 GB
10–50 MB	1,725	33.7 GB
1–10 MB	543	2.6 GB
< 1 MB	2,890	0.3 GB

Half the files are under a megabyte; 4 % of them hold 38 % of the bytes. The largest score by variant count is PGS003162 at 10,318,272 variants (five scores exceed 10 M); the largest by compressed size is PGS005197 — 9,503,677 variants in 185.8 MB. Those two rankings differ, which is a small but recurring trap: sizing work by file count, by bytes, and by rows gives three different answers. Across all 5,377 files the classifier examined 3,452,540,960 variant rows.

The split is historical, not arbitrary: pre-2015 scores are curated lists of genome-wide significant hits (tens to hundreds of variants, rsID-labelled), and post-2015 scores are LD-aware genome-wide posteriors (millions of variants, position-only). Any assumption that holds for one population fails for the other — which is why a 25-score smoke test once suggested weight_type was uniformly NR, a conclusion the full fleet demolished.

Cost model (measured, after optimization): ~0.98 s per MB of gzip → ~16.7 CPU-hours for the whole catalog. Sharded six ways on a 16-logical-core box that is ~2.9 h wall. Two non-obvious facts from profiling a 19.4 MB / 972k-row file:

The time is not I/O. gzip decompression plus readlines is 6 %. Tab-splitting is 8 %. Per-row dict construction plus PGSVariant allocation is 45 %, and feature accumulation another 41 %. The lever is Python object allocation, not decompression.
Therefore: shard by bytes, not by file count, and compute ETA from bytes. Under largest-first ordering, byte-progress and count-progress diverge wildly — at 1.2 % of scores the fleet run was already 12.1 % of bytes.

Cache integrity. Five of the 5,377 cached files were truncated by interrupted downloads (PGS002938, PGS002153, PGS000257, PGS003724, PGS004604); re-downloads were 25–62 % larger. A truncated gzip always fails its CRC/length trailer, so EOFError: Compressed file ended before the end-of-stream marker detects this reliably — a partial file can never be silently parsed as complete. Do not assume a cached scoring file is intact merely because it exists.

3.4 Weights and weight types

A weight is the per-allele effect size the study estimated. What scale it is on is declared (sometimes) in the #weight_type= header, and getting this wrong is a biological error the arithmetic cannot catch.

The scales you should expect:

Weight type	Meaning	Additive-summable?
`beta` (β)	Effect on a continuous trait, in the trait’s own units (or SD units).	✅ Yes — this is the natural additive scale.
`log(OR)`, `ln(OR)`	Natural log of the odds ratio, for a binary disease.	✅ Yes. Logs add; that is the point of taking them.
`log(HR)`	Log hazard ratio, from survival analysis.	✅ Yes, on the log-hazard scale.
raw `OR` / `HR`	The odds/hazard ratio itself.	❌ No. Ratios multiply. Summing an OR as if it were a β is meaningless — and produces a plausible-looking number.
`NR`	Not reported. The authors did not say.	⚠️ Unknown. Usually a β or log(OR), but that is an inference, not a fact.
`Unweighted`	Every variant counts equally — an allele count, not a weighted score.	⚠️ Computable, but it is a different object.

What the catalog actually contains (live from mods_prs_score_manifest, all 5,377 scores):

`weight_type`	Scores
`beta`	3,250
`NR`	1,390
`BETA`	441
`log(OR)`	199
`Inverse-variance weighting`	16
`Beta (corrected for sample overlap)`	12
`Odds Ratio over expected risk`	11
`Unweighted`	9
`ln(OR)`	9
`Beta`	8
`metaGRS`	6
`log(HR)`	4
`log(OR) …NOTE a EAF-normalised weight was used…`	3
`Dosage`, `metaPRS`, `posterior beta`, `MetaPRS Weight`, `PHS log(HR)`, `weighted by functional implication`, `weights.HC132/HC171/HC215/HC225/HC294/HC299`	1 each

Three things to take from that table:

Case normalization is mandatory. beta + BETA + Beta = 3,699 scores, or 69 % of the catalog — but only if you lowercase first. A case-sensitive match sees 3,250 and misclassifies 449.
NR is 26 % of the catalog and is not a scale. Roughly one score in four does not tell you what its numbers mean. Notably, 9 of the 13 scores in our own curated default set are NR (§4).
The tail is free text. weights.HC171 is not a unit. These need explicit per-score handling or explicit deferral, never a regex that shrugs.

The theoretical-bounds check. For any additive score, the raw total is mathematically bounded:

L = Σ min(0, 2·wᵢ)        U = Σ max(0, 2·wᵢ)

because every dosage is in {0, 1, 2}. Every computed score must satisfy L ≤ S ≤ U; a value outside that interval is a scorer defect, not a biological finding. It is free, it is computed for every score at classify time (stored as theoretical_lo/_hi), and it catches whole families of allele-orientation and arithmetic bugs.

It is also one-sided and scale-blind, which is worth internalizing: it cannot detect an odds-ratio-summed-as-a-beta, and it cannot detect a wrong-but-attainable score. We learned that concretely — a bounds check over the 40 legacy scores with |raw| > 1000 found all 40 inside their bounds, while every one of them was in fact wrong (§5.1). Bounds are a floor on correctness, not a ceiling.

3.5 Our internal classification

Given a scoring file, MODS assigns it exactly one compatibility class, plus the evidence for that assignment. This is Gate B of the rebuild, and it happens before any patient is involved — classification is a pure function of the file.

Class	Meaning	Outcome
A	Clean: additive, biallelic, autosomal SNVs, position-resolvable, explicit alleles, no complex-format columns.	Computed in production.
B	Conditionally supported: computable with care — indels, palindromic SNVs, chrX, inferred alleles, multiallelic sites, dominant/recessive models, some unplaceable rows.	Computed only if every one of its complications is on the promoted list.
C	Needs a dedicated calculator: haplotypes, diplotypes, interactions, per-genotype dosage weights (APOE/HLA/CYP star-allele style biology).	Deferred — “not calculated,” with a reason.
D	Unavailable: nothing in the file can be placed on the genome, or there are no usable rows.	”Not calculated.”

The signals the classifier reads, all structural:

Signal	Why it threatens the score
Unplaceable rows (no `chr`/`pos`)	Our genotype path is position-keyed; an unplaceable variant cannot be looked up.
Indel	The two alleles share an anchor base, so naive matching cannot orient them (§5.2).
Unresolvable allele	The allele names no single sequence — a harmonizer candidate set (`A/T`, `C/G/T`), an unknown base (`N`), presence/absence notation (`-`, `P`), or a star-allele name. Nothing to orient against the reference.
Palindromic SNV (A/T, C/G)	The strand complement looks identical, so a strand flip silently inverts the dosage. Resolvable only via allele frequency.
Multiallelic site	More than two alleles at one position; ordinary biallelic matching picks the wrong one.
chrX variant	Dosage is sex-dependent — a single-X person is hemizygous (0/1), not 0/1/2.
Inferred other-allele	The harmonizer guessed it. Not authoritative; trusting it can misassign reference vs. alternate and thus the dosage.
Missing other-allele	Cannot orient the variant against the reference at all, so a wild-type person’s dosage is undefined.
Dedicated-calculator columns	The score is not an additive sum.
Inheritance-model columns	The dosage → contribution mapping is not linear.

Promotion. Class B carries a second, signed-off decision: which complications we are actually willing to compute. A Class-B score is production-eligible iff every one of its categories is promoted — one un-promoted signal defers the whole score. The strict-AND is deliberate: it stops a mostly-simple score from sneaking through on the strength of its simple parts.

Class-B signal	Promoted?	Rationale
Indels	✅	The dense-callset builder emits normalized, reference-verified indel alleles.
chrX	✅	The current patient has two X chromosomes, so X dosage behaves like an autosome. A single-X patient would need hemizygous handling — this promotion is patient-specific.
Inferred other-allele	✅	Every allele is verified against the reference genome, so an inferred allele that matches is safe.
Palindromic SNV	✅ (resolvable subset)	Strand is resolved by allele frequency; genuinely ambiguous ones stay deferred at compute time.
Unresolvable allele (`ambiguous_other`)	❌	Added in classifier v3 (`GH #273`). Promoting it would mean deciding what `A/T` means — a mapping decision, not a threshold.
Multiallelic site	❌	Duplicate-position handling needs golden tests first.
Missing other-allele	❌	Cannot orient a wild-type dosage with no other allele.
Inheritance model	❌	Non-additive; out of scope for an additive scorer.
Unplaceable rows	❌	Partial coverage; handled by the coverage gate, not promoted here.

Re-classified under v3 on 2026-09-22. The numbers below are that run. Two fixes landed between the original 2026-07-21 classify and this one, and they move eligibility in the same direction for different reasons — see the attribution note after the blocker table. Eligibility fell 3,863 → 2,322.

Fleet results — the classifier has been run over the entire cached GRCh38 catalog, and these are live numbers from mods_prs_score_manifest:

Class	Scores	Share	Production-eligible	Variant rows
A	1,077	20.0 %	1,077	594,511,227
B	3,870	72.0 %	1,245	2,837,293,016
C	426	7.9 %	0	20,736,700
D	4	0.1 %	0	17
Total	5,377		2,322 (43.2 %)	3,452,540,960

So just over four in ten published scores are computable under the current policy, and the eligible set carries 894,528,461 variant rows (largest 7,082,943, mean 385,241) — the real size of the job facing the later gates. That is less than half the row count the v2 manifest implied, which roughly halves Gate E.

How common are the complications? n_palindromic > 0 in 3,124 scores, n_no_position > 0 in 2,629, n_indel > 0 in 1,709, n_chrx > 0 in 1,225, n_ambiguous_other > 0 in 996, n_inferred_other > 0 in 618, n_multiallelic > 0 in 567, n_missing_other > 0 in 509. Palindromes are the norm, not the exception — which is exactly why promoting them was the first policy decision worth making. Note n_indel fell from 2,300 to 1,709: 591 scores had no real indel at all, only unresolvable alleles a length test mistook for one.

What blocks the other 2,625, and — more usefully — what would be unblocked if each category alone were promoted:

Deferred category	Appears in	Sole blocker
Unplaceable rows	2,264	1,437
Unresolvable allele (`ambiguous_other`)	624	144
Multiallelic site	554	132
Missing other-allele	477	9
Inheritance model	3	1

Single-category promotions could unlock 1,723 of the 2,625 deferred scores (65.6 %).

⚠ This ranking inverted at v3, and the old advice is now wrong. The v2 table had multiallelic sole-blocking 588 and unplaceable_rows sole-blocking 18, so the standing recommendation was “promote multiallelic”. It is now the reverse: unplaceable_rows is the highest-leverage open decision by an order of magnitude (1,437 vs 132).

Attribution — two independent fixes moved these numbers, and conflating them would mislead:

GH #250 (merged 2026-07-26, five days after the original classify, so this is the first fleet run to include it) stopped the parser falling back from hm_pos to an original-build coordinate. Rows that used to be silently mis-placed are now honestly unplaceable. ~1,437 scores are sole-blocked by this. It is the dominant effect.
GH #273 (this change) split unresolvable alleles out of indel. ~144 scores are sole-blocked by this alone.

Note what “promote unplaceable_rows” would actually mean: computing a score while some of its loci cannot be placed — i.e. a partial-coverage policy, which belongs to the Gate F coverage decision, not to format handling. It is the biggest lever, but it is not the same kind of decision as the others in this table.

Class C is more homogeneous than expected: of the 426 scores, 378 are blocked by is_haplotype alone, 18 by is_interaction, 18 by the dosage_N_weight trio, and 12 by combinations. Class D is four scores, all tiny (1, 1, 5, and 10 variants — iritis, iridocyclitis, psoriatic arthritis, and SLE), where nothing at all could be placed.

3.6 The classification engine

Source: backend/mods_prs_compat.py · Tests: backend/tests/test_mods_prs_compat.py · Runner: worker task mods_prs_manifest_build, driven by backend/scripts/prs/build_prs_manifest.py.

Five design properties are worth understanding, because each was paid for:

1. One parser, one accumulator. There is exactly one row parser (pgs_catalog_connector._row_to_variant) and exactly one feature implementation (mods_prs_compat._FeatureAccumulator). Both the in-memory path (parse_scoring_file → score_features) and the streaming path (open_scoring_file → stream_score_features) drive that same accumulator, and both classifications land in the same classify_features. This is not tidiness — these counts decide production_eligible, so if the two read paths could disagree about a score’s features, they could disagree about whether it is safe to compute. Do not “optimize” by forking a second parse loop.

2. Conservative by presence. For complex-format and inheritance columns, the classifier flags on the presence of the column, not on confirming that rows use it. This over-defers a score carrying an unused column — the safe direction. Row-level refinement is a planned hardening step, not a bug.

3. Precedence is load-bearing. The order is D (unplaceable) > C (dedicated) > D (no usable rows) > B > A, and the placement of the second D check after C is deliberate: a dosage_N_weight file parses to zero effect_weight rows, and “needs a dedicated calculator” is the true explanation where “no usable rows” would falsely imply a broken file. Fifteen catalog scores hit exactly this path. Equally, that check must come before B and A — under classifier v1 an empty file had every feature count at zero, fell through to Class A, and was marked production-eligible: a score cleared for production that computes to nothing.

4. Versioned decisions. Every manifest row carries classifier_version and the scoring file’s sha256. A change in classifier meaning — or in what it extracts — bumps the version so a re-run re-classifies rather than skipping rows decided under different rules. v3 was the 2026-09-22 structural re-classify, v4 the unplaceable_rows threshold promotion (recomputed in place, no re-read), and v5 (2026-09-26, GH #377) adds the allele-convention counters below and therefore needs a full re-read.

A bump that is meant to change no decision should be made to prove it: scripts/prs/verify_classifier_bump.py snapshot before the re-run and compare after it diff every score’s class, categories and eligibility. v5 moves no counter and promotes no category, so its diff must be empty — and an empty diff is also how you know the re-run actually ran.

4b. The allele-convention profile (v5). n_snv_pair_top (A/C + A/G) and n_snv_pair_bot (C/T + G/T) split a file’s non-palindromic biallelic SNVs in two. On the forward strand those are the same variants seen from the two strands and occur about equally often, so top_fraction = top/(top+bot) sits at ~0.5 — measured 0.4991, 0.5017, 0.5018 on three ordinary scores. ~1.0 or ~0.0 means the alleles were written in an A/B convention (Illumina TOP or BOT), whose strand assignment is made per variant and recorded nowhere in the file, so about half its rows are reverse-complemented relative to GRCh38 with nothing saying which half.

Counts, never the ratio — mods_prs_compat.top_fraction derives it on read, the same rule Gate C’s coverage payload follows. And the profile alone never condemns a file: PGS005087/5088/5089 are top_fraction = 1.000 with a 0 % complement-only rate, their alleles being the forward-strand pair at 5,326 of 5,326 checkable loci. The guard needs the measured complement-only rate too (§5.4). Full reasoning: `prs_validation/palindrome_orientation_spike.md`.

5. Memory discipline at catalog scale. parse_scoring_file holds the decompressed text, the split-line list, and every PGSVariant at once — fine for one file, fatal across 59 GB. The streaming reader yields lazily. Two further optimizations, each verified behavior-identical against a verbatim reimplementation across 16 real files (72 → 1.1 M variants, classes A/B/C, zero mismatches):

_ColumnIndex resolves column positions once per file, so _row_to_variant indexes the split fields rather than building a dict per row — a ~1.3× speedup.
Allele pairs are interned to small integers before per-site storage. A scoring file has a handful of distinct allele pairs but millions of sites; keeping the pair strings per site dominated peak memory at 1,040 MB on the 9.5 M-row file. Interning is exact — distinct pairs get distinct ids — so unlike hashing it cannot under-count multiallelic sites, which would silently promote a score to eligible.

Cross-validation is still open. The classifier is deliberately structural and dependency-free, which means it is our reading of the format rather than the reference implementation’s. GH #237 tracks cross-validating it against the official pgscatalog libraries.

4. Aside: the curated 13, and why defaults are dangerous

SMOKE_TEST_PGS_SET in pgs_catalog_connector.py (named UKB_TRAIT_TO_PGS_DEFAULTS until 2026-09-28) maps 13 UK-Biobank trait codes to specific PGS IDs. It was assembled best-effort in 2026-05, before any genotype data existed, and its header comment claimed it was “not load-bearing.” It was in fact load-bearing in three places: it sized the KGIA/PDR2 disease axis, it supplied the PDR1/PDR2 export CUIs, and it was the fallback for mods_prs_compute.

Checking those 13 IDs against the actual scoring-file headers turns up something worse than “unvalidated” — three of them are for the wrong disease entirely:

Code	Label in our code	Actual `trait_reported` in the file
RA	Rheumatoid arthritis	Coeliac disease (`GRS-DQ2.5-CeD`)	❌ wrong trait
MS	Multiple sclerosis	Venous thromboembolism (`PRS_VTE`)	❌ wrong trait
PD	Parkinson’s disease	LDL cholesterol (`GLGC_2021_SAS_LDL_PRS_weights_PRS-CS`)	❌ wrong trait
HT	Hypertension	Systolic blood pressure (`GRS970_SBP`)	⚠️ proxy, not the same construct
CRC, AD, T2D, CAD, AF, ISS, BC, PC, AST	(as labelled)	matches	✅

The PD entry is doubly wrong: besides being a lipid score rather than a neurological one, GLGC_2021_SAS is trained on a South Asian cohort, which is the worst-case portability choice for a European-ancestry patient (§1.3).

Separately, 5 of the 13 are not production-eligible under the Gate B classifier (PGS000016, PGS000019, PGS000041, PGS000511, PGS000896), and 9 of the 13 carry weight_type=NR — so for most of the set we do not know what scale its weights are on.

The decision of record (GH #244, GH #95) is not to delete the set — a small fixture is genuinely needed, because you cannot do a test or export run against a 5,377-score manifest. The shape agreed is: keep it, rename it to say what it is (done: SMOKE_TEST_PGS_SET, 2026-09-28), and derive the production/training axis from the release instead: axis v2 is one position per released Gate F axis pick (#244). Never let a hand-picked literal be a silent default. Note also that these 13 are relevant-to-Maureen but not structurally representative of the manifest, so they make a fine clinical shortlist and a lopsided pipeline-stress fixture.

The general lesson, and the reason this section exists in a primer: a default that is never validated becomes a specification. This one silently defined the disease axis of a machine-learning model for two months.

5. The processing engine

5.1 The bug that defines the design

Everything below is shaped by one defect, found by external audit on 2026-07-12, that invalidated all 5,368 previously computed scores.

The old genotype store was variant-only: the VCF importer skipped reference blocks, so a position where the patient matches the reference has no row. The scorer then did the obvious thing — a scoring locus not found in the store was counted as missing and contributed nothing.

That is correct only when the effect allele is the alternate allele. When the PGS effect allele is the reference allele — and a person who is homozygous-reference therefore carries two copies of the effect allele — the correct dosage is 2 and the correct contribution is 2 × weight. The pipeline scored 0.

Measured on PGS000019, ~49 % of loci are REF-effect (25 REF vs 26 ALT). At those loci the dosage map became {0, 1, 0} instead of {2, 1, 0} — an “artificial heterozygote.” This is a bias, not a constant offset: it distorts different scores by different amounts depending on how many of their effect alleles happen to be reference.

Three lessons were extracted, and they are worth more than the fix:

“Absent” is ambiguous in a variant-only store. It can mean hom-ref or not sequenced, and those must never collapse to the same number.
snps_used / snps_in_scoring_file is not coverage. That ratio clustered at 42–52 % across the legacy fleet, and the retired UI colour-coded it as a data-quality verdict about the published score — presenting a symptom of our own defect as somebody else’s problem.
A self-check that shares your assumptions cannot fail. The internal reconstruction that “confirmed” the ~50 % coverage was benign assumed absent = dosage 0, so it reproduced the bug instead of detecting it. The audit caught it by asking the one question the check could not.

The response was to wipe rather than quarantine: 9,793 rows deleted on 2026-07-13 (5,368 distinct scores plus superseded re-runs), backed up first to backend/data/prs_backups/. Scores are fully recomputable from the retained genotypes plus the scoring files, so keeping known-wrong numbers around bought nothing and risked a leak. The 8,088,434-row mods_genotype table was untouched.

5.2 Dense callsets

Source: backend/connectors/wgs_dense_callset.py · Locus union: backend/connectors/pgs_locus_union.py · Worker task: genotype_dense_callset

The fix is to stop treating “absent” as “0” by making absence impossible at the sites we care about. The mechanism is a dense callset: a genotype record for every locus some target score needs, whether or not the patient carries a variant there.

“Dense” does not mean all 3.1 billion positions. It means dense over the union of PGS loci — a site-targeted reconstruction, typically a few million positions, built from the source gVCF rather than from the variant-only database. This is deliberate: build_locus_union() collects the (chrom, pos, alleles) sites the target scores actually reference, deduplicated, and the builder resolves exactly those.

Every target locus resolves to exactly one of four outcomes:

Outcome	Condition	Emitted
`CALLABLE_HOMREF`	A confident reference block spans the coordinate.	`GT 0/0` with a concrete REF (from the FASTA) and the scoring non-ref ALT — so a REF-effect locus becomes dosage-2 eligible.
`CALLABLE_VARIANT`	An explicit, quality-passing variant record.	Its genotype.
`EXPLICIT_MISSING`	A no-call, a FAILed record, below-threshold quality, or no record at all.	`./.` — routed to the missingness policy, never silently dropped.
`STRUCTURAL_FAIL`	REF disagrees with the FASTA, contig absent, orientation unresolvable.	Nothing. Logged, excluded, never a partial score.

And the Gate A exit invariant is a hard assertion, not a log line:

callable_homref + callable_variant + explicit_missing + structural_fail == union_size

i.e. zero residual implicit-zero. Every locus is accounted for by name. If the arithmetic does not reconcile, the build raises rather than returning a plausible result.

The callability policy, and the discovery that nearly re-introduced the bug. The patient’s gVCF is Sequencing.com 30× WGS, reblocked by GvcfReblocker — GATK-style (END/MIN_DP reference blocks, not DRAGEN RGQ), Ensembl contig naming. The provisional plan said hom-ref calls should require GQ ≥ 20, which is standard advice. Measuring the real file showed every PASS reference block carries GQ = 0 — 0 of 159,540 in the first 300k records had GQ > 0, because the reblocker does not preserve per-site GQ across collapsed blocks. Gating hom-ref on GQ would have rejected every hom-ref call and silently re-created the exact defect being fixed.

So the policy identifies the caller from the header first, then applies thresholds:

	Rule
hom-ref callable	ref block spans coord, `FILTER == PASS`, `GT == 0/0`, `MIN_DP ≥ 10`, `DP ≤ 3 ×` autosomal median, REF verified against the FASTA. No GQ gate.
variant callable	`FILTER == PASS`, `GQ ≥ 20`, `DP ≥ 10`, REF == FASTA. Indels stricter: `DP ≥ 15`, `GQ ≥ 30`.
explicit missing	`FILTER == FAIL` or `GT == ./.`

The depth ceiling rejects pile-up-inflated collapsed/repetitive regions. A strict() variant (MIN_DP ≥ 15, variant GQ ≥ 30) is emitted as a second callset for sensitivity analysis. The numeric floors remain provisional pending clinical sign-off.

The indel-orientation trap. For equal-length alleles you can ask the FASTA which allele is the reference — exactly one matches. For an indel you cannot. A deletion ATTCT → A has an ALT that is the anchor base of the REF, so a per-allele “does this equal the reference span of its own length?” test matches both alleles. From a FASTA plus an unordered allele pair, “deletion of TTCT” and “insertion of TTCT” are genuinely indistinguishable. The first implementation structural-failed all three real indel sentinels with a misleading reason, and would have quietly dropped every indel hom-ref locus in a catalog-scale run. The fix: TargetLocus carries explicit ref/alt when the source knows it, and the FASTA only verifies (ref == ref_bases(pos, len(ref))) rather than deriving. FASTA-derivation is now gated on len({len(a) for a in alleles}) == 1. Orientation-unknown indels defer honestly.

Real-data validation (2026-07-13, against the actual gVCF, 6 sentinel loci, reconciles 6/6, zero structural failures):

Autosomal median depth 32 → depth ceiling 96.
A known pathogenic variant, carried in the source report with no rsID at all — resolved correctly by coordinate, so the position-keyed builder catches exactly what the legacy rsID-keyed path dropped.
Five somatic tumor loci (BRAF V600E, NRAS Q61K, three indels) → germline hom-ref GT 0/0, hence REF-effect-eligible.

Both results are biologically correct: a carrier call at the germline site variant, wild-type in the germline at every somatic locus. That is the check that would have failed loudly if the reconstruction were wrong.

One boundary worth guarding. The dense-callset path does not read mods_genotype. It streams the local gVCF directly against the pinned local FASTA. That is the whole point — reading the source file is precisely what bypasses the variant-only store that caused the bug. Do not “helpfully” wire the builder onto mods_genotype; it would reintroduce the defect. A pleasant side effect is zero Cloud SQL egress, since the only non-patient input is a public reference genome on local disk.

5.3 Streaming versus loading in memory

This pipeline touches four things that do not fit in RAM, and each needed the same treatment. The pattern is uniform enough to state as a rule: hold one record, not one file.

Where	Naive approach	What it cost	The streaming form
gVCF ingest	`read_bytes` → `gzip.decompress` → `splitlines` → list	Multi-GB, OOM on a WGS gVCF	`iter_vcf_records()` yields lazily — peak 0.9 MB to walk ~5 M variants
Scoring-file classify	`parse_scoring_file` holds text + split lines + every `PGSVariant`	1,040 MB peak on one 9.5 M-row file; impossible across 59 GB	`open_scoring_file` yields `(meta, rows)` with `rows` a generator
Dense callset	Index the gVCF, random-access each target	Requires pysam/cyvcf2 (neither installs on the Windows dev box)	Streaming merge-join: both sides sorted by (contig-rank, pos), one linear pass, bounded memory, pure-Python gzip
Genotype cache	Re-stream 4.4 M rows from Cloud SQL per task	Billed egress, thousands of times, for a file that never changes	Pull once → pickle to shared local disk → every later task reads locally, zero egress

Two footguns that only show up at scale:

Chunk your commits. The first 4.4 M-row ingest accumulated everything into a single transaction and committed at the very end. A uvicorn --reload (triggered by staging a 200 MB VCF under the watched mount) killed the long-lived transaction and rolled back everything — hours gone, zero rows durable. The rewrite commits ~100 K records at a time, each its own short transaction, with per-chunk progress and a printed resume point. The upsert is idempotent, so a re-run resumes cheaply.

Streaming the file does not mean streaming the working set. The genotype cache solved the egress problem and created a memory one: each task unpickles the cache into a ~2 GB in-RAM dict. Five concurrent local workers on the 30 GB box hit ENOMEM, crashed, and churned their tasks; two workers ran clean and finished all 36 tasks. Fewer workers finished faster, because the only thing quicker than a fast worker is one that never has to restart. (2026-09-29: history. The cache, get_patient_genotype_dict_cached, served only the retired compute_prs_pure scorer and was deleted (GH #481). Gates A, C and D read the dense callset, never mods_genotype (§5.2), so do not revive it.)

One more, from skipping records: in the VCF importer, skip reference blocks by testing that ALT is a non-variant (., <NON_REF>, *) — not by testing GT == 0/0. A genuine hom-ref call at a real variant site has a real ALT and must be kept. Keying the skip on the genotype silently drops real variants; there is a regression test (test_parse_vcf_homozygous_ref_emitted) holding that line.

5.4 The scorer (Gate C)

Source: backend/mods_prs_scorer.py · Tests: backend/tests/test_mods_prs_scorer.py · Persistence: mods_prs_service.store_score_result · Worker task: mods_prs_compute

Gate A fixed how the genome is read; this is the arithmetic that reads it. The scorer never asks “is there a variant row here?” — the question that conflated homozygous reference with not sequenced. It asks, of the dense callset, which allele is the reference at this locus, and how many copies of the effect allele does the genotype name? A GT 0/0 at a REF-effect locus answers two, and that single line is the fix.

Orientation is resolved from the reference allele alone, not from REF+ALT. That is deliberate: Gate A emits ALT=. for an explicitly-missing locus and the scoring ALT for a hom-ref one, so a rule that needed the ALT would behave differently on the two.

Every scoring row lands in exactly one of eight outcomes, and the scorer raises rather than returning a number if they do not sum to the rows it read — the same exit gate Gate A applies to its four:

Outcome	Meaning
`scored`	Dosage resolved and added.
`not_in_callset`	No Gate A record at this coordinate — a gap in the targets.
`explicit_missing`	Record present, genotype is a no-call — a gap in the sequencing.
`ambiguous_palindromic`	A/T or C/G pair; strand unresolvable from alleles.
`ref_allele_mismatch`	The reference allele is neither scoring allele.
`unsupported_ploidy`	Genotype is not diploid (the chrX promotion assumes two X’s).
`unplaceable_row`	No harmonized position (§2.6).
`missing_other_allele`	Cannot orient a wild-type dosage with one allele.

The three “not scored for cause” outcomes are kept apart from the two “gap” outcomes because they mean different things to a reader; collapsing them into one “missing” bucket is how the retired coverage column came to mean nothing (§5.1).

Two conservative defaults, both recorded on the run rather than baked in.

Palindromic SNVs defer. The complement of an A/T pair is an A/T pair, so both strand readings fit the data and they give opposite dosages — the one class where a strand error is silent and inverting. Resolving it needs an allele frequency to compare against, i.e. a reference panel: Gate D. Palindromes appear in 1,197 of 2,322 scores (51.6 %) and account for 13,929,705 of 894,528,461 rows — 1.56 % overall.

But the aggregate badly understates the per-score cost, and this document said so for a day. Measured on the Gate E1 run (mods_polygenic_score, 2026-09-24), among the 1,197 affected scores the deferral takes a median of 10.24 % of their rows, p95 21.49 %, mean 11.20 % — and 1,088 of 2,321 scores lose more than 5 % of their rows. The aggregate 1.56 % was right; the per-score figures previously quoted here (median 2.47 %, p95 17.33 %) were computed on a different denominator and are 4× too low on the median.

The reason is visible the moment you split by score size, and it inverts the “concentrated in a tail” reading:

score size	scores	mean % of rows palindromic	median
< 1k rows	1,193	9.21 %	8.86 %
1k–100k	452	4.60 %	0.00 %
100k–1M	293	0.00 %	0.00 %
> 1M	383	0.87 %	0.00 %

The burden sits almost entirely on the small, curated, pre-2015 hit-list scores — which are precisely the interpretable ones — while the genome-wide LD-posterior scores that dominate the row count carry almost none. So “1.56 % of rows” is a statement about the big scores, not about the fleet. Gate D’s palindrome recovery is worth substantially more than a row-weighted average suggests.

Two earlier revisions got this wrong in opposite directions. One quoted “12.5 % of rows, median 7.6 % per score” from a 60-score sample; the correction to full-fleet aggregates then over-corrected the per-score claim. The durable lesson is that a row-weighted aggregate and a per-score median answer different questions, and for anything that is consumed per score — which a feature vector is — the per-score distribution is the one that matters.

Strand flips are measured, not applied. Both sides are GRCh38 and Gate A has already verified every record’s REF against the pinned FASTA, so a complement-only allele match is evidence of a problem — an author-strand submission, a harmonization error, a different variant — rather than a coordinate detail to repair in passing. Applying the flip would make the problem permanently undetectable: for a single-base pair, accepting complements resolves all four possible reference bases, so ref_allele_mismatch could never fire again and the scorer would lose its only structural check that the scoring file and the genome describe the same variant. So the rows a flip would have rescued are counted as n_complement_only_match and deferred. Near zero confirms the forward-strand premise; a large count refutes it. allow_strand_flip=True exists, and says so on the run.

Read that “near zero” per score, not fleet-wide (GH #371, 2026-09-26). E2’s aggregate is 0.0226 %, but 37 scores sit above 5 % and six near 50 % — the signature of a file written in Illumina’s A/B allele convention, where the strand is assigned per variant and recorded nowhere. PGS005393 is the clean case: 53,705 rows, every one A/C or A/G, not a single C/T, and 49.9 % of them reverse-complemented. For those files the premise is false, half their ordinary rows are being dropped as ref_allele_mismatch, and their palindromes are the least safe in the fleet.

Those rows are recoverable, per file and on evidence (GH #376). ScoringPolicy.strand_flip_recovery_scores names the scores whose complement-only rows may be flipped — a set of ids rather than a second boolean, because the global allow_strand_flip is an assumption about the whole catalog and this is a finding about named files, and a run’s provenance must not let the two look alike. A file qualifies on both signals: an A/B allele convention (§3.6 4b) and a measured complement-only rate above 5 %. scripts/prs/strand_recovery_candidates.py computes the set and refuses to call it final while any manifest row is still at an older classifier. Measured over the full v5 manifest (2026-09-26): 5 scores, 26,885 recoverable rows — 4 further files show the allele convention but match the reference forward, and 32 exceed the rate without the convention, and neither group is a candidate. Empty by default — it recovers coverage, not correctness, it cannot touch a palindromic row, and switching it on means re-scoring the affected scores, so it is a Gate F decision that can ride E3.

Weight scales are bucketed, not guessed (§3.4). beta/log(OR)/ln(OR)/log(HR) are additive and computed. A raw odds or hazard ratio is refused — summing ratios is a category error that produces a plausible number inside the theoretical bounds, so neither the arithmetic nor the bounds check would catch it. Everything else — NR (1,390 scores) and the free-text tail — computes with a weight_scale_unverified flag, which Gate F can gate on. The allow-list is exact strings rather than a regex, because a pattern loose enough to catch metaGRS would be inventing a scale for a quarter of the catalog.

What it writes. A coverage_v1 metadata payload carrying the three named quantities of §5.4’s read contract, plus a gate_c block with the outcome histogram, the dosage histogram, the retained and theoretical bounds, and n_homref_effect_dosage_2 — the count of loci a hom-ref genotype scored as dosage 2, every one of which the old scorer scored as 0. No ratio is ever stored: coverage fractions are computed on read from numerator and denominator, so storage and display cannot disagree. No percentile_v1 key is written at all, because Gate D owns it and an absent key renders as “the scorer did not record a percentile”, which is the true statement.

Rows are written at validity='quarantined', one commit per score, superseding any earlier unreleased score for the same (patient, PGS). Measured throughput on the largest catalog file (PGS005197, 9,503,677 rows): 220 s at a 2.2 MB peak heap — the streaming reader holds one record, where the in-memory parser peaked at 1,040 MB on that same file (§5.3).

Not yet run: the pgsc_calc cross-validation. scripts/prs/reconcile_pgsc_calc.py exports the bundle and diffs the two sums, but the middle step needs Nextflow and a container runtime.

It is not a one-time gate. The decision of record (Richard, 2026-09-22) is that both calculators run on every scoring pass and the delta is kept as a first-class output. A check that shares the scorer’s assumptions cannot fail (§5.1), so the structural answer is a standing second opinion rather than a snapshot — and the divergence is itself informative, because it decomposes into three named terms:

Term	Source	What it tells you
imputation	`pgsc_calc` adds `2p` for missing dosages; Gate C drops them	how much of the published score is invented rather than measured
palindrome	`pgsc_calc` resolves strand by allele frequency; Gate C defers	what the deferral in this section actually costs, per score
residual	neither policy explains it	should be zero — a non-zero residual is a real defect in one of the two

To isolate the residual, run pgsc_calc with --no_imputation --keep_ambiguous false; the remaining comparison is over the intersection of loci both tools scored, which is the part that tests the arithmetic. Running it without those flags is the other useful mode, because the difference between the two runs is exactly the imputation + palindrome terms.

MODS stays canonical for the released value — it is the only one that emits the coverage_v1 payload the read layer renders, and it reproduces in-repo without a container. pgsc_calc is the witness, and agreement within tolerance becomes a Gate F release criterion.

5.5 Current status

The rebuild is organized as lettered gates. A, B, C and the two E passes are done; D and F have run for the v0 set only (smoke run #1, 2026-09-28): the Gate F v0 release writer is backend/mods_prs_release.py, driven by backend/scripts/prs/release_gate_f.py (GH #296; --demote withdraws a release), and it released 590 rows tagged smoke_v0. The final, fleet-wide D and F have not started.** The single most important thing to understand is that “we computed scores” does not mean scores exist to a reader — every row is quarantined until Gate F promotes it, and the read layer filters on validity.

Gate	What it does	Status
A — dense callset	Rewrites how the genome is read, so a hom-ref site says `GT 0/0` instead of being absent. Fixes the input.	✅ Live, and the pinned GRCh38 FASTA is staged in the canonical checkout as of 2026-09-22 (§2.5). Migration 103; `mods_prs_run` holds one real run (`708c70ba…`, gate `gate_a_dense_callset`, reconciles 6/6, full provenance). Merged via `PR #235`.
B — compatibility manifest	Catalogs which published scores are computable, and why. Never touches the genome. Decides what is worth computing.	✅ Live. Migrations 104 + 107 + 108; all 5,377 scores classified, 3,753 production-eligible at v4. Re-classifying on `gate_b_structural_v5` (2026-09-26, GH #377) to add the allele-convention counters — decisions unchanged by construction, asserted with `verify_classifier_bump.py`.
B.5 — read layer + UI	Projects the manifest and runs into an API and five dashboard tabs (Fleet, Manifest, Policy, Runs, Scores).	✅ Live. `mods_prs_pipeline_service.py` / `mods_prs_pipeline_api.py`, `frontend/src/components/prs/`.
C — the scorer	Replaces `compute_prs_pure` with a scorer over the dense callset. Fixes the arithmetic.	✅ Built. `backend/mods_prs_scorer.py`; `compute_prs_pure` is deleted, not deprecated. `test_hom_ref_effect_allele_contributes_two_copies` is green, and the exit gate is a reconciliation assertion (every scoring row lands in exactly one of eight outcomes) plus the theoretical-bounds check. Cross-validation against `pgsc_calc` is exported and diffed by `scripts/prs/reconcile_pgsc_calc.py` but has not been run — it needs Nextflow + a container runtime.
D — ancestry & percentile	Normalize raw scores against a reference distribution, with continuous PCA adjustment that places exact coordinates rather than bucketing. Panel: `pgsc_HGDP+1kGP_v1` (user, 2026-09-24) — the same panel the `pgsc_calc` witness uses, so a MODS-vs-witness divergence stays attributable to arithmetic rather than to panel choice.	🚧 Built; first real pass 2026-09-28 on the v0 set (D1 `61aaccdf`, D2 `d63b3439` under the #499 fix, D3 against EUR). Not run fleet-wide; ancestry PCs not computed. Before that: built, not run. The archive is downloaded, md5-verified and extracted (§5.7): 15 GB of plink2 `.pgen`/`.psam`/`.pvar.zst`, 3,942 samples over 84,325,213 variants. D3 (PC fit, percentile, `run_gate_d.py`) was built in PR #386; D1 (`mods_prs_panel_af`) and D2 (`mods_prs_panel_score`) in GH #396, tested on synthetic plink2 fixtures (§5.7, “as built”). The real pass waits on the worker-image rebuild that pins `pgenlib` + `zstandard` (GH #339), and D2 on the E3 re-score (GH #422). Not optional, and not merely “after E” — it is what makes the feature vector usable at all (§1.4).
E1 — the measurement run	Compute numbers for the 2,322 v3-eligible scores, to settle policy.	✅ Done 2026-09-24 (§5.6). Union 11,400,292 loci; 2,321 scored, 1 refused by design; all reconcile, all within theoretical bounds. Quarantined. Superseded by E2.
E2 — the v4 fleet run	Recompute at v4 scope: 3,753 eligible scores, 2.166B rows.	✅ Done 2026-09-26 (§5.8). Union 19,215,298 loci → 16,964,243-record callset; 3,752 scored, 1 refused (the same `PGS000156`); all reconcile, none outside bounds, and the outcome counts sum to the manifest’s eligible rows exactly. All E1 rows now `superseded` — no mixed-callset blend remains. Quarantined.
F — biological checks + release	Weight-type sanity (§3.4), coverage gates, then the PDR feature manifest.	🚧 v0 writer built (#296) and run: 590 released as of 2026-09-28 (590 scores over 109 anchors, `smoke_v0`, #296’s rules with recorded waivers; `backend/mods_prs_release.py`, `scripts/prs/release_gate_f.py`). The full release is not started. It now has its inputs: 927 scores carry an unverified weight scale, and the palindrome/`NR` calls are the remaining policy.

The running order is E → D → E again, and the letters hide that. The gates are named in dependency order for contents, not for execution. Two things break the straight-line reading:

D is itself a fleet run. A percentile needs a scored reference distribution — thousands of panel genomes pushed through the same Gate C policy, or the deflation does not cancel between Maureen and the panel (§1.4). Gate E is the instrument Gate D uses, not its predecessor by alphabet. It is also a different data shape: the dense callset measured 8.45 GB for one genome over the 19,215,298-locus v4 union (498 B/record — and note that scaling E1’s measured 407 B/record would still have been 22 % low, §5.3), so a panel does not fit the Gate E layout and the architecture has to be designed against a measured single-genome run.
E is the measurement that settles the open policy decisions below. The unplaceable_rows threshold, the palindrome cost, and the NR weight-type policy were all judgement calls. E1 turned the first into a number and E2 confirmed it at scale; the palindrome question is now the only one that still gates a re-score, and it is answerable only from the panel. Normalising first would calibrate against a policy not yet chosen.

So the order of work is:

Pass	What	Released?
E1	Maureen, one pass over the 2,322 v3-eligible scores. A measurement run — its numbers were disposable by design.	✅ done 2026-09-24 — `quarantined`, now superseded
—	Settled `unplaceable_rows` on E1’s actual histograms → classifier v4, 3,753 eligible. Needed nothing from the panel.	✅ done 2026-09-24
—	Rebuild the union + callset at v4 scope. Must precede D, because D1 extracts frequencies over the union.	✅ done 2026-09-25
E2	Re-score at v4 scope. Ran before D rather than after it — so it is a second measurement run, not the “final policy” run the row below originally meant.	✅ done 2026-09-26 — `quarantined`
D1	Stage the panel; extract per-locus allele frequencies over the union. Built (GH #396), not run.	No
—	Decide palindromes on D1’s frequencies: leave deferred (E2 is the release candidate) or resolve them (Maureen needs an E3).	—
D2	Score the 3,330 unrelated panel samples under the final policy → the reference distribution. Built (GH #396), not run.	No
D3	Fit ancestry PCs + the mean/variance regression → z-score and percentile. Built (PR #386).	No
E3	Only if the palindrome decision changes the policy.	Candidate
F	Weight-scale sanity, coverage gates, then promote.	Yes

The trap is not doing E first. The trap is doing E once and treating its output as final. E1 → v4 → E2 is exactly that loop running as designed; whether it closes here or needs one more turn is the palindrome decision above.

The last hop to the model is now wired (GH #297, 2026-09-25). build_pdr2_input_pickle.py used to read float(row.raw_score) — not z_score, not percentile — so Gate D’s output would have had nowhere to go. resolve_prs_vector now reads the normalized field and records which one on the pickle, its .meta.json sidecar, the /prs/export/pdr2 headers and the /prs/preview/pdr-input payload. Two properties of that switch matter more than the field change. One field backs the whole vector: a z-score at one position and a 0–100 percentile at another are not one scale, and the loss reads the positions together, so a released score missing the chosen field drops out (its position takes the reference fill, measured says False) rather than contributing on a different axis. And the raw fallback is opt-in — with no normalized field anywhere the export raises UnnormalizedPRSError (409 from the API) unless --allow-raw-prs-fallback asks for raw sums, which are then flagged prs.raw_score_fallback.

Superseded 2026-09-28: 590 released rows (smoke_v0, release run c3be9df0). The paragraph below describes the state before the v0 release.

Therefore, right now: mods_polygenic_score has 0 released rows, and will keep having 0 until Gate F promotes — Gate C writes quarantined. mods_genotype has 8,088,434 rows and is no longer on the scoring path at all. The application correctly displays “no PRS computed,” and nothing downstream — API, MCP tool, or PDR export — can surface an unvalidated score, because the reader functions filter on validity == 'released'.

A Gate C row is rendered as withheld (“computed, but the result is not released”), which is a different and more informative statement than the awaiting_computation the same screen showed before. That distinction is the whole reason the read layer discriminates five states instead of rendering a blank cell.

The release lifecycle that enforces that: migration 102 added a validity column plus provenance columns and a partial unique index (uq_prs_released), and all_scores_for / latest_score_for / _build_prs_vector / _build_kgia_inputs_summary all filter on it. A rebuilt-but-unvalidated score cannot leak into a model input.

How the UI states this honestly. mods_prs_pipeline_service.py discriminates five score states — released, withheld, deferred, awaiting_computation, unclassified — because “we deliberately will not compute this, here is why” and “Gate C has not run yet” are opposite claims that a blank cell would render identically. Numeric keys are physically absent, not null, whenever a number is not authorized, so a caller that forgets to check state gets undefined and renders nothing rather than a misleading zero. And coverage is three separate named quantities (structural overlap / genotype callability / effective variance-weighted information), never one scalar — the single “Cov.” column is what encoded the REF-effect bug as a verdict.

Open decisions, ranked by leverage:

~~Decide the partial-coverage policy for unplaceable_rows~~ — done 2026-09-24, classifier v4 (user, this session). It turned out to be barely a decision. The 1,437 sole-blocked scores are blocked by a median of 0.005 % of their rows (p95 0.138 %, worst 5.0 %), while carrying 1.27 billion variant rows — more than the entire v3 eligible set. A threshold of 0.1 % would admit 1,338 of them, 1 % admits 1,431, 5 % admits all 1,437: the curve is flat, so the threshold buys no coverage and exists only to refuse a score with a genuinely meaningful hole. MAX_UNPLACEABLE_FRACTION = 0.01 (a score must be ≥ 99 % placeable). Eligible: 2,322 → 3,753; eligible rows: 894.5M → 2.166B. No score lost eligibility.

Two things make this cheaper and safer than it looks. Unplaceability is a property of the scoring file, not of the patient’s genome, so it deflates the patient and the reference panel identically and Gate D’s normalization — a monotone transform — absorbs it exactly (§1.4); the cost is variance, not bias. And because v4 moved no feature extraction, eligibility is a pure function of counts the manifest already stores, so backend/scripts/prs/reclassify_v4_eligibility.py recomputes it without re-reading 59 GB. That shortcut is verified, not assumed: --verify N re-parses N scoring files end to end and checks a real classification agrees (30/30, 0 mismatches).

Promoting multiallelic is the next format decision and is worth only 132 scores.
~~Close the hm_pos → chr_position fallback~~ — done (GH #250, §2.6). The live successor, and now the only open decision that can still force a re-score: resolve palindromic SNVs. Gate C defers them. Restated at E2’s v4 scope (2026-09-26, over the 3,752 computed scores): 1,687 scores carry at least one palindromic row, and among those the deferral takes a median 10.41 % of the score’s rows (mean 11.05 %, p95 19.70 %), with 1,575 scores losing more than 5 % — up from 1,088 of 2,321 at E1 scope, because v4 admitted exactly the small curated pre-2015 scores where palindromes concentrate. The fleet-wide row figure is 64,810,321 / 2,166,261,879 = 2.99 %, and as always that row-weighted aggregate describes the big genome-wide scores, not the per-score burden (§5.4). Their strand cannot be read off the alleles, and only 8.3 % of scores carry an effect allele frequency of their own — so the scoring files cannot resolve this themselves and Gate D’s reference panel is the only route. Decide this on D1’s frequencies, before D2: resolving palindromes changes the scoring policy, and therefore obliges an E3.

Scoped 2026-09-26 by GH #371 (§5.7b): there is no second route — hm_source carries no strand information — and the frequency route yields 8.48 % of the deferred rows, not 9.5 %, once frequencies within 0.08 of 0.5 are excluded as undecidable. So what is left to decide is where to set that margin; the residual ~91 % becomes a Gate F coverage gate rather than a pending question.
~~Re-classify the fleet under gate_b_structural_v3~~ — done 2026-09-22 (GH #273, 6 shards, ~103 min, 0 failures, all 9 manifest invariants pass).
~~Retire or rename UKB_TRAIT_TO_PGS_DEFAULTS~~ — renamed SMOKE_TEST_PGS_SET 2026-09-28 (#244); it no longer sizes the disease axis. Its three wrong-trait entries are documented on the constant, not fixed, since it is now a named fixture (#95) (§4).
Normalize weight_type case, and decide the NR policy for 1,390 catalog scores (§3.4). E2 gives this its operative number: 927 of the 3,752 computed scores carry a weight_scale_unverified flag (E1: 607). Gate F’s call; it gates release, not re-scoring, so it does not block D.
Cross-validate the classifier against the official pgscatalog libraries (#237).
Clinical sign-off on the callability floors (§5.2), currently provisional.
Wishlist, not scheduled (2026-09-28): give the measurement-trait scores a place. The PDR2 loss is keyed by disease: each PRS value sits at a position that is a disease node in the graph. Only about 1,620 of the 3,753 eligible scores map to a MONDO disease. The other ~2,033 predict measurements or attributes: LDL/HDL cholesterol, BMI, height, blood pressure, HbA1c, C-reactive protein, kidney function, protein levels, smoking. They have no position, so the model never sees them. That runs against §1.4’s “maximise usable data points”. Options, each a mapping judgement:
- the related disease (LDL → hypercholesterolaemia);
- a phenotype node (HP “increased LDL”);
- the gene or protein node, for protein-level scores, entering beside the expression channel;
- a per-feature-type route in the Phase 12 criterion (docs/pdr2_phase12_per_relation_criterion_planning.md §5).
Not needed for smoke run #1 (#439), which uses the disease-mapped subset.

5.6 What Gate E1 actually measured (2026-09-24)

The first full pass. Ran as three worker tasks joined to one mods_prs_run (6658f656…), all rows written quarantined.

Gate A — the dense callset. 19.1 minutes over the 11,400,292-locus union:

outcome	loci
`CALLABLE_HOMREF`	7,509,057
`CALLABLE_VARIANT`	3,134,326
`EXPLICIT_MISSING`	517,661
`STRUCTURAL_FAIL`	239,248
sum	11,400,292 — exact

The exit invariant held, and the artifact’s sha256 matches its mods_prs_run row.

Gate C — the scores. 2,321 computed, all 2,321 reconciling across the eight outcomes and all 2,321 inside [theoretical_lo, theoretical_hi] — zero violations of the free check that §3.4 calls a floor on correctness.

outcome	rows
`scored`	853,148,795
`explicit_missing`	24,747,219
`ambiguous_palindromic`	13,929,705
`not_in_callset`	1,877,982
`ref_allele_mismatch`	115,560
`unplaceable_row`	0
`missing_other_allele`	1
`unsupported_ploidy`	0

Four of these settle questions that were judgement calls the day before:

The palindrome estimate was right. 13,929,705 deferred rows against the 13,930,267 the manifest predicted — the §5.4 figures stand, and the cost is confirmed as tail-concentrated rather than broad.
unplaceable_row is 0, necessarily: it is an eligibility blocker, so E1 cannot measure what promoting it would cost. Open decision 1 is still open and this run does not touch it. Measuring it needs a deliberate run over the 1,437 sole-blocked scores.
not_in_callset = 1,877,982 is not a defect: the callset emits no row for a STRUCTURAL_FAIL locus, and each of those 239,248 loci is referenced by ~8 scores.
607 of 2,321 scores carry an unverified weight scale — the NR policy (§3.4) now has a number attached to it for Gate F.

The REF-effect fix, quantified. 104,650,679 scored rows sit at a locus where the effect allele is the reference allele, and 29,033,853 of them were resolved from a hom-ref genotype to dosage 2. Every one of those contributed 0 under the variant-only store. That is §5.1 closed in data rather than in code.

n_complement_only_match = 113,313 across 893.8M rows — 0.013 %. Near zero, which is what confirms the forward-strand premise behind deferring strand flips (§5.4) rather than refuting it.

One refusal, and it is the right one. PGS000156 declares weight_type = "Odds Ratio over expected risk", and the scorer refused it: “ratios multiply and must not be summed.” A summed odds ratio would have landed comfortably inside the theoretical bounds and looked entirely plausible — exactly the failure the bounds check cannot catch.

Raw sums span −19,145,748 to +6,182,128. The twelve-order-of-magnitude spread §1.4 predicted from the theoretical bounds, now confirmed on computed values. This is the empirical case that raw sums cannot feed the PDR2 loss, and therefore that Gate D is load-bearing rather than cosmetic.

Operational findings worth carrying into Gate D, all of which cost time here:

The dense callset is 4.55 GB resident, 407 B/locus — 2.5× the 163 B/locus assumed. The --shards guard demands a measured figure for exactly this reason; the old estimate would have authorised 9 concurrent shards on a 30 GB box.
Slicing the work list by contiguous pgs_id put 77.5 % of the variant rows on one shard (2.8 / 19.7 / 77.5 %), because accession correlates with era and era determines size. Wall times were 25.7 / 61.8 / 172.7 min against a possible ~87. Fixed with longest-processing-time-first on n_variants (GH #311), which balances to 1.00×.
The union build (~3 h) and the merge-join must both be covered by _start_pdr_heartbeat_thread. Covering only the phase that looks slow let orphan recovery kill a healthy task three times over.
worker/requirements.txt is a hand-maintained mirror of backend/requirements.txt. pyfaidx had been missing from it since 2026-07-13 and only surfaced after a three-hour union build. scripts/prs/check_worker_env.py now checks the whole path.

5.7 Gate D — the reference panel (staged 2026-09-24; D1–D3 built, not run)

The panel is 3,942 samples, not the 4,150 this section used to claim (verified 2026-09-26 by opening the artifact: .psam has 3,943 lines including its header, and PgenReader.get_raw_sample_ct() agrees). The catalogued HGDP+1kGP release is 4,150; the PGS Catalog’s packaged plink2 build carries 3,942, and the 208 missing are not the relatedness exclusions. Those are a separate list — GRCh38_HGDP+1kGP.king.cutoff.out.id, 612 IIDs — which is not pre-applied, so fitting PCs and the mean/variance regression on unrelated samples leaves 3,330. Super-population split of those 3,330: EAS 735, AFR 684, CSA 678, EUR 667, AMR 409, MID 157. Maureen is essentially entirely European (§5.5), so 667 unrelated EUR samples is the stratum her normalisation leans on hardest even under a continuous-PC model.

Same lesson as the pgsc_all+hm3 retraction below: a documented headline figure is a hypothesis about an artifact, not the artifact. Read the file.

Gate D is not Gate E run four thousand times. This is the single most important design point, because the obvious reading of “a fleet run over thousands of panel genomes” is 3,942 × the E2 run — 8.45 GB of dense callset each and hours of scoring each, so roughly tens of terabytes and ten thousand CPU-hours. That reading is wrong, and it is wrong because a reference panel is not shaped like the patient’s data:

	Maureen	the panel
source	single-sample gVCF with reference blocks	multi-sample matrix, one row per site
the hard problem	telling hom-ref from not-sequenced — the whole point of Gate A	already explicit; every sample’s genotype is on the row
passes needed	one per genome	one, for every sample at once
accumulator	—	`n_samples × n_scores` = 3,942 × 3,753 × 8 B ≈ 118 MB

So the panel pass is one stream over the union carrying a per-(sample, score) accumulator, which factors naturally as a sparse matrix product — dosages (samples × loci) against oriented weights (loci × scores), chunked by locus block so neither operand is ever resident. The per-row Python overhead that set E1’s pace (~28.8k rows/s) disappears into a C-level matmul.

MODS scores the panel, not pgsc_calc. This is the reason the panel is staged locally rather than the whole job delegated to the witness. A percentile means something only if the patient and the panel were deflated the same way, and Gate C deliberately drops missing dosages and defers palindromic SNVs where pgsc_calc mean-imputes to 2p and resolves palindromes by frequency (§5.4). Scoring the panel under a different policy than the patient reintroduces exactly the bias Gate D exists to remove. The witness still runs — on the same panel — and the delta is still a first-class output.

The circular dependency, and how it dissolves. Palindrome resolution needs allele frequencies, which need the panel; but the final scoring pass needs palindromes already resolved. The way out is that the panel yields frequencies before it yields a score distribution, so Gate D splits: D1 stage the panel and extract per-locus allele frequencies over the union → D2 score the panel under the now-final policy → D3 fit PCs and the mean/variance regression. E2 then re-scores Maureen under that same final policy, and only then do z-score and percentile exist.

Panel coverage of the union: 98.93 % at v3 — 90.7 % at v4

Re-measured at v4 scope (2026-09-27), and it is lower. Building the D1 index against the v4 union (read-only, on the real .pvar.zst) found 15,200,494 of 16,755,919 union positions on panel contigs carried by the panel: 90.7 %. The 98.93 % below is a v3 figure. The v4 promotion added ~5.4M positions, and only ~73 % of them are in the panel. This is a count of positions; the row-weighted, per-score figure is what D2 records for every score (outcomes.not_in_panel). Either way, the patient-aligned column below matters more than the v3 number suggested.

Measured against the panel’s own .pvar (2026-09-24): 11,238,703 of 11,359,699 distinct union positions = 98.93 %, uniform at 98.3–99.2 % on every autosome and 99.17 % on chrX. Only 120,996 positions are missing, and only 8 contigs are absent from the panel altogether — alt and random scaffolds carrying 77 positions between them.

(The union has 11,400,292 rows over 11,359,699 distinct positions; the 40,593-row difference is extra allele-sets at a shared position, which the union deliberately keeps separate. Positions is the right denominator when asking what a position-keyed panel carries.)

Superseded, and worth reading as a cautionary note. An earlier revision of this section reported 94.77 %, measured against pgsc_all+hm3_variants_hg38_latest.txt.gz — a 153 MB file sitting in the panel’s own FTP directory, chosen on the strength of its name and a size estimate that seemed to corroborate it. That file is not the panel’s variant list. It holds 33.8M variants; the panel’s .pvar holds 84,325,213. Roughly three quarters of the “missing” loci in that measurement were an artifact of comparing against the wrong file, and so was the per-score analysis built on top of it.

The generalisable part: a filename, a directory location and a plausible size estimate are hypotheses about contents, not contents — and two weak signals agreeing is not verification. What caught it was inspect_reference_panel.py, which reads the actual artifact, on its first run, in one line of output. Check the thing itself before building analysis on it.

The idea behind the cheap pre-check was still right — the packaged panel is a variant subset of unknown extent, and learning that from a small file beats learning it from a 14.9 GB one. The mistake was not verifying which small file answered the question.

Per-score coverage is uniform, and there is no outlier. Measured against the real .pvar (chr22, 2026-09-24):

score	chr22 loci in panel	coverage
`PGS005353`	34,762 / 34,883	99.65 %
`PGS005354`	24,079 / 24,138	99.76 %
`PGS005359`	19,758 / 19,801	99.78 %
`PGS005373`	20,530 / 20,586	99.73 %
`PGS005374`	21,025 / 21,076	99.76 %
`PGS005369`	20,590 / 20,639	99.76 %

Every score sits at ~99.7 %, matching the union-wide 98.93 %. Nothing needs special handling and no score’s reference distribution is meaningfully thinner than any other’s.

Retracted. An earlier revision of this section reported an elaborate finding: that nine scores in the accession band PGS005353–PGS005377 accounted for 96 % of the “missing” loci, that PGS005353 was only 77 % covered, and that Gate D therefore needed a per-score panel-coverage flag because one score’s distribution would rest on three quarters of its loci while looking normal in the manifest.

All of it was an artifact of measuring against the wrong file (see the note above). PGS005353 is 99.65 % covered. The band is real only in the sense that those scores do share a variant list; it has no consequence.

Three hypotheses were falsified along the way — positional clustering, concentration in large genome-wide scores, and %.2f rounding — and a fourth, the one I settled on, was also wrong. The hypothesis never tested was the cheapest one: that the reference I was measuring against was wrong. An unexplained gap is at least as likely to be an instrument error as a finding, and “check the denominator” should come before “explain the discrepancy.”

One genuine lesson survives independent of the artifact, because the scan itself worked: attribution is a search, not an estimate. A 44-of-2,322 sample had ~98 % odds of missing a 9-score cluster, and enumerating the missing items and walking their sources found it in minutes. That technique is sound; it was simply pointed at a phantom.

Is a per-score panel-coverage record still worth keeping? Yes, but as hygiene rather than because the data demands it. Gate C already emits coverage_v1 on the patient side, and a reader should be able to see what a normalized score rests on without re-deriving it. What is not justified is any claim that particular scores need flagging or exclusion — measured coverage is uniform.

Consequence for Gate D, which does not depend on closing that. The per-score number is the one that predicts how well a given score normalizes, and at 100 % for essentially every score sampled it is excellent. The union number governs a different question — which loci enter the reference distribution at all — and the discipline it implies is the same either way: Maureen and the panel must be scored over the same locus set, so the deflation cancels by construction rather than by hope.

D1 and D2, as built (GH #396; not yet run)

backend/mods_prs_panel.py, run as two worker tasks and driven by scripts/prs/run_gate_d_panel.py (--step d1, then --step d2 --shards N, then --step merge). Neither task has run against the real panel: that waits on the worker-image rebuild that pins pgenlib and zstandard (GH #339), and D2 additionally on E3 (below). The one real-data step exercised so far is D1’s index build, read-only on the host (2026-09-27). It streamed all 84,325,213 .pvar records in 221 s at a ~1.4 GB peak and kept 16,055,496 records: 14,802,796 SNV sites, of which 255,021 are split multiallelic, plus 993,947 indel records. It found no multi-ALT records and no REF conflicts, and the position-adjacency guard passed.

D1 (mods_prs_panel_af) md5-checks every panel file against the panel’s own checksums.txt, streams the .pvar.zst keeping the records at union positions (the index D2 scores through), then reads those records from the .pgen in position blocks and counts alleles over the 3,330 unrelated samples, overall and per super-population. The product is union_panel_af.tsv.gz: one line per union target, with af_effect, af_other and an per group, read forward exactly as the union wrote the alleles. That is the panel half of the palindrome comparison. Both alleles are given because at a split multiallelic site one is not 1 − the other. D1 also measures the autosomal missing-call rate over the union’s own sites (the #404 spike saw ~2.4 % on random panel variants), and it records its provenance, including the archive checksum, on a gate_d1_panel_frequencies run row.

D2 (mods_prs_panel_score) resolves every scoring row once into an oriented weight on an allele slot, then streams the panel records it touches in sorted blocks, adding one sparse product per block into the n_scores × n_samples accumulator. Three properties are enforced, not assumed:

Same policy. The policy is her stored Gate C dict, rebuilt with ScoringPolicy.from_dict, which refuses a dict whose fingerprint it cannot reproduce. The E2 rows are such dicts (they predate two policy fields), so a real D2 waits for E3. A golden test also scores every fixture panel sample through score_variants and requires agreement, under the default policy and every per-file variant.
Same missingness rule. A missing panel hardcall is dropped, as her no-calls are, and never imputed. Every matrix row records the mean matched rows dropped per panel sample beside her own count, so “the two deflations match” is a per-score measurement. Smoke run #1 measured that they did not: over the 594 v0 scores her dropped share sat a median 2.3 points above the panel’s (GH #499).
Same locus set. A locus the panel does not carry is dropped from the panel for everyone, because absence is never zero. Her Gate C sum still includes it, so she would be compared on loci the panel never saw. At a REF-effect row that is a systematic +2w on her side only, which is the REF-effect bug’s shape one level up. D2 therefore scores her in the same pass from her Gate A callset, as one more column under the same rules. run_gate_d.py places that panel-aligned score and names it in percentile_v1.score_basis, beside her unchanged gate_c_raw_score. Since gate_d2_panel_scorer_v2 (GH #499) the rule also runs from her side: a locus she has no call at is dropped from every panel sample too. v1 dropped it from her column alone, which moves her z by roughly gap × panel mean / sd, a predicted −5.5 SD on PGS005253. The header says patient_missing_masked; run_gate_d.py refuses a patient column without it unless --allow-unmasked-patient-missingness is passed, which percentile_v1 then records.

The site rules are as follows. Split biallelic SNV records at one position form a site, whose REF count is ploidy − Σ ALT. An allele never observed at a called site counts 0, since it is called absent. An SNV row with no SNV record at its position is not_in_panel, and so is any indel without an exact record match, because indel representation is not canonical. Males outside the chrX PARs, and on chrY, are dropped in D2 (Gate C refuses a non-diploid call) and counted once in D1.

Scoring rows are parsed once. The panel is re-read once per batch of scores that fits --max-batch-rows, at ~12 B per matched row held and ~36 B per row peak while its CSR is built, so batch size trades memory for panel passes. A pass is bound by the Windows bind mount, which the #404 spike measured at ~16.6k records/s in order and ~1k/s scattered, so the real run should stage the .pgen on container-local storage first. Shards split the eligible list LPT on n_variants. --step merge refuses parts that disagree on panel, policy, samples or patient column.

Operational notes

ftp.ebi.ac.uk throttles per connection (~180 KB/s), the same as ftp.ensembl.org; a single stream would take ~30 h for 14.9 GB. Over 8 ranged connections it runs at ~940 KB/s (~4.6 h). More connections is not obviously better and this is a free public mirror.
The archive is md5-pinned against the catalog’s published md5s.txt (2bd3586e…dcee49), and the verified archive is copied to second storage before any cleanup can reclaim it — re-downloading from a throttled mirror is hours. (2026-10-06: that copy is G:\My Drive\maple_archive\reference\panel\pgsc_HGDP+1kGP_v1.tar.zst, with its SHA-256 in the folder’s MANIFEST.tsv. The local archive was deleted to free disk; the extracted panel stays in backend/data/reference/panel/.)
Staging peaks at ~2× the archive, because ranged chunks and the reassembled file coexist. prepare_reference_panel.py refuses to start without that headroom plus a 5 GB floor, and a DiskGuard thread watches free space during the transfer — a pre-flight check is not enough when a locus-union rebuild is writing concurrently.
pgenlib works — read the .pgen directly, no plink2 conversion needed. The earlier note here said to plan on the standalone plink2 binary because pgenlib has no Windows wheel. True on the host, irrelevant to the job: the workers are Linux containers, and import pgenlib succeeds there (also scipy 1.17.1 / scikit-learn 1.9.1, which D3 needs — GH #339). Verified by opening the real artifact, not by importing the module: PgenReader reports 3,942 samples × 84,325,213 variants and read(0, buf) returns dosages. Keeping the read in-process avoids materialising an intermediate dosage dump of a 15 GB panel, and keeps all scoring policy in MODS.

5.7b The road from here — what each remaining stage is for

Approved by Richard 2026-09-26: no release candidate is declared until the palindrome question is settled, and an E3 follows D1. This section is the plain-language version of that plan, including a correction to what “settled” can mean.

The correction, measured 2026-09-26. Resolving a palindromic SNV means comparing the scoring file’s declared effect-allele frequency against the panel’s frequency at that locus: concordant (both ≈ 0.15) means same strand, complementary (0.15 vs 0.85) means flipped. The panel supplies one side of that comparison. The scoring file has to supply the other, and most do not.

across the 1,687 eligible scores with ≥1 palindromic row	scores	palindromic rows	median % of that score’s rows
carry `allelefrequency_effect` — orientable	167 (9.9 %)	6,186,687 (9.5 %)	11.71 %
do not carry it — no anchor exists	1,520	58,647,492 (90.5 %)	10.20 %

So D1 can orient at most 9.5 % of the deferred palindromic rows, and the 90.5 % that stay deferred are not a light tail — their per-score burden is the same size. “Resolve the palindromes” cannot mean all of them; nothing in a reference panel can orient a variant whose file never said what frequency it expected. Treat the residual as a measured coverage property to be gated in F, not as a decision still pending.

The stages, and what each is actually for:

Stage	What it does	Why it exists
D1	One pass over the panel, emitting per-locus allele frequencies across the union. Built (GH #396); not run.	Produces the one number that can orient a palindrome. Nothing here touches scores.
—	Implement: orient palindromes in the 167 EAF-carrying scores; leave the rest deferred and record why. Some of the 167 will have EAF ≈ 0.5 and stay ambiguous — that is expected, not a failure.	This is the policy change. It must be made before D2, because D2 has to commit to a policy.
E3	Re-score the patient under the final policy.	Only the 167 change value, but re-running the fleet keeps one policy label across every row — the alternative leaves a fleet whose provenance has to be argued rather than read (§5.8, GH #360). ~4 h unattended.
D2	Score the 3,330 unrelated panel samples, and the patient beside them, under that same policy. Built (GH #396); not run.	The reference distribution. Same policy or the coverage deflation stops cancelling between patient and panel, and that bias is unrecoverable downstream (§1.4).
D3	Fit ancestry PCs, then regress each score’s mean and variance on them. The regression and percentile are built (PR #386); computing the PCs themselves is not.	D2 gives the cloud of numbers; D3 gives the function that places one person in it. Without it a percentile mostly measures ancestry rather than risk.
F	Weight-scale sanity, per-score coverage gates, then promote.	Where the residual palindrome deferral is finally adjudicated — per score, against a threshold, rather than fleet-wide.

E3 and D2 are independent once the policy is fixed, so they can run in either order or together; D3 needs D2; F needs everything.

Closed 2026-09-26 — GH #371 answered “no”. The open question was whether the 1,520 EAF-less scores could be oriented by some other signal, and the lead was hm_source: if its ENSEMBL value meant “matched to a strand-defined variant record”, ~45 % of deferred rows would need no frequency at all. It does not. hm_source is the provenance of the coordinate, written by the harmonizer’s HmPOS pass, which never reads or writes an allele; strand lives in a separate HmVCF pass whose output the Catalog does not publish, and which labels palindromes ambiguous rather than resolving them. Four lines of evidence, the blunt one being that in PGS005393 hm_source = ENSEMBL labels 26,901 forward-strand rows and 26,801 reverse-complemented ones in the same file.

Two consequences carry forward. The EAF route is the only route, and it delivers 8.48 % of the fleet’s deferred palindromic rows once frequencies within 0.08 of 0.5 are excluded as undecidable — not the 9.5 % headline. And scoring files coded in Illumina’s A/B convention are real and identifiable: their unambiguous SNVs are all A/C or A/G, their strand is per-variant and unrecorded, and a blanket forward-strand assumption is refuted by counterexample rather than merely unproven.

Full answer, the fleet-wide hm_source distribution and the recommended policy: `prs_validation/palindrome_orientation_spike.md`.

But the EAF route is not the only one left, and the spike is what opened the other (GH #404, 2026-09-27). Since harmonization never touches alleles, a file’s non-palindromic rows are a direct readout of its author’s strand convention — and 99.899 % of the 64.8M deferred palindromic rows sit in files that match the reference forward at ≥ 99.9 %, against 8.48 % reachable by declared frequency. Reading those palindromes forward is an inference about that file from thousands of measured rows, not the blanket assumption §5.4 refutes. Its one failure mode — a pipeline that orients unambiguous variants and leaves palindromes alone, which real sumstats harmonisers do — is what the 167 EAF-carrying scores exist to measure. Decide it with D1, before D2.

Measured the same day, without waiting for D1 — the licence holds per class of file, not for the 99.899 %. The spike read the panel directly and tested the 167 EAF-carrying scores. Flipped share: ≤ 0.96 % in small curated files with zero complement-only rows (61 files, 49 publications) and ≤ 0.12 % in genome-wide files. But the failure mode is real: PGS000020 (7,502 rows, 0.175 % complement-only) has ~40 % of its palindromes reversed, so the rate is necessary and not sufficient. The recommendation licenses the small class (527 scores, median 11.3 % of each score’s rows). It puts the genome-wide class (98.8 % of the rows, 87 % of them in three publications with no EAF to test) to Richard as an extrapolation, and refuses 1k–100k. The mechanism, ScoringPolicy.palindrome_forward_scores, ships empty and rides E3: `prs_validation/palindrome_forward_rule.md`.

5.8 What Gate E2 measured (2026-09-26)

The v4-scope fleet run. mods_prs_run 029be6e3, scorer gate_c_scorer_v1, against the Gate A run 64a90788 callset (19,215,298-locus union → 16,964,243 records, sha256 d76e1c65…, verified against the run row and the bytes on disk before scoring). All rows written quarantined. Wall clock 3 h 51 m.

It reconciles at the fleet level, which E1 could not demonstrate. Per-score reconciliation is the eight-outcome assertion the scorer already makes. E2 adds the aggregate:

outcomes summed over the 3,752 computed scores   2,166,261,879
PGS000156, refused, n_variants                               6
                                                 -------------
                                                 2,166,261,885  = manifest eligible rows

Exact. Every scoring row in the eligible set is accounted for by name. reconciles=True on all 3,752 and zero scores outside [theoretical_lo, theoretical_hi].

outcome	rows	share
`scored`	2,036,222,864	94.00 %
`ambiguous_palindromic`	64,810,321	2.99 %
`explicit_missing`	60,837,692	2.81 %
`not_in_callset`	2,738,356	0.126 %
`ref_allele_mismatch`	1,561,651	0.072 %
`unplaceable_row`	90,994	0.0042 %
`missing_other_allele`	1	—
`unsupported_ploidy`	0	—

Five things this settles:

The v4 unplaceable_rows promotion cost is now measured, and it was free. E1 could not measure it — the category was an eligibility blocker, so its outcome was necessarily 0 (§5.6). E2 is the first run with those scores in scope: 90,994 rows, which equals the union manifest’s skipped_no_position exactly, and is 0.0042 % of the fleet. The flat-curve argument behind MAX_UNPLACEABLE_FRACTION = 0.01 holds at full scale. It bought 1,431 scores and ~1.27 B rows for four thousandths of a percent.
The REF-effect fix at full scope. 433,652,126 scored rows sit where the effect allele is the reference allele, and 184,331,602 of them resolved from a hom-ref genotype to dosage 2. Every one contributed 0 under the variant-only store. E1’s figures were 104,650,679 / 29,033,853, so this is 4.1× / 6.3× larger — §5.1 closed in data again, on a bigger sample.
The forward-strand premise still holds. n_complement_only_match = 490,602 = 0.0226 % of rows (E1: 0.013 %). Near zero at 1.9× the scope is what licenses continuing to measure strand flips rather than apply them (§5.4).
The palindrome cost rose, for a known and predicted reason — 2.99 % of rows against E1’s 1.56 %, because v4 admitted the small curated pre-2015 scores where palindromes concentrate. Per-score restatement is in Open decision 2 of §5.5.
Raw sums span −1.91457 × 10⁷ to +6.18213 × 10⁶. Twelve orders of magnitude, now on computed values across the full eligible set rather than on theoretical bounds. This is the empirical case that raw sums cannot feed the PDR2 loss (§1.4) and therefore that Gate D is load-bearing rather than cosmetic.

No blend remains. All 2,321 E1 rows are superseded; zero quarantined rows carry the old callset. The cheap invariant that proves it — and the one that would have caught the hazard earlier — is the gate_a.n_records fingerprint on mods_polygenic_score.metadata:

select run_id, validity, metadata->'gate_a'->>'n_records', count(*)
from mods_polygenic_score group by 1,2,3;

Anything other than one n_records value across the live quarantined rows means a mixed-callset fleet. That is a real hazard, not a hypothetical: the Gate C resume predicate keys on SCORER_VERSION alone and is blind to the callset (GH #360). SCORER_VERSION was unchanged between E1 and E2 while the callset changed completely, so a default resume=True would have skipped every score that already had an E1 row and silently produced exactly that blend — both halves reconciling, both inside bounds, nothing downstream distinguishing them. E2 was run with resume=False over an explicitly computed work list for this reason.

Operational notes worth carrying into Gate D:

The callset measured 8.45 GB resident, 498 B/record — not the 407 B/record E1 measured. Scaling the previous measured figure would have been 22 % low, so B/record is a property of the particular union, not of DenseCallset (§5.3).
Budget concurrency against the Docker VM (30.7 GB), not the Windows host (62.8 GB total / 35.5 GB free at the time). The host figure nearly authorised 4 shards where 2 was right.
Shard count and concurrency are separate knobs, and run_gate_e.py --shards conflates them by design (it assumes N shards = N workers). E2 wanted 20 shards for restart resilience — a killed shard restarts from zero, so small slices bound the loss — and 2-way concurrency for memory. The lever that expresses both is a pipeline_config row (stage_name='mods_prs_compute', enabled=1, max_concurrent=2), which the claim path enforces race-free under pg_advisory_xact_lock. Delete the row afterwards or the task type stays throttled forever. LPT balance over n_variants came out at 1.00×, and all 20 shards completed first try with zero retries.

6. Where the code lives

Concern	Path
Reference FASTA pin + base lookups	`backend/connectors/reference_fasta.py`
PGS Catalog REST/FTP + the one row parser	`backend/connectors/pgs_catalog_connector.py`
Compatibility taxonomy, features, classifier	`backend/mods_prs_compat.py`
Score-aware locus union	`backend/connectors/pgs_locus_union.py`
Dense-callset builder (merge-join)	`backend/connectors/wgs_dense_callset.py`
Scorer (Gate C) — dosage, outcomes, coverage payload	`backend/mods_prs_scorer.py`
Panel (Gate D1/D2) — index, frequencies, panel scoring, patient column, matrix	`backend/mods_prs_panel.py`
Normalization (Gate D3) — PC fit, rank, percentile, `percentile_v1`	`backend/mods_prs_normalization.py`
Score persistence + released-only readers	`backend/mods_prs_service.py`
Read/projection layer + score states	`backend/mods_prs_pipeline_service.py`
REST API	`backend/mods_prs_api.py` (legacy), `backend/mods_prs_pipeline_api.py` (rebuild)
Dashboard tabs	`frontend/src/components/prs/`
Worker tasks	`worker/task_handlers.py` — `genotype_ingest_wgs`, `genotype_dense_callset` (A), `mods_prs_manifest_build` (B), `mods_prs_compute` (C), `mods_prs_panel_af` (D1), `mods_prs_panel_score` (D2)
Gate D drivers	`backend/scripts/prs/run_gate_d_panel.py` (D1, D2 shards, merge), `run_gate_d.py` (D3)
Offline dry run / pgsc_calc reconciliation	`backend/scripts/prs/gate_c_dry_run.py`, `reconcile_pgsc_calc.py`
Scripts (mixed live/deprecated — read the README)	`backend/scripts/prs/`
Tables	`mods_genotype`, `mods_polygenic_score`, `mods_prs_run` (mig. 103), `mods_prs_score_manifest` (mig. 104)

Further reading in this repo: `genomics_glossary.md` · `prs_validation/compatibility_classification.md` · `prs_validation/audit_handoff_2026-07-12.md` · `prs_validation/frontend_requirements.md` · `pdr_personalization_writeup.md` · `essays/whole_catalog_prs.md` (the narrative version of §5.3).

MODS Notebook

Explorer