get_sparse_pca.R

get_sparse_pca <- function(gbm, n_pcs = 10, fastpath = FALSE) {
  x <- as.matrix(t(Biobase::exprs(gbm)))
  mu <- colMeans(x)
  s <- apply(x, 2, sd, na.rm = TRUE)
  s[s == 0] <- min(s[s > 0])
  svd_res <- irlba::irlba(x, nv = n_pcs, center = mu, scale = s, fastpath = fastpath)
  n <- dim(x)[1]
  variance_sum <- sum(apply(x, 2, var, na.rm = TRUE)/(s^2))
  var_pcs <- svd_res$d^2/(n - 1)/variance_sum
  pca_with_outliers <- list(
    x = svd_res$u %*% diag(svd_res$d), 
    rotation = svd_res$v,
    sdev = svd_res$d/sqrt(n - 1), 
    tot_var = variance_sum,
    var_pcs = var_pcs
  )
  pca_with_outliers
}