PheWeb Beta Matrix Builder
A reusable, config-driven pipeline that ingests any PheWeb instance and outputs a sign-consistent variant × phenotype β-matrix and its low-dimensional embeddings — the data foundation for downstream multi-phenotype analysis.
Overview
Most GWAS analyses look at one trait at a time. This pipeline builds the infrastructure for looking at all of them together — turning per-trait summary statistics into a single sign-consistent variant × phenotype matrix and ready-to-use embeddings.
The Genome-Wide Association Study (GWAS) field has spent two decades publishing per-trait summary statistics, and several large research consortia (MGI-BioVU, UKB-TOPMed, FinnGen) now publish them through PheWeb — a public viewer that exposes the data through a standard API. The naive way to use that data is one trait at a time; the more interesting question is what happens when you align many traits into a single matrix and look for shared structure. That's a non-trivial engineering problem: each PheWeb instance has dozens to thousands of phenotypes, allele orientations aren't guaranteed consistent across traits, palindromic SNPs are ambiguous, and the resulting matrix needs careful preprocessing before SVD. The pipeline handles all of it as four separable stages — parallel downloads with schema detection, sign-consistent allele alignment, truncated SVD into a shared embedding space, and similarity analysis — so the matrix-construction work is done once, carefully, and the downstream research doesn't have to redo it.
Stack
The Problem
Public PheWAS data is everywhere. A reusable, sign-correct preprocessing pipeline for getting it into a form where you can study cross-phenotype structure isn't.
Every group that wants to use PheWAS data for cross-phenotype analysis ends up writing the same one-off script: download N phenotypes, intersect rsIDs, align alleles, build the matrix. The same fragile preprocessing gets rewritten across labs, the same edge cases (palindromic SNPs, duplicate rsIDs, missing p-values, gzip vs. plain TSV, column-naming inconsistencies) get rediscovered, and the resulting matrices aren't directly comparable between studies. Two specific failures recur: (1) sign errors from inconsistent allele orientation silently corrupt every downstream analysis; (2) inconsistent schemas between PheWeb deployments mean code written for one biobank breaks on the next.
Methodology
Four separable stages, each tunable independently. The hardest one is stage 2 — allele orientation. Getting it wrong silently corrupts every downstream analysis.
| Stage | What it does | Why it matters |
|---|---|---|
| 1 · Download | Parallel per-phenotype download via PheWeb's /api/phenotypes.json. Schema auto-detection handles column-name variants (rsids/rsid, pval/p_value/p). Resumable, column-pruned for speed. |
The same PheWeb data is published with slightly different schemas across deployments. Detecting columns at runtime is what makes the pipeline portable across biobanks. |
| 2 · Align & Build | Intersect rsIDs across all phenotypes. For each phenotype beyond the first, compare ref/alt allele orientation against the reference: match → keep β · flip → negate β · mismatch → NaN. Apply p-value threshold to retain only significant associations. | GWAS effect sizes depend on which allele is called "reference." If two files disagree on ref/alt, the sign of β must be flipped for values to be comparable. This is the step that silently corrupts every downstream result if done wrong. |
| 3 · Embed | TruncatedSVD on the assembled β-matrix → variant embeddings (U·Σ, rows = rsIDs) and phenotype embeddings (V, rows = phenotypes). Saved as CSV alongside 2D/3D PCA scatter and variance-explained plots. | Project a sparse, ~thousand-dimensional phenotype space and a ~million-row variant space into a shared low-dimensional latent space where downstream similarity work becomes possible. |
| 4 · Analyze | Cosine similarity across the embeddings: phenotype ↔ phenotype, variant ↔ variant, and variant ↔ phenotype cross-similarity. Distributions + top-pair tables. | The point of constructing the matrix is to ask cross-phenotype questions. Stage 4 is where those questions become queryable; everything before is plumbing. |
Why each case in stage 2 is named explicitly. The first version of the alignment logic collapsed match / flip / mismatch into a single transformation and produced sign errors in ~5% of variants. Rewriting with explicit, named cases made the failure modes legible and the bug went away. Palindromic SNPs (A/T, C/G — strand-ambiguous) currently fall into the mismatch → NaN bucket; the cleaner long-term fix is allele-frequency matching (the MungeSumstats / LDpred2 convention).
Pipeline
Four numbered scripts, run in order. Threshold-independent intermediate cache (beta_full.parquet, pval_full.parquet) means re-running at a new significance threshold is fast — no re-download, no re-alignment.
Stage 2 memory discipline. The matrix can be large (variants on the order of 10⁵–10⁶, phenotypes in the 10²–10³ range depending on biobank). Stage 2 streams one phenotype at a time during alignment to keep peak memory bounded, then writes the assembled β / p-value matrices to parquet. Stage 3 reads the cached matrices, applies the configured threshold, and runs SVD.
Switching biobanks is a one-line config change. The pipeline currently has presets for MGI-BioVU and UKB-TOPMed; pointing it at FinnGen requires the same one-line dataset switch in config.py (with a caveat that FinnGen is mostly disease endpoints, so the analysis questions you'd ask of it are different).
Findings (MGI-BioVU · 69 laboratory-test GWAS)
The constructed matrix passes the basic sanity checks for a defensible embedding space — reconstruction holds under both Pearson and Spearman, the SVD captures real shared structure, and the embeddings cluster known-related phenotypes together.
| Sanity check | Result |
|---|---|
| Reconstruction (top-k SVD vs. original β-matrix) | ~90% Pearson and Spearman correlation at ~25 components — not a few-outliers artifact. Spearman exceeds Pearson at low k, meaning the embedding recovers rank ordering of effects before it recovers magnitudes. |
| Σ-scaling on the phenotype embedding | Required. Unscaled V collapses at high k because noise components get weighted equally with signal. Σ-scaled stays robust. |
| Embedding structure captures genetic similarity | Validated against an external ground truth of 24 known-correlated lab pairs + 10 negative controls. Known pairs cluster at the high end of the cosine-similarity distribution. |
Validation (brief)
Matrix construction is the work this project is about. Validation is the harness that proves it's research-grade — there's a separate validation/ half of the repo organized as three methods with a hierarchy.
- Method 1 — reconstruction (diagnostic). Confirms SVD is implemented correctly. Cannot rank choices (reconstruction is monotonic in k).
- Method 2 — LD pruning vs. baseline (marker). Tests whether SVD already absorbs LD redundancy. Re-scored with Method 3.
- Method 3 — known lab-pair separation (decision metric). Independent labels (24 positive + 10 negative-control lab pairs, each with literature citation) → AUROC + Wilcoxon. This is the metric that legitimately ranks choices like p-value threshold, k, LD handling, and scaling.
The full validation framework — including the threshold sweep, AUROC findings, and the open question about whether to replace threshold-then-zero with z-score filling — lives in docs/THRESHOLD_ANALYSIS.md and docs/CITATIONS.md. Linked in the buttons below if you want the full story.
Limitations
- NaN → 0 imputation is implicit in SVD. Current zero-imputation conflates "no signal" with "p ≥ threshold," which silently biases the embedding. Replacing with softImpute or iterative SVD-with-missingness is in the open work.
- Palindromic SNPs are dropped to NA when flipped. A/T and C/G strand-ambiguous variants need allele-frequency matching (MungeSumstats / LDpred2 convention) or upfront exclusion. Currently dropped.
- Single-biobank deployment so far. Most extensively run on MGI-BioVU. UKB-TOPMed preset exists; FinnGen would need additional handling because it's mostly disease endpoints, not lab tests.
- The downstream embedding analyses are observational. The pipeline produces the matrix; the matrix is a substrate for hypothesis generation, not a substitute for causal-genetics designs (MR, fine-mapping).
Lessons Learned
The two unglamorous decisions — how to handle allele orientation and how to handle missing data — decide whether the matrix is research-grade or research-garbage. Most of the engineering effort went into them, and most of the open work still left is in them.
What worked. Writing the allele-orientation logic with explicit, named cases (match / flip / mismatch → NA) rather than collapsing them into a single transformation. The first version conflated cases and produced sign errors in ~5% of variants before I rewrote it. Naming each case explicitly made the failure modes legible.
What also worked. The threshold-independent intermediate cache (beta_full.parquet, pval_full.parquet) — caching the aligned matrices once and reapplying thresholds later turned threshold experiments from hours-long re-runs into seconds-long queries.
What didn't. The current p-value threshold (5e-5) is hard-coded as a default rather than tested against a sensitivity sweep from day one. I knew it mattered — different thresholds produce different embedding structures — but I didn't yet have the apparatus to compare embeddings across thresholds rigorously. The validation framework (Method 3) is what fixed this later.
What I'd do differently. Build the sensitivity-analysis harness alongside the core construction pipeline, not after. Specifically, parameterize the threshold and add a comparison step that quantifies how much the downstream embedding cosine-similarity ranking changes when the threshold moves. Same principle as building model evaluation before model training.