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 work with two simulated data sets:

  • A small data set with 10 SNPs, suitable for manual calculations and easy manipulation in R.
  • A large data set with ~300k SNPs across 22 chromosomes simulated from UK Biobank, representing a more realistic scenario.

We will use the small data set to illustrate the principles and calculate PGS by hand in R. Then, the larger data set will be analysed with a pipeline using PLINK (https://www.cog-genomics.org/plink/1.9/).

Note: the small data set can be freely downloaded. The UK Biobank data set cannot be downloaded because it contains real individual genotypes, although individual IDs are masked.

2 Data

2.1 Small data set

This toy data set includes:

  • A training population (discovery GWAS) of 325 individuals, genotyped for 10 SNPs.
  • A validation population of 31 individuals (hold-out sample).

The trait was simulated such that SNP 1 has an effect size of 2 and SNP 5 has an effect size of 1. The trait heritability is set at 0.5.

Data location:

ls /data/module5/toy/

Example R code to load the data:

# data for training population
data_path="/data/module5/toy/"
X <- as.matrix(read.table(paste0(data_path,"xmat_trn.txt")))
y <- as.matrix(read.table(paste0(data_path,"yvec_trn.txt")))

# data for validation population
Xval <- as.matrix(read.table(paste0(data_path,"xmat_val.txt")))
yval <- as.matrix(read.table(paste0(data_path,"yvec_val.txt")))

You can download this small data set using the following Linux command in the Terminal (or Powershell for Windows):

scp -r username@hostname:/data/module5/toy your/local/path

2.2 Large data set

This data set is based on real genotypes at 273,604 SNPs from the UK Biobank, with a simulated trait affected by 100 causal variants and heritability 0.5. The samples are split into three independent cohorts:

  • Discovery population for conducting GWAS (n = 3,000).
  • Testing population for finding the best prediction model (n = 500).
  • Target population for prediction (n = 100).

Data location:

ls /data/module5/ukb/

The GWAS summary statistics are stored in gwas.ma:

head /data/module5/ukb/gwas.ma
## SNP A1 A2 freq b se p N
## rs10000010 C T 0.477993 -0.2796 0.3719 0.4522 2999
## rs10000030 A G 0.135302 -0.09771 0.5496 0.8589 2997
## rs1000007 C T 0.275333 -0.1548 0.4153 0.7094 3000
## rs10000121 G A 0.468092 -0.5717 0.3729 0.1253 2993
## rs10000141 A G 0.087725 0.6495 0.6649 0.3287 2998
## rs1000016 G A 0.068402 0.6553 0.7509 0.3829 2997
## rs10000169 C T 0.250833 -0.1971 0.4317 0.648 3000
## rs1000022 T C 0.415138 0.07813 0.3815 0.8377 2999
## rs10000272 C T 0.056985 0.8443 0.8181 0.3024 2992

3 Principle illustration using a small data set

Let’s start with using R to calculate PGS in the small data set. Start R and load the data as shown above.

Perform GWAS in the training population by regressing the phenotype on each SNP one at a time. This gives marginal SNP effect estimate (bhat), ignoring linkage disequilibrium (LD) between SNPs.

fit = apply(X, 2, function(x){summary(lm(y ~ x))$coef[2,]})
bhat = fit[1,]

Now, using the validation genotype matrix Xval and the SNP effect estimates bhat,

Question: Can you write a R code to predict the PGS for the 31 validation individuals? [Hint: PGS is a weighted sum of allele counts across SNPs.]

# Given the SNP effect vector
print(bhat)
##         V1         V2         V3         V4         V5         V6         V7 
##  2.6026705  0.6301389  0.6301389  2.2707780  2.0718035  2.2161266  0.5976793 
##         V8         V9        V10 
## -2.2707780  0.8994498 -2.2707780
# and validation genotype matrix
head(Xval)
##      V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
## [1,]  1  1  1  1  1  1  2  1  0   1
## [2,]  1  0  0  1  1  1  2  1  0   1
## [3,]  1  0  0  1  1  1  2  1  0   1
## [4,]  1  0  0  1  1  1  2  1  0   1
## [5,]  0  0  0  0  0  0  1  2  0   2
## [6,]  1  1  1  1  1  1  2  1  0   1


Take a moment to think about it before continuing.




Suppose we predict PGS using all SNPs with their marginal effect estimates from GWAS:

PGS = Xval %*% bhat

Let’s check the prediction accuracy by taking the square of correlation between phenotype and PGS.

print(cor(yval, PGS)^2)
##         [,1]
## V1 0.3696879

Question: Is this prediction accuracy high? [Hint: What value do you expect if PGS can perfectly predict the true genetic value?]


Take a moment to think about it before continuing.




Answer: Not too high, as the upper bound of prediction accuracy is the trait heritability (0.5), which determines the maximum proportion of variation that can be explained by the genetic predictor.

Let’s also check the bias of PGS by regressing the phenotype on the PGS.

# center PGS and phenotype to remove the intercept because the intercept is irrelevant
PGS = PGS - mean(PGS)
yval = yval - mean(yval)

# plot phenotype against PGS
plot(PGS, yval, xlab="PGS using all SNPs", ylab="Phenotype")
abline(a=0, b=1)  # add y=x line
abline(lm(yval ~ PGS), col="red")  # add regression line

# you can check the value of the regression slope by
print(lm(yval ~ PGS))
## 
## Call:
## lm(formula = yval ~ PGS)
## 
## Coefficients:
## (Intercept)          PGS  
##   2.152e-17    1.675e-01

You may have noticed that the regression slope of phenotypes on PGS is substantially smaller than one.

Question: What does it mean? [Hint: slope < 1 means that one unit increase in x leads to less than an unit increase in y. Does PGS tend to be inflated or deflated relative to the true phenotype value?]


Take a moment to think about it before continuing.




Answer: Ideally, we want the slope to be equal to one, which would indicate that an unit change in PGS leads to an unit change in phenotype. This is often referred to as “unbiasedness”. When the slope is below one, it means the PGS is upward biased, i.e., the absolute predicted value is larger than the true value. You can also see this from the difference in the range of values between x and y axes.


Question: What could be the cause(s) of such a relative low prediction accuracy and a large bias?


Take a moment to think about it before continuing.




3.1 LD pruning

To answer the question above, let’s check the LD correlations between SNPs.

R = cor(X)
print(R)
##             V1         V2         V3         V4          V5          V6
## V1   1.0000000  0.1338359  0.1338359  0.7284017  0.57898634  0.64984918
## V2   0.1338359  1.0000000  1.0000000  0.3789184  0.25731616  0.36895246
## V3   0.1338359  1.0000000  1.0000000  0.3789184  0.25731616  0.36895246
## V4   0.7284017  0.3789184  0.3789184  1.0000000  0.81429611  0.81303937
## V5   0.5789863  0.2573162  0.2573162  0.8142961  1.00000000  0.90141987
## V6   0.6498492  0.3689525  0.3689525  0.8130394  0.90141987  1.00000000
## V7   0.1023013  0.1726592  0.1726592  0.2053259  0.10824685  0.18231275
## V8  -0.7284017 -0.3789184 -0.3789184 -1.0000000 -0.81429611 -0.81303937
## V9   0.3062689  0.4491011  0.4491011  0.4820692  0.07706692  0.05818721
## V10 -0.7284017 -0.3789184 -0.3789184 -1.0000000 -0.81429611 -0.81303937
##              V7         V8          V9        V10
## V1   0.10230128 -0.7284017  0.30626887 -0.7284017
## V2   0.17265922 -0.3789184  0.44910112 -0.3789184
## V3   0.17265922 -0.3789184  0.44910112 -0.3789184
## V4   0.20532585 -1.0000000  0.48206922 -1.0000000
## V5   0.10824685 -0.8142961  0.07706692 -0.8142961
## V6   0.18231275 -0.8130394  0.05818721 -0.8130394
## V7   1.00000000 -0.2053259  0.07758962 -0.2053259
## V8  -0.20532585  1.0000000 -0.48206922  1.0000000
## V9   0.07758962 -0.4820692  1.00000000 -0.4820692
## V10 -0.20532585  1.0000000 -0.48206922  1.0000000

You can see that some SNPs are in high LD correlations with others. This may result in double counting of SNP effects because the SNP effects estimated from GWAS do not account for LD between SNPs.

Below is a simple R code to remove one of the SNPs in pairwise LD \(r^2\) greater than a threshold (e.g., 0.5). This procedure is known as LD pruning.

rsqThreshold = 0.5
nsnp = ncol(X)
removed = c()
for (i in 1:(nsnp-1)) {
  if (! i %in% removed) { 
    for (j in (i+1):nsnp) {
      if(R[i,j]^2 > rsqThreshold) removed = c(removed, j)
    }
  }
}
removed = unique(removed)
all = 1:nsnp
keep = all[!all %in% removed] # remaining SNPs after LD pruning
print(R[keep, keep])  # LD correlations among LD pruned SNPs
##           V1        V2         V5         V7         V9
## V1 1.0000000 0.1338359 0.57898634 0.10230128 0.30626887
## V2 0.1338359 1.0000000 0.25731616 0.17265922 0.44910112
## V5 0.5789863 0.2573162 1.00000000 0.10824685 0.07706692
## V7 0.1023013 0.1726592 0.10824685 1.00000000 0.07758962
## V9 0.3062689 0.4491011 0.07706692 0.07758962 1.00000000

Now, let’s use the LD-pruned SNPs only to calculate PRS and check the prediction accuracy and bias.

PGS2 = Xval[,keep] %*% bhat[keep]
print(cor(yval, PGS2)^2)
##         [,1]
## V1 0.4361839
PGS2 = PGS2 - mean(PGS2)

plot(PGS2, yval, xlab="PGS using LD-pruned SNPs", ylab="Phenotype")
abline(a=0, b=1)
abline(lm(yval~PGS2), col="red")

print(lm(yval ~ PGS2))
## 
## Call:
## lm(formula = yval ~ PGS2)
## 
## Coefficients:
## (Intercept)         PGS2  
##   2.390e-16    5.428e-01


Question: Are the prediction accuracy and bias improved? If so, why?


Take a moment to think about it before continuing.




Answer: Yes, both prediction accuracy and bias are improved. This improvement occurs because LD pruning helps to avoid double counting of SNP effects.


3.2 P-value thresholding

Checking the P-values of the LD-pruned SNPs finds that some SNP effects are statistically insignificant different from zero.

pval = fit[4,]
print(pval[keep])
##           V1           V2           V5           V7           V9 
## 4.491829e-49 6.581417e-02 8.597549e-23 5.433139e-02 5.818065e-04

Including SNPs that do not affect the phenotype in the PGS may add noise to the prediction. But where is the cutoff line? Let’s try to keep genome-wide significant (P-value < 5e-8) SNPs only in the equation. This step of filtering is called P-value thresholding.

pvalThreshold = 5e-8
keepFinal = keep[keep %in% which(pval < pvalThreshold)]
PGS3 = Xval[,keepFinal] %*% bhat[keepFinal]

print(cor(yval, PGS3)^2)
##         [,1]
## V1 0.4794083
PGS3 = PGS3 - mean(PGS3)
plot(PGS3, yval, xlab="PGS using LD-pruned and P-value thresholded SNPs", ylab="Phenotype")
abline(a=0, b=1)
abline(lm(yval~PGS3), col="red")

print(lm(yval ~ PGS3))
## 
## Call:
## lm(formula = yval ~ PGS3)
## 
## Coefficients:
## (Intercept)         PGS3  
##   3.378e-16    6.651e-01


Question: Are the prediction accuracy and bias improved? If so, why?


Take a moment to think about it before continuing.




You may try other significance thresholds and see how the result changes.


In conclusion, LD pruning and P-value thresholding improves prediction accuracy and unbiasness by avoiding double counting of SNP effects and excluding possible false positives. However, in practice, it is often not clear where to draw the line between the false and true positives. In addition, LD clumping is usually better than LD pruning, because the clumping algorithm takes P-values into account when filtering SNPs in LD. In the next session, we will use PLINK to carry out the clumping + P-value thresholding (C+PT) method for the analysis of the simulated UKB data.