Fast linear regression (single predictor with intercept) for GWAS

fast_gwas.R

run_gwas = function(X, y) {
  a_y = sd(y)
  b_x = apply(X, 2, sd)
  y_tilde = y / a_y
  # message(1)
  x_tilde = sweep(X, 2, FUN = '/', b_x)
  # message(2)
  xtx = colSums(x_tilde ^ 2)
  xbar = colMeans(x_tilde)
  ybar = mean(y_tilde)
  n = nrow(x_tilde)
  # message(dim(x_tilde), '  ', length(y_tilde))
  xty = colSums(sweep(x_tilde, 1, FUN = '*', y_tilde))
  # message(3)
  denom = xtx - n * xbar ^ 2
  mu0 = 1 / denom * (xtx * ybar - xbar * xty)
  bhat = 1 / denom * (-n * xbar * ybar + xty)
  ypred = sweep(sweep(x_tilde, 2, FUN = '*', bhat), 2, FUN = '+', mu0)
  # message(4)
  sigma2 = 1 / (n - 2) * colSums((sweep(ypred, 1, FUN = '-', y_tilde)) ^ 2)
  # message(5)
  bhat_se = sqrt(sigma2 / (xtx - n * xbar ^ 2))
  bhat = bhat * a_y / b_x
  bhat_se = bhat_se * a_y / b_x
  return(list(bhat = bhat, bhat_se = bhat_se))
}