R/run_silence.R

In tgrimes/dnapath2: Differential Network Analysis Using Gene Pathways

#' Wrapper for silencer method
#' 
#' See Barzel and Barabasi (2013) for details.
#' @param x The n by p matrix of counts.
#' @param threshold Cutoff for significant associations. If NULL, all correlations
#' are returned. Otherwise, correlations of magnitude at or below this threshold are 
#' set to zero.
#' @param verbose If true, update messages are printed.
#' @return A p by p matrix of association scores.
#' @export
run_silencer <- function(x, method = "spearman", verbose = FALSE, ...) {
  Norm <-  0.5
  # First calculate correlation matrix.
  if(method == "pcor") {
    G <- corpcor::pcor.shrink(x)
    G <- matrix(as.numeric(G), nrow(G), ncol(G))
  } else {
    G <- cor(x, method = method)
  }
  # diag(G) <- 0
  
  # Preprocessing: check singularity of G.
  N <- nrow(G)
  R <- qr(G)$rank
  
  alpha <- 1
  
  if(R < N) {
    G <- G * 0.5 # Initial Norm = 0.5.
    diag(G) <- 1
    alpha <- alpha * Norm
  }
  
  MaxValue <- 2
  Maximum <- 0.9
  Minimum <- 0.5
  
  Go <- G - diag(N)
  S <- (Go + diag(diag(Go %*% G))) %*% solve(G)
  S <- Go
  
  n <- 0
  while((MaxValue <= Minimum || MaxValue >= Maximum) && n < 20) {
    Go <- G - diag(N)
    
    # Perform the transormation.
    S <- (Go + diag(diag(Go %*% G))) %*% solve(G)
    MaxValue <- max(abs(eigen(S)$values))
    
    if(MaxValue <= Minimum || MaxValue >= Maximum) {
      Norm <- 1 / (sqrt(MaxValue / 0.5 * (Maximum - Minimum)))
      # If Norm is close to 1, calculate again.
      if(abs(1 - Norm) < 0.01) {
        Norm <- 1 / (sqrt(MaxValue / 0.5 * (Maximum + Minimum)))
      }
      G <- G * Norm
      diag(G) <- 1
      alpha <- alpha * Norm
    }
    
    n <- n + 1
    
    if(verbose) {
      cat(paste("n =", n, 
                  "- alpha =", round(alpha, 4), 
                  "MaxValue =", round(MaxValue, 4),
                  "Norm =", round(Norm, 4), "\n"))
    }
  }
  
  # S <- cov2cor(S)
  diag(S) <- 0
  S <- S / max(abs(S))
  S <- 0.5 * (S + t(S))
  
  colnames(S) <- colnames(x)
  
  return(S)
}