R/run_glasso.R

In tgrimes/dnapath2: Differential Network Analysis Using Gene Pathways

# Required:
#   edgeR - for calcNormFactors() and cpm().
#   glasso - for adjacency.fromSimilarity()
# library(edgeR)
# library(huge)

#' Wrapper for glasso method
#' 
#' Conducts co-expression analysis using glasso.
#' @param x The n by p matrix of counts.
#' @param method Graph estimation method. Should be either "mb", "glasso", or "ct".
#' Default is "glasso".
#' @param criterion Model selection criterion. Should be either "ric" or "stars".
#' Default is "ric".
#' @param verbose boolean indicating whether trace information during estimation 
#' and selection should be printed. Default is FALSE.
#' @return A p by p matrix of association scores.
#' @export
run_glasso <- function(x, 
                       method = c("glasso", "mb", "ct"), 
                       criterion = c("ric", "stars"), 
                       verbose = FALSE, ...) {
  method <- tolower(method[1])
  criterion <- tolower(criterion[1])
  if(!(method %in% c("mb", "glasso", "ct"))) {
    stop('method should be one of "mb", "glasso", or "ct".')
  }
  if(criterion != "ric" & criterion != "stars") {
    stop("criterion must be either \"ric\" or \"stars\".")
  }
  
  # Use unsigned to treat positive and negative associations equally.
  x_huge <- huge::huge(x, method = method, verbose = verbose)
  result <- huge::huge.select(x_huge, criterion = criterion, verbose = verbose)
  
  # Choose largest lambda at or below the optimal value.
  index <- which(result$lambda <= result$opt.lambda)[1]
  if(length(index) == 0) {
    index <- 10 # = length(x_huge$icov)
  } else if(is.na(index)) {
    # If no lambdas were smaller than optimal value, choose the smallest.
    index <- 1
  }
  
  if(verbose) {
    cat("Optimal lambda:", result$opt.lambda, "\n")
  }
  
  scores <- as.matrix(x_huge$icov[[index]])
  scores <- cov2cor(scores)
  diag(scores) <- 0
  scores <- (scores + t(scores)) / 2
  
  colnames(scores) <- colnames(x)
  
  return(scores)
}