Build SVD embeddings from PheWeb GWAS summary statistics and validate that the implementation choices (p-value threshold, component count, scaling) are defensible. Configured for MGI-BioVU (69 laboratory-test GWAS).

The goal is not just to produce embeddings but to show, with external evidence, that the chosen approach yields embeddings that capture real shared genetic architecture and beats the alternatives — the prerequisite for publishing and for the downstream variant-attribution work.

config.py                    shared config (dataset, p-value threshold, paths)
requirements.txt
pipeline/                    build the data + embeddings — run in order
  1a_download_mgi_biovu_phenotypes.py
  2_build_matrix.py
  3_create_embeddings.py
  4_analyze_similarity.py
validation/                  is the embedding good, and how was the threshold chosen?
  lab_pairs.py               ground-truth lab pairs + literature keys (reference data)
  5_reconstruction.py        reconstruction diagnostic (Pearson + Spearman vs k)
  6_threshold_sensitivity.py p-value threshold: coverage + known-biology basis for 5e-5
  7_component_coherence.py  DeGAs-style validation: components map to physiological panels
docs/                        write-ups (see "Documentation" below)
  archive/                   superseded deliverables, kept for reference
data/{DATASET}/              all outputs + caches (git-ignored)

Always run scripts from the repo root (not from inside pipeline/ or validation/), e.g. python validation/7_component_coherence.py. Each script adds the root to its import path so config.py is found.

Run in order. Set the dataset and parameters in config.py first.

python pipeline/1a_download_mgi_biovu_phenotypes.py
python pipeline/2_build_matrix.py
python pipeline/3_create_embeddings.py
python pipeline/4_analyze_similarity.py

Note: steps 3–4 plot only the first 2–3 components, which capture ~21% of the signal (see validation/5_reconstruction.py). Cosine similarity on the full embeddings is the right tool; the 2D plots are for intuition only.

python validation/5_reconstruction.py
python validation/6_threshold_sensitivity.py
python validation/7_component_coherence.py

How the threshold was chosen and the embedding validated is documented end-to-end in docs/THRESHOLD_AND_VALIDATION.md. In brief: no metric computed on the matrix can select the p-value threshold (it is a data-quantity vs signal-quality trade-off, so in-sample and held-out reconstruction favor opposite extremes). Following matrix-decomposition precedent (DeGAs used p<10⁻³), we use a permissive cutoff (5×10⁻⁵) justified by coverage, and validate the embedding by the biological coherence of its components.

• Embedding captures real, physiologically-organized genetic structure. DeGAs-style component validation: the latent components map onto coherent physiological panels (lipid, hepatic, renal, erythroid, leukocyte, iron …), with variance-weighted panel purity ~0.16 above a permuted-label null. This is the main validation.
• Reconstruction is faithful. ~90% at ~25 components for Pearson and Spearman alike, so the matrix is genuinely low-rank and not an artifact of a few large-beta outliers (diagnostic, validation/5_reconstruction.py).
• p-value threshold: 5×10⁻⁵ (decided). No matrix-internal metric can select it — it's a data-quantity vs signal-quality trade-off. The choice rests on: method precedent (a permissive inclusion threshold, more conservative than DeGAs's 10⁻³ — not a named GWAS convention), coverage (only 1/69 phenotypes starved vs 7/11/18/20 at stricter cutoffs), and biological coherence of the components. Honest caveat: finer similarity structure is threshold-dependent, so interpretation rests on the strongest relationships. Full story + paper-ready text in docs/THRESHOLD_AND_VALIDATION.md.
• Embeddings used throughout are Σ-scaled phenotype vectors (V·Σ), ~20–30 components.

Some known-correlated pairs are recovered weakly (ALT/AST, CRP/ESR, amylase/lipase, CK/CK-MB): expected, because these are clinically co-elevated through tissue damage but have weak common-variant genetic correlation, and the embedding only sees genetics. See docs/CITATIONS.md.

validation/lab_pairs.py holds 24 positive pairs (documented shared loci, positive direction) and 10 negative-control pairs, each with a rationale and a literature reference key. Inverse relationships (TSH/FT4, TIBC/iron) are excluded because signed cosine would point them apart. docs/CITATIONS.md gives the per-pair evidence and full references — 16 of the 24 positives are covered by two consortium blood-cell GWAS (Astle 2016, Vuckovic 2020) — and flags the four weak-genetic pairs kept as an honest stress test. Verify PMIDs before submission.

• DATASET_NAME, PHEWEB_BASE — switch PheWeb instance (UKB-TOPMed / MGI presets included).
• PVAL_THRESHOLD (default 5e-5) — significance filter; also used by get_latest_matrix() to pick the canonical matrix to validate on (see docs/THRESHOLD_AND_VALIDATION.md).
• N_COMPONENTS (50) — components for the exploratory step 3 only; the validation scripts use ~20–30.
• VARIANT_SAMPLE_SIZE — subsampling for variant scatter plots.
• PHENO_DATA_PATH is legacy/unused (step 1 writes per-phenotype parquet instead).

• raw/ — downloaded summary stats (cache); positions.parquet — rsID → chrom/pos.
• beta_full.parquet, pval_full.parquet — aligned matrices (threshold-independent cache).
• beta_matrix_p*.csv — beta matrix at a given threshold.
• svd_validation/ — reconstruction-diagnostic metrics (reconstruction_metrics.csv), figures (reconstruction_vs_k.png, significant_scatter_k*.png), and SUMMARY.txt.
• method3/threshold_sensitivity.{png,csv} — coverage + known-biology basis for 5e-5; method3/component_coherence.{png,csv,_SUMMARY.txt} — DeGAs-style component validation.

• docs/THRESHOLD_AND_VALIDATION.md — start here for the threshold + validation story: full procedure (background → why selection fails → DeGAs precedent → decision → validation → conclusions) with paper-ready text.
• docs/CITATIONS.md — literature support for the lab pairs.
• docs/archive/ — superseded/historical material kept for reference: the earlier meeting memo and Korean progress summary (both predate and partly describe since-removed approaches), and the May-2025 Word report.

Optionally adopt the DeGAs construction more fully — a z-score fill (z = sign(β)·Φ⁻¹(1−p/2)) with a standard-error filter — which removes the hard threshold entirely and settles the NaN→0 question; replicate on FinnGen/UKB (check allele coding and lab-trait availability first — FinnGen is mostly disease endpoints). The research contribution beyond validation is variant attribution (which variants drive the signal; near-significant variants → gene → rare-disease paths), where LD handling becomes essential.

pip install -r requirements.txt

(requests, pandas, scikit-learn, matplotlib, numpy, scipy, seaborn.)