1 Objectives

In this practical, we will learn how to perform polygenic prediction using a multi-component mixture Bayesian model using individual-level data (BayesR) or summary statistics (SBayesR). As discussed in the lecture, BLUP can be regarded as a Bayesian model assuming a normal prior distribution for the SNP effects. With different prior assumptions, especially with a mixture prior, Bayesian methods are flexible in accounting for the genetic architecture of the trait. We will explore the principle of Bayesian methods using a small data set and R script (as we did in BLUP session), and then will employ GCTB (https://cnsgenomics.com/software/gctb/#Overview) to perform the simulated data analysis.

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

2 Optional reading: Understanding the algorithm

Show optional reading

This section helps you understand the R code used in the method.

Different Bayesian methods vary only in the prior specification for SNP effects. Here, we focus on two Bayesian methods: BayesC and BayesR. BayesC is a special case of BayesR with a two-component mixture. The algorithm for BayesR (i.e., Markov chain Monte Carlo sampling scheme) is implemented in the R script bayesr.R.

cat /data/module5/bayesr.R

Do not run the R code below in this section.

BayesR assumes each SNP effect follows a mixture of normal distributions. The input parameters include

# Input parameters for bayesr. Do not run this.
bayesr = function(X, y, niter, gamma, startPi, startH2){
  ...
}
  • niter is the number of iterations for MCMC sampling
  • X is the genotype matrix
  • y is the vector of phenotypes
  • gamma is a vector of scaling factors for the variance of the mixture components
  • startPi is a vector of starting values for the mixture proportions \(\pi\)
  • startH2 is the starting value of SNP-based heritability.

The number of elements in gamma and startPi determines the number of mixture components, and they must match.

Let’s dive into more details to understand what the code does.

The first step is to declare and initialise variables. This step also adjusts y for the population mean and SNP effects (which are initially set to zero), so that ycorr represents the residuals.

  # Step 1. Declaration and initialisation of variables
  n = nrow(X)          # number of observations
  m = ncol(X)          # number of SNPs
  ndist = length(startPi)  # number of mixture distributions
  pi = startPi         # starting value for pi
  h2 = startH2         # starting value for heritability
  vary = var(y)        # phenotypic variance
  varg = vary*h2       # starting value for genetic variance
  vare = vary*(1-h2)   # starting value for residual variance
  sigmaSq = varg/(m*sum(gamma*pi))    # common factor of SNP effect variance
  nub = 4              # prior degrees of freedom for SNP effect variance
  nue = 4              # prior degrees of freedom for residual variance
  scaleb = (nub-2)/nub*sigmaSq  # prior scale parameter for SNP effect variance
  scalee = (nue-2)/nue*vare     # prior scale parameter for residual variance
  beta = array(0,m)    # vector of SNP effects
  beta_mcmc = matrix(0,niter,m) # MCMC samples of SNP effects
  mu = mean(y)         # overall mean
  xpx = apply(X, 2, crossprod)  ## diagonal elements of X'X
  probDelta = vector("numeric", ndist)
  logDelta = array(0,2)
  keptIter = NULL      # keep MCMC samples
  
  # adjust/correct y for population mean and all SNP effects (which all have zero as initial value)
  ycorr = y - mu

Then the Gibbs sampling begins. Gibbs sampling iteratively draws samples from the full conditional distribution of each parameter (i.e., posterior distribution conditional on the values of other parameters). Over time, this converges to the joint posterior distribution.

The first step within each iteration is to sample the population mean (fixed effect) from its full conditional distribution, assuming a flat prior. This distribution is centered at the BLUE (best linear unbiased estimator) solution, with variance equal to the inverse of the coefficient matrix in the “mini” mixed model equation for the mean.

  # MCMC begins
  for (iter in 1:niter){
    # Step 2. Sample the mean from a normal posterior
    ycorr = ycorr + mu   # unadjust y with the old sample
    rhs = sum(ycorr)     # right hand side of the mixed model equation
    invLhs = 1/n         # inverse of left hand side of the mixed model equation
    muHat = invLhs*rhs   # BLUE solution
    mu = rnorm(1, muHat, sqrt(invLhs*vare))   # posterior is a normal distribution
    ycorr = ycorr - mu   # adjust y with the new sample

The next step is to sample the mixture component indicator variable and effect size for each SNP, conditional on the effects of all the other SNPs. We first sample the indicator variable delta for the mixture distribution membership (multinomial posterior), and then sample the SNP effect beta from the corresponding normal distribution or a point mass at zero.

    # Step 3. Sample each SNP effect from a multi-normal mixture distribution
    logPi = log(pi)
    logPiComp = log(1-pi)
    invSigmaSq = 1/(gamma*c(sigmaSq))
    logSigmaSq = log(gamma*c(sigmaSq))
    nsnpDist = rep(0, ndist)  # number of SNPs asigned to each mixture component
    ssq = 0  # weighted sum of squares of SNP effects 
    nnz = 0  # number of nonzero effects
    ghat = array(0,n)  # individual genetic value
    # loop over SNPs
    for (j in 1:m){
      oldSample = beta[j]
      rhs = crossprod(X[,j], ycorr) + xpx[j]*oldSample  # right hand side of the mixed model equation
      rhs = rhs/vare
      invLhs = 1.0/(xpx[j]/c(vare) + invSigmaSq)        # inverse of left hand side of the mixed model equation
      betaHat = invLhs*c(rhs)                           # BLUP solution
      
      # sample the mixture distribution membership
      logDelta = 0.5*(log(invLhs) - logSigmaSq + betaHat*c(rhs)) + logPi  # log likelihood + prior for nonzero effects
      logDelta[1] = logPi[1]                                              # log likelihood + prior for zero effect
      for (k in 1:ndist) {
        probDelta[k] = 1.0/sum(exp(logDelta - logDelta[k]))   # posterior probability for each distribution membership
      }
      
      delta = sample(1:ndist, 1, prob = probDelta)   # indicator variable for the distribution membership
      nsnpDist[delta] = nsnpDist[delta] + 1
      
      if (delta > 1) {
        beta[j] = rnorm(1, betaHat[delta], sqrt(invLhs[delta]))  # given the distribution membership, the posterior is a normal distribution
        ycorr = ycorr + X[,j]*(oldSample - beta[j])              # update ycorr with the new sample
        ghat = ghat + X[,j]*beta[j]                              
        ssq = ssq + beta[j]^2 / gamma[delta]
        nnz = nnz + 1
      } else {
        if (oldSample) ycorr = ycorr + X[,j]*oldSample           # update ycorr with the new sample which is zero
        beta[j] = 0
      }
    }   
    beta_mcmc[iter,] = beta

Next, we sample the other parameters in the model, including mixing probabilities (\(\pi\)), SNP effect variance (\(\sigma^2_{\beta}\)), and residual variance (\(\sigma^2_e\)). The full conditional distribution for \(\pi\) is a beta distribution. For \(\sigma^2_{\beta}\) and \(\sigma^2_e\), they are scaled inverse chi-square distributions.

    # Step 4. Sample the distribution membership probabilities from a Dirichlet distribution
    pi = rdirichlet(1, nsnpDist + 1)
    
    # Step 5. Sample the SNP effect variance from a scaled inverse chi-square distribution
    sigmaSq = (ssq + nub*scaleb)/rchisq(1, nnz+nub)
    
    # Step 6. Sample the residual variance from a scaled inverse chi-square distribution
    vare = (crossprod(ycorr) + nue*scalee)/rchisq(1, n+nue)

Given the genetic values and residual variance, we can compute the SNP-based heritability. This is an appealing property of the MCMC approach – given the sampled values of SNP effects, you can easily compute any statistics of interest.

    # Step 7. Compute the SNP-based heritability
    varg = var(ghat)
    h2  = varg/(varg + vare)

The final step in each iteration is to store the MCMC samples for the parameter values

    keptIter <- rbind(keptIter,c(pi, nnz, sigmaSq, h2, vare, varg))

After completing all MCMC iterations, we compute posterior means as point estimates.

  colnames(keptIter) <- c(paste0("Pi", 1:length(pi)),"Nnz","SigmaSq","h2","Vare","Varg")
  postMean = apply(keptIter, 2, mean)  # posterior mean of MCMC samples is used as the point estimate for each parameter
  cat("\nPosterior mean:\n")
  print(postMean)
  return(list(par=keptIter, beta=beta_mcmc))

3 Illustration of the principle using BayesR

Load 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"))

4 BayesR with a point-normal prior (two components)

We will first analyse this small data set using a two-component BayesR with a point-normal prior, which is also known as BayesC in the literature.

\[\beta_j \begin{cases} \sim N(0, \sigma^2_\beta) & \text{ with probability } \pi \\ \\ = 0 & \text{ with probability } 1 - \pi \end{cases} \] where parameter \(\pi\) is estimated from the data. When \(\pi = 0\), the prior collapses to a normal distribution, i.e., \(\beta_j \sim N(0, \sigma^2_\beta)\), which matches the assumption in BLUP.

BayesR has been implemented in BayesR.R. Let’s load the script in R.

source("bayesr.R")

Specify the model parameters:

gamma = c(0, 1)
startPi = c(0.8, 0.2)

This assumes that 20% of SNPs have nonzero effects drawn from a normal distribution, while 80% have zero effect. Now run the MCMC:

res_bayesc = bayesr(X = X, y = y, gamma = gamma, startPi = startPi)
## 
##  iter  100, nnz =    4, sigmaSq =  2.120, h2 =  0.477, vare =  3.792, varg =  3.454 
## 
##  iter  200, nnz =    3, sigmaSq =  1.253, h2 =  0.506, vare =  3.719, varg =  3.816 
## 
##  iter  300, nnz =    2, sigmaSq =  4.430, h2 =  0.543, vare =  3.921, varg =  4.651 
## 
##  iter  400, nnz =    4, sigmaSq =  0.869, h2 =  0.466, vare =  3.908, varg =  3.412 
## 
##  iter  500, nnz =    3, sigmaSq =  1.411, h2 =  0.541, vare =  3.899, varg =  4.590 
## 
##  iter  600, nnz =    5, sigmaSq =  0.810, h2 =  0.450, vare =  4.463, varg =  3.647 
## 
##  iter  700, nnz =    2, sigmaSq =  1.779, h2 =  0.495, vare =  3.927, varg =  3.848 
## 
##  iter  800, nnz =    2, sigmaSq =  1.120, h2 =  0.496, vare =  3.875, varg =  3.810 
## 
##  iter  900, nnz =    3, sigmaSq =  1.149, h2 =  0.450, vare =  3.683, varg =  3.010 
## 
##  iter 1000, nnz =    4, sigmaSq =  0.810, h2 =  0.494, vare =  3.963, varg =  3.862 
## 
## Posterior mean:
##       Pi1       Pi2       Nnz   SigmaSq        h2      Vare      Varg 
## 0.5606956 0.4393044 4.2410000 1.5557804 0.4938728 3.9202915 3.8386095

The output includes sampled values of key parameters every 100 iterations. For example, Nnz is the number of nonzero effects fitted in the model. Unlike BLUP, which assumes all SNPs contribute equally to the trait, BayesC allows some SNPs to have zero effect.

After MCMC, you can find the sampled values for the model parameters and SNP effects for each iteration in the result list. For example, you can summarise the posterior mean and standard deviation for each parameter by

colMeans(res_bayesc$par)
##       Pi1       Pi2       Nnz   SigmaSq        h2      Vare      Varg 
## 0.5606956 0.4393044 4.2410000 1.5557804 0.4938728 3.9202915 3.8386095
apply(res_bayesc$par, 2, sd)
##        Pi1        Pi2        Nnz    SigmaSq         h2       Vare       Varg 
## 0.20817722 0.20817722 1.93769803 1.17203485 0.03554869 0.30531519 0.44132465

The posterior mean gives a point estimate, and the standard eviation gives an estimate of uncertainty (posterior standard error).

You can plot parameter traces over iterations to check convergence. For example:

# Trace plot for SNP 1, which is the causal variant with effect size of 2
plot(1:nrow(res_bayesc$beta), res_bayesc$beta[,1], xlab="Iteration", ylab="beta[,1]") 
abline(h = 2, col="red")

Q1: Use this code to investigate each parameter. Do they fluctuate around the true value?

# Trace plot for SNP 2, which is a null SNP
plot(1:nrow(res_bayesc$beta), res_bayesc$beta[,2], xlab="Iteration", ylab="beta[,2]") 
abline(h = 0, col="red")

# Trace plot for SNP 5, which is a causal variant with a smaller effect of 1 and in LD with other SNPs
plot(1:nrow(res_bayesc$beta), res_bayesc$beta[,5], xlab="Iteration", ylab="beta[,5]") 
abline(h = 1, col="red")

Show answer
  • For SNP 1, the trace plot fluctuates around 2, consistent with its true effect.
  • For SNP 5, the plot fluctuates around 1, with more variability due to LD with nearby SNPs.
  • For SNP 2 (null SNP), the plot fluctuates around 0.

This indicates good mixing and convergence.


We can also plot the posterior distribution for each SNP effect. We discard the first 100 iterations of the program as “burn in”:

# Posterior distribution of SNP 1 effect
plot(density(res_bayesc$beta[100:1000,1]))

Q2: Does the distribution appear to be normal? What about the distributions of the other parameters?

# Posterior distribution of SNP 2 effect
plot(density(res_bayesc$beta[100:1000,2]))

Show answer
  • For causal SNPs, the posterior distribution of effects is approximately normal but may be slightly skewed or truncated due to the mixture prior.
  • For null SNPs, the posterior is often a spike at zero, reflecting the point-mass component of the prior.
  • For other parameters like heritability and variance components, the posteriors are not necessarily symmetric, especially with small sample sizes.


The posterior distribution of the SNP effect can be summarised to quantify the evidence of the SNP associated with the trait. From the posterior samples, we can calculate the Bayesian 90% credible interval, which represents 95% probability that the true parameter would lie within the interval, given the evidence from the data.

quantile(res_bayesc$beta[100:1000,1], probs=c(0.05, 00.95))
##       5%      95% 
## 1.711408 2.396162

Q3: What can we conclude from the Bayesian credible interval?

Show answer

The 90% credible interval shows the range of plausible effect sizes for a SNP, given the data and prior. If the interval excludes 0, there’s strong evidence the SNP is associated with the trait. For null SNPs, the interval will typically include 0 or be sharply peaked at zero.


Another way to measure SNP-trait association is to compute the posterior inclusion probability (PIP), the probability of the SNP being fitted as a non-zero effect in the model across MCMC samples. It has been widely used in the fine-mapping analysis.

For example, PIP for SNP 1 can be calculated as

mean(res_bayesc$beta[100:1000,1] != 0)
## [1] 1

This could be generalised for all SNPs:

pip = colMeans(res_bayesc$beta[100:1000, ] != 0)
plot(pip, type="h", xlab="SNP", ylab="Posterior inclusion probability")

Q4: Which SNPs are likely to be the causal variants based on PIP?

Show answer

SNPs with high PIP (close to 1) are likely to be causal (e.g., SNP 1 and SNP 5 in your simulation). A plot of PIP across SNPs will highlight those with strong support. SNPs with PIP near 0 are unlikely to be associated.


To get the SNP effect estimates, we calculate the posterior means of SNP effects:

betaMean = colMeans(res_bayesc$beta)


Q5: How to predict PGS using BayesC effect estimates? How does joint modelling of genome-wide SNPs in BayesR differ from ranking SNPs by P-value and clumping in C+PT?

Show answer 
ghat_bayesc=Xval %*% betaMean
summary(lm(yval ~ ghat_bayesc))$r.squared
## [1] 0.5357362

BayesR fits a genome-wide Bayesian model on GWAS summary statistics and re-estimates joint SNP effects accounting for LD, whereas C+PT selects a subset of SNPs by marginal P-value and LD clumping and builds the PGS from those SNPs’ marginal GWAS effects. The \(R^2\) from BayesR reflects posterior mean effects and typically outperforms BLUP or C+PT in sparse genetic architectures (few large-effect SNPs). The difference is due to BayesC allowing effect sizes to be exactly zero, avoiding overfitting noise.


5 Large simulated data analysis using GCTB

We now apply SBayesR to the simulated data based on real genotypes at 273,604 SNPs from the UK Biobank using GCTB.

6 LD reference

SBayesR need LD information among SNPs, not just marginal GWAS results. The original SBayesR paper (2019) used chromosome-wide shrunk LD matrices. Here we use eigen-decomposed block LD matrices. They tend to be more stable and often converge more reliably GCTB SBayesRC.

GCTB expects a folder such as ldm_eigen/ with binary eigenvalue/eigenvector files per LD block plus text files snp.info and ldm.info. A ready-made ldm_eigen/ folder is provided for this practical.

You cannot read the binary eigen files directly, but you can inspect the text summaries:

TERMINAL

head ldm_eigen/snp.info
head ldm_eigen/ldm.info

For your own projects: the teaching reference uses fewer than 300k SNPs. GCTB also supports much larger LD builds for real studies eigen LD. To build LD from genotypes, see the tutorial; blocks often follow ref4cM_v37.pos (Berisa & Pickrell). Example commands (no need to run this for the practical):

bfile="proper_genotype_file"
blockinfo="ref4cM_v37.pos"
gctb --bfile $bfile \
--make-block-ldm \
--block-info $blockinfo \
--thread 4 \
--out ldm

gctb --ldm ldm --make-ldm-eigen --thread 4 --out ldm

You can change --thread to match your machine. On an HPC cluster, people can run one LD block per array job (details).

7 Impute summary statistics, then fit SBayesR

Eigen LD requires a GWAS row for every SNP in the LD reference. If any SNPs in the LD reference are missing in your summary statistics file, GCTB can impute them using LD and the SNPs that do match. The –out flag sets the output prefix.

TERMINAL

gctb \
--ldm-eigen ldm_eigen \
--gwas-summary gwas.ma \
--impute-summary \
--thread 4 \
--out gwas

You should obtain gwas.imputed.ma. During imputation, GCTB also does basic QC (allele matching and dropping SNPs whose per-SNP N differs from the median by more than three standard deviations).

On an HPC cluster: you can impute one block per job with --block $i.

When imputation is done, run SBayesR:

gctb \
--ldm-eigen ldm_eigen \
--gwas-summary gwas.imputed.ma \
--sbayes R \
--thread 4 \
--out sbayesr

Increase --thread if your machine has spare cores. --sbayes R selects the SBayesR mixture model (C would select SBayesC, a spike-and-slab model, equivalent to a two-component SBayesR).

SBayesR defaults to five mixture components (a point mass at zero plus four normal components with increasing variance). That usually works well. You can specify an alternative model. For example, to run four components, you can add flags (no need to run this for the practical):

--pi 0.95,0.02,0.02,0.01 \
--gamma 0.0,0.01,0.1,1

Here --pi gives prior probabilities for each component (comma-separated, must sum to 1) and --gamma sets the scale of each component as a fraction of SNP-based heritability. The number of values must match in both flags. While GCTB runs, the screen prints MCMC summaries every 100 iterations, including numbers of SNPs per component (NumSnp1-5), variance explained per component (Vg1-5), SNP effect variance (SigmaSq), SNP-based heritability (hsq), residual variance (ResVar), and NumSkeptSnp (SNPs whose shrunk joint effect looks inconsistent with the marginal effect and therefore would be removed), and TimeLeft (predicted time to finish).

Q6: Do the MCMC traces for key parameters look stable? What might make them noisier?

Show answer

In an ideal run, SNP-based heritability, pi, and variance components should settle after burn-in without strong drift. Traces look noisier with shorter chains, weaker GWAS signal, or a prior that is too flexible for the data (many mixture components fighting for the same signal).


By default GCTB runs 3,000 iterations and discards the first 1,000 as burn-in. For this practical it takes about three minutes. You can change length and burn-in (chain length must exceed burn-in), for example:

--chain-length 1000 \
--burn-in 500

Main outputs are sbayesr.parRes (parameter posteriors) and sbayesr.snpRes (per-SNP results).

sbayesr.parRes gives posterior means and standard deviations for model parameters (posterior mean gives point estimate and SD measures uncertainty).

Q7: What is the estimated SNP-based heritability (hsq)? Is it different from the simulation value (0.5)? What could pull it away from 0.5?

Show answer

The SNP-based heritability estimate is slightly over 0.5. That can happen with limited GWAS sample size, LD mismatch, and because the default five-component prior is more flexible than the true simulation (100 causal variants and rest SNPs with zero effect).


sbayesr.snpRes lists, for each SNP: ID, alleles, frequency, joint effect of A1, posterior SE, variance explained (i.e., per-SNP heritability estimate), PEP (posterior probability of the SNP with greater per-SNP heritability than a random SNP), mixture probabilities Pi1-5, and PIP (= 1 − Pi1). PIP, posterior inclusion probability, is useful for fine-mapping (probability of a non-zero effect given all other SNPs).

Q8: How many SNPs have high PIP (e.g. PIP > 0.8)? (you can read sbayesr.snpRes in R and check (sum(snpRes$PIP) > 0.8) Should you keep only those SNPs for prediction? Why or why not?

Show answer

You may see roughly tens of SNPs with PIP > 0.8 (many near the 100 causal variants, plus some false positives). For prediction, do NOT drop SNPs by PIP alone: SBayesR spreads signal across the genome via shrinkage, and PGS uses joint effects for all SNPs in snpRes. PIP is mainly for interpretation and fine-mapping, not for building the prediction equation.