1 Objectives

In this practical exervise, we will learn how to compute polygenic scores (PGS) for complex traits, or polygenic risk scores (PRS) for complex diseases, using a basic approach: clumping + P-value thresholding (C+PT) based on GWAS results. We will use PRSice to perform C+PT analysis. We will work with GWAS summary statistics from a simulated discovery sample based on real genotypes and individual-level genotype and phenotype data for a target cohort. The exercise is to build PGS from the GWAS summary statistics and predict the phenotypes in the target.

Note: R code is shown with a light blue background, while terminal commands are shown with a light orange background.

2 How were the practical data produced (context only)

The GWAS summary statistics come from a trait simulation built on real genotypes at 273,604 SNPs in 3,000 unrelated individuals of European ancestry in UK Biobank. To simulate a complex trait, 100 SNPs were chosen at random as causal variants, with their effects drawn from a normal distribution, and an environmental residual was added so that the simulated trait had a heritability of 0.5. The simulation was done using GCTA. GWAS summary statistics were then produced with PLINK 1.9 using a linear model with sex, age, and the first 10 principal components (PCs) as covariates.

For validation, we used 494 individuals of European ancestry from the 1000 Genomes Project. Phenotypes were simulated using the same 100 causal variants and effects as in the discovery simulation, so meaningful prediction could be observed even with a small GWAS sample size. We conducted standard GWAS quality control (QC) and removed related individuals in the target cohort.

3 Data location

ls /data/module5/

We prepared the following files for this practical:

Role File or folder name
GWAS summary statistics file gwas.ma
Target genotypes (PLINK prefix) 1000G_phase3.eur.QC.unrel (.fam, .bim, .bed)
Target phenotype target_phenotypes.txt
Covariates (sex, PCs) target_covariates.txt
Tuning / held-out ID lists ind_94.txt, ind_400.txt
PRSice executable PRSice_linux
PLINK executable plink

We have moved the executables to the system-wide path so we can run them directly by typing their names (PRSice_linux, plink, and R).

4 Set up

Copy the data files to your working directory. In the terminal, run:

mkdir PGS
cd PGS
cp -r /data/module5/* .

4.1 Check GWAS summary statistics

For this practical, GWAS summary statistics are provided in GCTA-COJO .ma format: eight columns in a fixed order with the per-SNP fields that PRSice and GCTB expect. Use exactly this order. If you include a header row, the text labels can differ as long as each column’s meaning matches what you declare in your command-line flags.

SNP A1 A2 freq b se p N

Example lines:

SNP A1 A2 freq b se p N
rs3934834 T C 0.146653 0.7677 0.5384 0.1541 2973
rs3737728 A G 0.286691 -0.4979 0.4167 0.2321 2998
rs6687776 T C 0.159326 0.3382 0.5198 0.5154 2997
rs9651273 A G 0.271667 -0.1792 0.4275 0.6751 3000
rs4970405 G A 0.103638 -0.1893 0.6209 0.7605 2996
rs12726255 G A 0.137896 -0.07247 0.5549 0.8961 2995
rs9660710 A C 0.063837 0.08227 0.7632 0.9142 2992

What the columns mean:

  • SNP: usually an rsID. Some public releases use chromosome–position IDs (chr:bp). If your pipeline expects rsIDs (e.g., for SBayesRC), map to rsIDs before matching to LD panels or target genotypes.
  • A1 / A2: A1 is the effect allele (the allele that b refers to). Harmonise with your LD reference and target genotypes. If GWAS A1 is the same physical allele as the LD reference’s A2 (label swap relative to the reference), flip the sign of b. If neither GWAS allele appears in the reference (not just a swap of A1 and A2), the SNP is usually dropped from prediction. Tools such as GCTB will usually perform harmonisation for you, but check outputs when you prepare files by hand or combine cohorts in a meta-analysis.
  • freq: frequency of A1 in the GWAS sample. It can be useful for QC. Large mismatches with the LD reference often lead to exclusion.
  • b / se: marginal effect of A1 and its standard error from the GWAS.
  • p: association P-value.
  • N: GWAS sample size for that SNP (often constant across SNPs but can vary when summary statistics come from a meta-analysis as not every study typed every variant).

4.2 Check target files

PRSice needs:

  • Target genotypes in PLINK 1.9 binary format: prefix.bed, prefix.bim, prefix.fam (you pass prefix to –target). A PLINK2 .pgen setup can also be used; see PRSice file formats for options.
  • Phenotype file with three columns: FID, IID, PHENO (here supplied as target_phenotypes.txt without a header row).
  • Covariate file: FID, IID, then columns such as sex and PCs (used in the PRSice command below via --cov).

Important: the same individuals must appear with the same FID and IID in genotype, phenotype, and covariate files.

Run head on the gwas.ma and target_phenotypes.txt to verify contents.

head gwas.ma
head target_phenotypes.txt

5 What PRSice will do in one pass

PRSice will - Clump SNPs in the target LD structure (default windows and LD \(R^2\) unless you override them). - Build PGS at several P-value cutoffs --bar-levels. - Regress phenotype on each PGS and report \(R^2\), coefficient, SE, number of SNPs, etc.

Default-style clumping parameters (for reference; you do not need to add these unless you want to change them). These are typical PRSice defaults for --bar-levels and clumping:

--bar-levels 0.001,0.05,0.1,0.2,0.3,0.4.0.5,1 \
--clump-kb 250kb \
--clump-p 1.000000 \
--clump-r2 0.100000

Q1: Why can’t we include highly correlated SNPs in the score without clumping?

Show answer

Highly correlated SNPs likely tag the same causal variant. Including highly correlatedSNPs without clumping will upbias polygenic scores.


6 Run PRSice

Run from the directory that contains PRSice_linux. Copy–paste as a block. The line-ending backslashes mean “same command continues on next line”.

PRSice_linux \
--a1 A1 \
--base gwas.ma \
--target 1000G_phase3.eur.QC.unrel \
--pheno target_phenotypes.txt \
--cov target_covariates.txt \
--beta \
--pvalue p \
--stat b \
--bar-levels 1e-8,1e-7,1e-6,1e-5,3e-5,1e-4,3e-4,0.001,0.003,0.01,0.03,0.1,0.3,1 \
--binary-target F \
--fastscore \
--out output
Flag Meaning
--base File of GWAS summary statistics (effect allele and beta must match --a1 and --stat).
--target PLINK prefix for target genotypes (no .bed suffix).
--pheno Phenotype file for the same individuals.
--beta/--stat b Betas are in column b; not odds ratios.
--pvalue p P -values are in column p.
--bar-levels List of P -value ceilings to try (SNPs with GWAS P below each level are candidates after clumping).
--binary-target F Tells PRSice the target phenotype is quantitative (F); required here with –beta (PRSice command reference).
--fastscore Faster scoring option (appropriate for this teaching scale).
--out output All output files start with output.

The covariate file includes sex and 10 PCs (PRSice will adjust for these when you pass --cov):

Check: when the run finishes, your should see new files such as output.prsice, output.summary, output.best, output.log, and possibly output.mismatch.

Verify that your execution directory generates: output.prsice, output.summary, output.best, output.log, and possibly output.mismatch.

7 Read the main outputs

cat output.prsice

One row per P-value threshold tried. Example (abbreviated conceptually from a real run):

Pheno Set Threshold R2 P Coefficient Standard.Error Num_SNP
- Base 1e-08 0.135658 2.91564e-17 18.8386 2.14823 16
- Base 1e-07 0.139878 8.61955e-18 22.4530 2.51520 21
- Base 1e-06 0.150095 4.40527e-19 26.0771 2.80289 25
- Base 1e-05 0.167423 2.62128e-21 31.6102 3.18339 31
- Base 3e-05 0.168322 2.00365e-21 41.9624 4.21231 44
- Base 0.0001 0.148659 6.70352e-19 52.4289 5.66731 63
- Base 0.0003 0.120005 2.54931e-15 77.2140 9.44760 122
- Base 0.001 0.0959434 2.13904e-12 123.3440 17.1120 272
- Base 0.003 0.0460463 1.58737e-06 152.7360 31.4341 648
- Base 0.01 0.0296673 0.000125141 232.3250 60.0858 1797
- Base 0.03 0.0135440 0.00986061 292.3550 112.8470 4864
- Base 0.1 0.00389772 0.167306 299.6170 216.6510 13674
- Base 0.3 9.59294e-06 0.945437 28.8941 421.9800 32840
- Base 1 0.000631923 0.578519 -448.029 805.9280 68751

Q2: Locate the row with the largest \(R^2\). Which Threshold is it? How many Num_SNP are there in that score? Then discuss why \(R^2\) across P-value thresholds often change non-linearly, and what factors might shape that pattern in real data.

Show answer

The best row is 3e-05 with 44 SNPs. \(R^2\) often peaks at an intermediate threshold (a hump-shaped curve): it rises as more signal SNPs enter the score, then falls as too many weak or null SNPs dilute it. You are trading off including more causal variants versus adding noise. In real data, the shape is influenced by polygenic architecture, LD and clumping, signal-to-noise in the SNP set, and sample size.


output.summary

head output.summary

Often one line per analysis configuration, which includes PRS.R2 and related fields. Example:

Phenotype Set Threshold PRS.R2 Full.R2 Null.R2 Prevalence Coefficient Standard.Error
- Base 3e-05 0.17103 0.184156 0.0158338 - 41.9624 4.21231

You may notice that PRS.R2 in output.summary is slightly higher than Full.R2 - Null.R2. For the same threshold, \(R^2\) in output.prsice matches Full.R2 - Null.R2. You will see in Part A4 that incremental R² from R should agree with Full.R - Null.R2, not necessarily with PRS.R2.

output.best

cat output.best

This file shows per-individual PGS from the best P -value threshold. Example:

FID IID In_Regression PRS
HG00097 HG00097 Yes 0.276613636
HG00099 HG00099 Yes 0.142772727
HG00100 HG00100 Yes 0.0610681818
HG00101 HG00101 Yes 0.297465909
HG00102 HG00102 Yes 0.321113636
HG00103 HG00103 Yes 0.258840909
HG00105 HG00105 Yes 0.19825
HG00106 HG00106 Yes 0.0906931818
HG00107 HG00107 Yes 0.0655681818
HG00108 HG00108 Yes 0.207181818

8 Check in R (covariate-adjusted incremental \(R^2\))

Although PRSice provides R², it would be useful to know how to calculate it by hand.

In R, set the working directory to where your files live (or start R in the working directory). The phenotype column is PHENO. The covariate file must list the same FID/IID rows (here: sex and PC1-PC10).

RStudio

# Per-individual PGS from PRSice (best threshold for this run)
pgs <- read.table("output.best", header = TRUE)
pgs$PGS <- pgs$PRS # PRSice names the score column PRS

# Phenotype and covariates for the same target individuals
pheno <- read.table("target_phenotypes.txt", header = FALSE)
names(pheno) <- c("FID", "IID", "PHENO")
covar <- read.table("target_covariates.txt", header = TRUE)

# One row per sample: PGS + phenotype + covariates (individuals have non-missing values in all three files)
d <- merge(merge(pgs[, c("FID", "IID", "PGS")], pheno, by = c("FID", "IID")),
covar, by = c("FID", "IID"))

# covariate names for the formula (sex, PC1-PC)
covar_cols <- setdiff(names(covar), c("FID", "IID"))
null_model <- reformulate(covar_cols, response = "PHENO") # PHENO ~ covariates
full_model <- reformulate(c(covar_cols, "PGS"), response = "PHENO") # PHENO ~ covariates + PGS

fit_null <- lm(null_model, data = d)
fit_full <- lm(full_model, data = d)

r2_null <- summary(fit_null)$r.squared # matches Null.R2 in output.summary
r2_full <- summary(fit_full)$r.squared # matches Full.R2 in output.summary
incr_r2 <- r2_full - r2_null # incremental R²; compare to R2 in output.prsice

c(R2_covariates_only = r2_null, R2_covariates_plus_PGS = r2_full, incremental_R2 =
incr_r2)

Q3: Why use incremental R² (full minus covariate-only model) instead of simply correlating PHENO with PGS in the validation sample (or fitting PHENO ~ PGS with no covariates)?

Show short answers

Covariates such as PCs capture population structure. If you correlate PHENO with PGS without adjusting for them, the correlation can be inflated by between-group differences (ancestry and environment that track ancestry) rather than causal genetic effects on the trait. Incremental \(R^2\) (full model with covariates + PGS minus covariates-only model) measures how much extra variance the PGS explains after those adjustments, so it is fairer and easier to compare across methods.


Q4: We have calculated the prediction accuracy in the validation sample, now think about this question: If you applied the same P -value threshold in a new cohort with the same genetic architecture but new environmental noise, would you expect identical R²? Why or why not?