bpred

Bayesian polygenic score (PGS) training from GWAS summary statistics, with automatic inference of the genetic architecture.

bpred fits an additive effect-size model to GWAS summary statistics and a reference LD panel using a Gibbs sampler, jointly inferring:

• h² — SNP heritability
• p — polygenicity (the fraction of causal variants)
• α — the MAF-dependent effect-size scaling parameter

The model is closely related to LDpred2-auto (Privé et al., 2023, "Inferring disease architecture and predictive ability with LDpred2-auto"). Beyond the single-trait model, bpred includes two-trait extensions that borrow strength across genetically correlated traits:

• a shared model (a common causal component plus trait-specific components), and
• a directed model (one trait's effects feed into another).

Models

All three models place a spike-and-slab prior on per-allele effects and infer the posterior with a multi-chain Gibbs sampler (default 30 chains, 500 burn-in + 500 sampling). LD enters as a pruned, windowed block-diagonal correlation matrix, so the per-SNP scan parallelizes over (chain × block). The slab variance is MAF-dependent: a variant with reference allele frequency f has per-allele SD S = sqrt(2·f·(1−f)), and a causal effect is drawn from N(0, σ²·S^(2α)). The exponent α controls how effect size scales with allele frequency (α = −1 ≈ effects independent of MAF on the standardized scale).

Hyperparameters are, by default, point-estimated each iteration — α, σ² by L-BFGS-B MLE, p by its posterior mean — rather than sampled, making the sampler a stochastic-EM hybrid. Setting sample_hyperparams=True instead draws them from their full conditionals (a Laplace/Gaussian draw for (α, σ²), a Beta draw for p). An optional QC step filters variants whose GWAS-implied SD is inconsistent with the reference panel (the LDpred2-auto consistency checks).

Single-trait (bpred.py)

The base model, closely related to LDpred2-auto. Each SNP is causal with probability p; causal effects follow the MAF-dependent slab above. The sampler jointly infers h² (SNP heritability), p (polygenicity), and α, and returns posterior-mean PGS weights, per-SNP posterior causal probabilities, credible intervals for each hyperparameter, and the implied PGS variance.

Shared two-trait (bpred_shared.py)

A symmetric model for two genetically correlated traits. Each SNP carries three independent effect components, each with its own spike-and-slab (p, α, σ²):

• γ_s — a shared effect that acts on both traits,
• γ₁ — a trait-1-specific effect, and
• γ₂ — a trait-2-specific effect.

The observed marginal effects are β̂₁ = S⁻¹RS(γ_s + γ₁) + noise₁ and β̂₂ = S⁻¹RS(λ·γ_s + γ₂) + noise₂, where λ rescales the shared component into trait 2 (the trait-1 loading is fixed to 1 for identifiability). Because each component is independently spike-and-slab, a SNP can take any of 8 causal configurations. A second cross-trait parameter ρ captures the correlation between the two GWAS error terms induced by overlapping samples. The shared component lets the two traits borrow strength from each other where their genetic architectures overlap.

Directed two-trait (bpred_directed.py)

An asymmetric model in which trait 1's effects feed into trait 2 but not the reverse (a causal / Mendelian-randomization-style structure). With cross-trait matrix [[1, 0], [λ, 1]], the latent per-trait effects relate to the observed ones as β₁ = γ₁ and β₂ = λ·γ₁ + γ₂. Each trait keeps its own h², p, α, σ², and causal indicators; λ is the directed effect of trait 1 on trait 2 and ρ again absorbs sample-overlap error correlation.

Repository layout

bpred/
  core/         model + data plumbing
    bpred.py            single-trait Gibbs sampler
    bpred_shared.py     two-trait shared model
    bpred_directed.py   two-trait directed model
    ldref.py            loads the LD panel (snps.parquet + blocks/)
    gwas.py             GWAS harmonization onto the panel (shared helpers)
    block_ld.py         windowed (banded) block-diagonal LD container
    config.py           reads config.json
  gwas/         one loader module per GWAS (source, column map, load())
    iq_savage_jansen_2018.py
    ea_okbay_2022.py
  external/     wrappers for baseline tools
    ldpred2.py, run_ldpred2.R   LDpred2-auto baseline
    mtag.py                     MTAG baseline
  simulation/   GWAS simulation under the model
  scripts/      runnable entry points (see Usage)
  utils/        LD-panel construction utilities
  tests/        unit tests
config.example.json   template for your local config.json

Installation

Python 3.10+. Install the dependencies into a virtual environment:

python3 -m venv .venv
source .venv/bin/activate
pip install -r requirements.txt

Some utilities additionally require R. Install these packages (any recent version; Matrix already ships with R):

Rscript -e 'install.packages(c("jsonlite", "bigsnpr", "bigreadr", "data.table"))'

• the LD-panel build utilities (bpred/utils/) use Matrix and jsonlite;
• the LDpred2 baseline (bpred/external/) uses bigsnpr, bigreadr, data.table, and Matrix.

The R scripts do not self-install: if a package is missing they exit with the exact install.packages(...) command to run. Installing bigsnpr from source needs a C++/Fortran toolchain (gcc, gfortran, make).

The MTAG baseline requires a separate MTAG install (see External baselines). These are only needed if you run those specific comparisons/utilities.

Most entry points resolve config.json relative to the repository, so run commands from the repository root.

Configuration

Copy the template and edit it for your machine:

cp config.example.json config.json

config.json is git-ignored. Any key you omit falls back to the defaults in bpred/core/config.py.

Key	Meaning	Default
working_dir	Harmonized-GWAS caches and intermediate files	working_dir
ld_dir	Base directory for LD data (see below)	working_dir/ld
gwas_dir	Raw GWAS summary statistics, laid out as {gwas_dir}/{file}	working_dir/gwases
mtag_path	Path to mtag.py (only for the MTAG baseline)	mtag/mtag.py
mtag_python	Interpreter used to run MTAG (e.g. its conda env's python)	python

ld_dir holds two subdirectories: {ld_dir}/bpred_ldref is the built bpred panel (snps.parquet + blocks/) the models consume, and {ld_dir}/ldref_hm3_plus is the raw LDpred2 HapMap3+ reference (panel-build input and LDpred2 baseline).

Data setup

1. Get the LD reference panel

bpred consumes a pruned, windowed block-diagonal LD panel (bpred_ldref) built from the published LDpred2 HapMap3+ reference.

Recommended — download the prebuilt panel (~4.8 GB) from Hugging Face into {ld_dir}/bpred_ldref. The fetch is resumable and sha256-verified:

python -m bpred.utils.download_ldref

(Source: https://huggingface.co/datasets/bkorpan/bpred_ldref, CC BY 4.0.)

Alternative — reproduce the panel from scratch. Only needed to rebuild it from the raw LDpred2 reference. Step (b) runs in R, so install the R dependencies first (see Installation):

# (a) Download the raw LDpred2 hm3_plus reference (~14 GB) into {ld_dir}/ldref_hm3_plus
python -m bpred.utils.download_ldpred2_ldref

# (b) Prune + window per chromosome (writes transient band files)
Rscript bpred/utils/build_bpred_ldref.R

# (c) Finalize into snps.parquet + blocks/chr{c}.npz under {ld_dir}/bpred_ldref
python -m bpred.utils.build_bpred_ldref

Steps (a)–(c) all derive their paths from ld_dir in config.json, so they stay consistent: (a) downloads to {ld_dir}/ldref_hm3_plus, (b) reads that and writes the built panel to {ld_dir}/bpred_ldref, and (c) finalizes that same panel. The result is byte-identical to the prebuilt download above.

2. GWAS summary statistics

Each GWAS has a loader module under bpred/gwas/ that records its source and column map and exposes load(). The raw summary-statistics file lives under gwas_dir; on first load() it is harmonized onto the panel and cached as parquet under working_dir/harmonized_gwas/.

The two datasets used by the training scripts:

• IQ — Savage-Jansen 2018 (bpred.gwas.iq_savage_jansen_2018). Publicly downloadable; load() fetches it automatically on first use. To pre-fetch it explicitly:

  python -m bpred.gwas.iq_savage_jansen_2018   # downloads + extracts into gwas_dir

• EA — Okbay 2022 / EA4 (bpred.gwas.ea_okbay_2022). Requires an SSGAC data-access account, so it can't be fetched automatically — request the EA4 results (additive, excluding 23andMe) from the SSGAC and place the file at {gwas_dir}/EA4_additive_excl_23andMe.txt.gz. See the module docstring for details.

To add another GWAS, copy one of these modules: set its source/filename, COLUMN_MAP (canonical fields → file columns), SEPARATOR, and any per-dataset preprocessing, then call gwas.harmonized_gwas(...) from load().

Usage

All entry points are Python modules run from the repository root:

# Single-trait training on real GWAS, with held-out validation R²
python -m bpred.scripts.train_bpred

# Two-trait training
python -m bpred.scripts.train_shared       # shared model (e.g. IQ + EA)
python -m bpred.scripts.train_directed     # directed model (e.g. IQ -> EA)

# Simulation studies (validate against known causal effects)
python -m bpred.scripts.simulate_bpred
python -m bpred.scripts.simulate_shared
python -m bpred.scripts.simulate_directed

# Comparisons
python -m bpred.scripts.compare_with_ldpred2      # bpred vs LDpred2-auto
python -m bpred.scripts.compare_two_trait         # shared vs directed vs MTAG

Each training function defaults to the GWAS loader modules in bpred/gwas/ (e.g. train_bpred(training_gwas=iq_savage_jansen_2018, validation_gwas=None)); pass different loader modules to use other datasets.

External baselines

• LDpred2-auto — wrapped by bpred/external/ldpred2.py, which calls run_ldpred2.R. The R script reads the LD reference directory from the BPRED_LDREF_HM3_PLUS_DIR environment variable; the Python wrapper sets this automatically from config.ldpred2_dir() ({ld_dir}/ldref_hm3_plus). To run the R script directly, export that variable yourself.

• MTAG — wrapped by bpred/external/mtag.py. MTAG is a Python 2.7 tool, so it needs its own interpreter. Create the pinned conda environment, clone MTAG, and point config.json at them (conda required):

  conda env create -f bpred/external/mtag-environment.yml
  git clone https://github.com/JonJala/mtag.git mtag

  Then set in config.json:

  "mtag_path":   "mtag/mtag.py",
  "mtag_python": "<conda envs dir>/mtag/bin/python"

  (conda run -n mtag which python prints the interpreter path.) Any working MTAG install will do — the env file just pins a combination known to run it.

Tests

python -m pytest bpred/tests

License

MIT — see LICENSE.