PRS-CS: a polygenic prediction method that infers posterior SNP effect sizes under continuous shrinkage (CS) priors

This function is a wrapper for the PRS-CS method implemented in C++. It takes marginal effect size estimates from regression and an external LD reference panel and infers posterior SNP effect sizes using Bayesian regression with continuous shrinkage priors.

Usage

prsCs(
  bhat,
  R,
  n,
  a = 1,
  b = 0.5,
  phi = NULL,
  maf = NULL,
  nIter = 1000,
  nBurnin = 500,
  thin = 5,
  verbose = FALSE,
  seed = NULL
)

Arguments

bhat

A vector of marginal effect sizes.

R

The LD correlation matrix (a single matrix over the analysed window), as in susieR::susie_rss().

n

Sample size of the GWAS.

a

Shape parameter for the prior distribution of psi. Default is 1.

b

Scale parameter for the prior distribution of psi. Default is 0.5.

phi

Global shrinkage parameter. If NULL, it will be estimated automatically. Default is NULL.

maf

A vector of minor allele frequencies, if available, will standardize the effect sizes by MAF. Default is NULL.

nIter

Number of MCMC iterations. Default is 1000.

nBurnin

Number of burn-in iterations. Default is 500.

thin

Thinning factor for MCMC. Default is 5.

verbose

Whether to print verbose output. Default is FALSE.

seed

Random seed for reproducibility. Default is NULL.

Value

A list containing the posterior estimates: - betaEst: Posterior estimates of SNP effect sizes. - psiEst: Posterior estimates of psi (shrinkage parameters). - sigmaEst: Posterior estimate of the residual variance. - phiEst: Posterior estimate of the global shrinkage parameter.

Examples

# Generate example data
set.seed(985115)
n <- 350
p <- 16
sigmasq_error <- 0.5
zeroes <- rbinom(p, 1, 0.6)
beta.true <- rnorm(p, 1, sd = 4)
beta.true[zeroes] <- 0

X <- cbind(matrix(rnorm(n * p), nrow = n))
X <- scale(X, center = TRUE, scale = FALSE)
y <- X %*% matrix(beta.true, ncol = 1) + rnorm(n, 0, sqrt(sigmasq_error))
y <- scale(y, center = TRUE, scale = FALSE)

# Calculate sufficient statistics
XtX <- t(X) %*% X
Xty <- t(X) %*% y
yty <- t(y) %*% y

# Set the prior
K <- 9
sigma0 <- c(0.001, .1, .5, 1, 5, 10, 20, 30, .005)
omega0 <- rep(1 / K, K)

# Calculate summary statistics
b.hat <- sapply(1:p, function(j) {
  summary(lm(y ~ X[, j]))$coefficients[-1, 1]
})
s.hat <- sapply(1:p, function(j) {
  summary(lm(y ~ X[, j]))$coefficients[-1, 2]
})
R.hat <- cor(X)
var_y <- var(y)
sigmasq_init <- 1.5

# Run PRS CS
maf <- rep(0.5, length(b.hat)) # fake MAF
out <- prsCs(b.hat, R.hat, n, maf = maf)
# In sample prediction correlations
cor(X %*% out$betaEst, y) # 0.9944553
#>           [,1]
#> [1,] 0.9936077