Polygenic Risk Scores Explained: From GWAS Summary Stats to PRS

A genome-wide association study hands you a table with millions of rows: one variant per line, each with an effect size and a p-value, and almost none of them individually meaningful. A polygenic risk score (PRS) is the trick that turns that sprawling summary-statistics file into a single number per person — an estimate of how much genetic liability an individual carries for a trait. This post explains what that number actually is, why it is just a weighted sum of risk alleles, and how two very different tools, PRSice-2 and LDpred2, compute it in practice.

What a Polygenic Score Actually Is

At its core a PRS is arithmetic, not magic. For each variant you know two things from a GWAS: which allele raises the trait (the effect allele) and by how much (the effect size, beta for continuous traits or log(OR) for binary ones). For a given person you also know their genotype at that variant — a dosage of 0, 1, or 2 copies of the effect allele. Multiply dosage by effect size, do that for every variant, and add it all up:

PRS_i = sum over variants j of ( beta_j * dosage_ij )

That is the entire idea. The score is high when someone carries many copies of many small-effect risk alleles, and low when they don't. No single variant dominates; the signal is spread across thousands or millions of loci, which is exactly why it is called polygenic. Everything sophisticated about PRS methods is about choosing which variants to include and what weight to give each one — because the raw GWAS beta values are noisy and correlated.

Why You Can't Just Sum Every Variant

Two problems make the naive sum a bad idea. First, linkage disequilibrium (LD): nearby variants are inherited together, so a real association at one causal SNP produces dozens of correlated, near-identical signals around it. Summing all of them double-counts the same underlying effect. Second, winner's curse and noise: most variants with a marginal p-value have effect sizes indistinguishable from zero, and including them just adds random noise that dilutes the true signal.

Every PRS method is a different answer to the same question: how do you keep the real signal while discarding the LD redundancy and the noise? The two dominant strategies are clumping and thresholding (prune correlated SNPs, then keep those below a p-value cut) and Bayesian shrinkage (keep everything but shrink each weight toward zero based on how much the data supports it). PRSice-2 is the reference implementation of the first; LDpred2 is the leading example of the second.

Getting Your Inputs in Order

Before any tool runs you need two things aligned. The base dataset is the GWAS summary statistics: a whitespace- or tab-delimited file with columns for SNP ID, effect allele, other allele, effect size, and p-value. The target dataset is your own genotyped cohort — array data, or short reads run through variant calling — held in PLINK .bed/.bim/.fam format or the newer .pgen triple.

The single most common source of wrong scores is allele mismatch. If the base file reports beta for the A allele but your target .bim codes that variant with T as the counted allele, the sign of the contribution flips. Strand ambiguity (A/T and C/G SNPs) makes this worse because you cannot resolve them by flipping alone. Good tools harmonize alleles automatically, but you should still filter to variants with matched IDs, remove ambiguous SNPs when frequency information is unavailable, and confirm your genome build matches between base and target.

Clumping and Thresholding with PRSice-2

PRSice-2 implements the classic C+T (clumping + thresholding) pipeline and, crucially, automates the search for the best p-value threshold. Clumping is the LD step: PLINK's --clump algorithm walks the summary stats from most to least significant, keeps the lead SNP, and removes every nearby variant in LD with it above an r^2 cutoff (commonly 0.1) within a window (commonly 250 kb). What survives is a roughly independent set of index SNPs.

Then comes thresholding. Rather than guess one p-value cut, PRSice-2 builds many scores — one at each of a grid of thresholds like 5e-8, 1e-5, 0.05, 0.1, up to 1.0 — and tests which threshold's score best predicts the phenotype in the target data. A typical run looks like this:

Rscript PRSice.R --dir . \
  --prsice ./PRSice_linux \
  --base gwas_sumstats.txt \
  --target target_cohort \
  --snp SNP --A1 A1 --A2 A2 --stat BETA --pvalue P \
  --clump-kb 250 --clump-r2 0.1 \
  --bar-levels 5e-8,1e-5,0.05,0.1,0.5,1 \
  --fastscore \
  --binary-target T \
  --out prsice_run

PRSice-2 outputs the best-performing threshold, the variance explained (R^2), and per-individual scores. One caveat worth internalizing: because it optimizes the threshold against the target phenotype, you are fitting a parameter to your test data. To report an honest predictive R^2 you must evaluate the final score in a held-out sample the threshold was never tuned on, or the number is optimistically biased.

Bayesian Shrinkage with LDpred2

LDpred2, part of the bigsnpr R package, takes the opposite philosophy: don't throw variants away, reweight all of them. It models each SNP's true effect as drawn from a prior distribution and updates every beta using a reference LD matrix — so instead of clumping away correlated SNPs, it explicitly accounts for their correlation and redistributes the effect across them. The result is a set of posterior effect sizes that are then summed exactly like any other PRS.

LDpred2 comes in a few flavors. LDpred2-inf assumes an infinitesimal model where every variant is causal with a tiny effect. LDpred2-grid searches over a grid of two hyperparameters — the proportion of causal variants (p) and the SNP heritability (h2). LDpred2-auto estimates both directly from the data with no validation set required, which makes it especially attractive when you can't spare samples for tuning. In practice you build the LD reference from a subset of your data or a matched panel, run the model, and get posterior weights.

Applying the Weights with PLINK2 --score

Whichever method you use, the final step is identical: take a three-column file of variant ID, effect allele, and (possibly reweighted) effect size, and apply it to genotypes. PLINK2 --score is the workhorse for this and handles allele matching, missing-genotype imputation to the mean, and dosage from imputed data:

plink2 --pfile target_cohort \
  --score weights.txt 1 2 3 header cols=+scoresums \
  --out final_prs

The three numbers 1 2 3 name the columns holding the variant ID, effect allele, and weight. By default PLINK2 reports the average per-allele score (dividing by the number of scored alleles), which correctly handles missing data; cols=+scoresums adds the raw sum. Both PRSice-2 and LDpred2 ultimately produce weight files you could feed straight into this command — the methods differ only in how they compute column 3.

Choosing a Method

There is no universally best PRS method; the right choice depends on data, ancestry, and how much tuning you can afford.

As a rule of thumb, LDpred2 and other genome-wide Bayesian methods tend to outperform C+T for highly polygenic traits with large training GWAS, while C+T remains a fast, transparent, and surprisingly competitive baseline — especially when the true genetic architecture is sparse. Running both and comparing held-out R^2 is entirely reasonable and often expected.

Evaluating and Interpreting the Score

A PRS is only as good as its out-of-sample performance. Fit a regression of phenotype on the score plus covariates — age, sex, and the top genetic principal components to control for population structure — and report the incremental R^2 (for continuous traits) or the area under the ROC curve (for binary ones) that the score adds beyond covariates alone. The genetic principal components matter: without them, a score can look predictive simply because it tracks ancestry rather than biology, a confounding trap that also haunts upstream genotype quality control.

Interpretation demands humility. A PRS gives a probabilistic, population-relative estimate of liability, not a diagnosis. Someone in the top percentile of a coronary-artery-disease score has elevated average risk, but many such people never develop the disease and many in the bottom percentile do, because environment and rare variants matter too. Treat the score as one predictor among many, and always state the training GWAS's ancestry.

Practical Takeaways

• A PRS is a weighted sum of risk-allele dosages: sum(beta * dosage) over many variants, nothing more exotic than that.
• The hard part is not the sum but choosing weights — LD and noise mean you cannot use raw GWAS betas directly.
• Use PRSice-2 for a fast, interpretable clumping-and-thresholding baseline; use LDpred2 for genome-wide Bayesian shrinkage on highly polygenic traits.
• Harmonize alleles between base and target obsessively — a flipped effect allele silently corrupts the whole score.
• Always evaluate in a held-out sample with genetic principal components as covariates, and report percentile, not raw score.
• A European-trained PRS transfers poorly across ancestries; state the training population and never assume portability.

If you're setting up the tooling, both workflows lean on PLINK and R — start by getting comfortable with the PLINK genotype file formats that every step reads and writes, and make sure your R environment is reproducible with conda environments and Bioconda channels before installing bigsnpr.