Practical 2 Best Linear Unbiased Prediction

1 Objectives

In this practical, we will first investigate the overfitting issue in C+PT when target sample is used for SNP selection. We will then learn how to perform Best Linear Unbiased Prediction (BLUP). BLUP is a method to estimate the effects of all SNPs simultaneously by treating them as random effects. Unlike ordinary least squares (OLS), a standard method used in GWAS, BLUP does not require SNP selection or filtering.

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

2 Overfitting the threshold and a tuning split for C+PT

Q1: 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?

Hint

PRSice’s default workflow tries many P-value thresholds and reports R² for each using the same target sample. You then pick whichever threshold looks best. That “winner” reflects real genetic architecture, but random environmental noise in this validation sample can also favour one threshold over another. Is that a problem when you interpret the reported R²?

Show answer

No — you would not expect the same R². PRSice’sdefault workflow tries many P-value thresholds on the same target sample and picks the best. The same random noise in that sample contributes to both (1) choosing the threshold(which SNPs enter the PGS) and (2) evaluating prediction accuracy (the reported R²), sothe winning row’s R² is usually optimistically high compared with fixing the threshold in advance and testing once on new samples (winner’s curse / overfitting the tuning step). In anew cohort with the same genetic architecture, environmental residuals are new, so the noise that favoured a particular threshold in the first cohort does not recur and R² will differ.

Tuning vs held-out target: One solution is to use a small tuning subset of individuals only to pick a P-value threshold (e.g. best R² among bars). Then apply only that threshold on a held-out target subset for a cleaner out-of-sample comparison.

The examples use PLINK --keep files named ind_94.txt and ind_400.txt in your current working directory (two columns: FID and IID, no header).

3 Tuning in an indepdent sample

TERMINAL

PRSice_linux \
--a1 A1 \
--base gwas.ma \
--target 1000G_phase3.eur.QC.unrel \
--pheno target_phenotypes.txt \
--cov target_covariates.txt \
--keep ind_94.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 tune

Open tune.prsice and note which threshold gives the best R² on this subset (record the threshold – you will need it in the next step B2).

4 Held-out target validation at the chosen threshold only

Run PRSice again on the held-out individuals (ind_400.txt), but only at the single threshold chosen in B1. The command below uses 3e-5 as a concrete example (it matches the prepared teaching solution); if your best threshold from B1 differs, change --bar-levels to that value and keep everything else the same.

TERMINAL

PRSice_linux \
--a1 A1 \
--base gwas.ma \
--target 1000G_phase3.eur.QC.unrel \
--pheno target_phenotypes.txt \
--cov target_covariates.txt \
--keep ind_400.txt \
--beta \
--bar-levels 3e-5 \
--pvalue p \
--stat b \
--binary-target F \
--fastscore \
--out heldout_target

Q2: Read R2 from the .prsice files (not PRS.R2 from .summary — see Part A3). Open output.prsice (Part A, whole target), tune.prsice, and heldout_target.prsice. Each sample can have a different P-value threshold (e.g. best on the whole target in Part A, best on the tuning subset in B1, and your B1 choice applied in B2 on the held-out subset). Record the threshold and R^2^ you are comparing in the table below.

	Whole target (output.prsice)	Tuning (tune.prsice, N=94)	Held-out (heldout_target.prsice, N=400)
P-value
R²

Show answer

	Whole target (output.prsice)	Tuning (tune.prsice, N=94)	Held-out (heldout_target.prsice, N=400)
P-value	3e-5	3e-5	3e-5
R²	0.168322	0.246327	0.152673

Q3: If the validation dataset was much larger, would you expect the tuning vs held-out gap to shrink, grow, or stay the same, and why? What happens if you let many thresholds compete again on a very small tuning set?

Show answer

The gap would usually shrink. Larger samples give more stable R² estimates and less room for noise to make the tuning subset look unusually good. Some difference can remain (sampling variation, subtle population differences). When the sample size of the tuning set is very small, R² becomes unstable and overfitting worsens: noise can pick a threshold that looks excellent on ~94people but fails on held-out data. That is why a separate tuning set and a single threshold on held-out individuals are used.

5 Toy data for BLUP

We will use a small toy data set to illustrate the basic principles and manually compute BLUP estimates in R. This data set is available to download. You can also use R on our server to do this exercise.

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.

6 Estimate SNP effects using BLUP

Let’s start by loading the data in R.

# data for training population
X <- as.matrix(read.table("xmat_trn.txt"))
y <- as.matrix(read.table("yvec_trn.txt"))

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

BLUP estimates SNP effects using the following system of equations:

\[\begin{equation} \begin{bmatrix} 1_n'1_n & 1_n'X \\ \\ X'1_n & X'X+I\lambda \end{bmatrix} \begin{bmatrix} \mu \\ \beta \end{bmatrix} = \begin{bmatrix} 1_n'y \\ X'y \end{bmatrix} \end{equation}\]

where \(1_n\) is a vector of ones \((325 \times 1)\), \(X\) is the genotype matrix \((325 \times 10)\), \(I\) is an identity matrix, \(\lambda = \frac{1 - h^2}{h^2/m}\) is the shrinkage parameter, \(\mu\) is the overall mean, and \(\beta\) is the vector of SNP effects. In the literature, this system of equations is known as mixed model equations (MME), and the big matrix on the left is called coefficient matrix.

First, we need to know the value of \(\lambda\).

Q4: What’s the interpretation of \(\lambda\)? Can you write a code to compute it?

Show answer

\(\lambda\) is the ratio of residual variance to the random effect variance, which controls the degree to which random effect estimates shrunk toward zero. The higher the \(\lambda\), the more uncertain we are about the SNP effects, so we shrink them more strongly toward zero. It can be estimated as \((1-h^2)/(h^2/m) = 10\) where \(m=10\) is the number of SNPs.

h2 = 0.5           # trait heritability
nind = nrow(X)     # number of training individuals
nsnp = ncol(X)     # number of SNPs
lambda = (1-h2)/(h2/nsnp)
lambda

## [1] 10

We already have matrix X and vector y. We can use the following code to solve the system of BLUP equations:

blup = function(X, y, lambda){
  nsnp = ncol(X)
  D = diag(c(0, rep(lambda,nsnp))) # first element is 0 because we do not add shrinkage parameter to the intercept
  W = cbind(1, X)
  coeff = crossprod(W) + D
  rhs = crossprod(W, y)
  return(solve(coeff, rhs))
}

solution_blup = blup(X, y, lambda)
solution_blup

##              V1
##      0.64651875
## V1   1.90752815
## V2  -0.14131984
## V3  -0.14131984
## V4   0.16483065
## V5   0.22968939
## V6   0.27095334
## V7   0.09624804
## V8  -0.16483065
## V9  -0.07477323
## V10 -0.16483065

Q5: What is the solution for the mean? Which SNP has the largest effect? Which SNP has the second largest effect? Does the result make sense to you?

Hint

The causal variants are SNP 1 and SNP 5 with causal effect size of 2 and 1, respectively. What’s the LD correlation between SNP 5 and 6?

cor(X)

##             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

Show answer

The 1st SNPs is a causal variant with a large effect (true effect size = 2). The BLUP estimate of its effect is 1.9. Nearly spot on! The 5th SNP is the another causal variant with a smaller effect (true effect size = 1). However, the second largest BLUP estimate is at the 6th SNP (BLUP estimate = 0.27) and the third largest is at the 5th SNP (BLUP estimate = 0.23). This is because SNP 5 and 6 are in high LD (LD correlation 0.9) and SNP 6 has a larger effect estimate due to sampling. Besides SNP 1 and 5, the other SNPs are null SNPs, but their BLUP estimates are not zero. This is because BLUP assumes that all SNPs contribute to the trait with an equal effect variance. As a result, the effect of causal variant is smeared over SNPs in LD with it.

7 Predict polygenic scores with BLUP

Q6: The genotypes for the validation individuals are stored in Xval. Can you write a R code to calculate their polygenic scores?

Show answer

We use all SNPs with their BLUP solutions to predict PGS. No SNP selection or filtering is needed.

PGS_BLUP = Xval %*% solution_blup[-1]  # the first element is the intercept estimate

The prediction accuracy is calculated as the squared correlation between phenotypes and PGS in the validation sample. The prediction bias is represented by the regression slope of validation phenotype on PGS.

# prediction accuracy
print(cor(PGS_BLUP, yval)^2)

##             V1
## [1,] 0.5450557

# prediction bias
yval = yval - mean(yval)
PGS_BLUP = PGS_BLUP - mean(PGS_BLUP)

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

print(lm(yval ~ PGS_BLUP))

## 
## Call:
## lm(formula = yval ~ PGS_BLUP)
## 
## Coefficients:
## (Intercept)     PGS_BLUP  
##   5.709e-17    1.124e+00

Practical 2 Best Linear Unbiased Prediction

Genetics & Genomics Winter School | Module 5 Polygenic Prediction

2026-07-03

1 Objectives

2 Overfitting the threshold and a tuning split for C+PT

3 Tuning in an indepdent sample

4 Held-out target validation at the chosen threshold only

5 Toy data for BLUP

6 Estimate SNP effects using BLUP

7 Predict polygenic scores with BLUP