R/helper_function.R

In ressomlab/INDEED: Interactive Visualization of Integrated Differential Expression and Differential Network Analysis for Biomarker Candidate Selection

## This file contains all the helper functions needed to 
## properly run non_partial_cor(), select_rho_partial(), partial_cor(), and network_display().

#' @title Compute the correlation
#' @description This function computes either the pearson or spearman correlation coefficient.
#' @param data_group_1 This is a n*p matrix.
#' @param data_group_2 This is a n*p matrix.
#' @param type_of_cor If this is NULL, pearson correlation coefficient will be calculated as 
#'     default. Otherwise, a character string "spearman" will calculate the spearman correlation
#'     coefficient.
#' @return A list of correlation matrices for both group 1 and group 2.
#' @importFrom stats cor 

compute_cor <- function(data_group_1, data_group_2, type_of_cor) {
    if (is.null(type_of_cor) || type_of_cor == "pearson") {
        cor_group_1 <- cor(data_group_1, method = "pearson")
        cor_group_2 <- cor(data_group_2, method = "pearson")

    } else if (type_of_cor == "spearman") {
        cor_group_1 <- cor(data_group_1, method = "spearman")
        cor_group_2 <- cor(data_group_2, method = "spearman")
    }
    cor <- list("Group1" = cor_group_1, "Group2" = cor_group_2)
}


#' @title Compute the partial correlation
#' @description This function computes the partial correlation coefficient.
#' @param pre_inv This is an inverse covariance matrix.
#' @return A \eqn{p*p} partial correlation matrix.
#' @importFrom utils tail

compute_par <- function(pre_inv) {
  p <- nrow(pre_inv)

  i <- rep(seq_len(p-1), times=(p-1):1)
  k <- unlist(lapply(2:p, seq, p))

  pre_inv_i <- vapply(seq_len(p-1), function(x) pre_inv[x,x], numeric(1))
  pre_inv_i <- rep(pre_inv_i, times=(p-1):1)

  pre_inv_j <- vapply(2:p, function(x) pre_inv[x,x], numeric(1))
  pre_inv_j <- unlist(lapply(seq_len(p), function(x) tail(seq_len(p), -(x))))

  pc_value <- pre_inv[upper.tri(pre_inv)]
  pc_calc <- -pc_value / sqrt(pre_inv_i * pre_inv_j)

  pc <- matrix(0, p, p)
  pc[upper.tri(pc)] <- pc_calc
  pc[lower.tri(pc)] <- t(pc)[lower.tri(t(pc))]
  return(pc)
}


#' @title Permutations to build a differential network based on correlation analysis
#' @description A permutation test that randomly permutes the sample labels in distinct
#'     biological groups for each biomolecule. The difference in each paired biomolecule
#'     is considered statistically significant if it falls into the 2.5% tails on either end of the 
#'     empirical distribution curve. 
#' @param m This is the number of permutations.
#' @param p This is the number of biomarker candidates.
#' @param n_group_1 This is the number of subjects in group 1.
#' @param n_group_2 This is the number of subjects in group 2.
#' @param data_group_1 This is a \eqn{n*p} matrix containing group 1 data.
#' @param data_group_2 THis is a \eqn{n*p} matrix containing group 2 data.
#' @param type_of_cor If this is NULL, pearson correlation coefficient will be calculated as 
#'     default. Otherwise, a character string "spearman" will calculate the spearman correlation 
#'     coefficient.
#' @return A multi-dimensional matrix that contains the permutation result.
#' @importFrom utils txtProgressBar setTxtProgressBar
#' @importFrom stats cor 

permutation_cor <- function(m, p, n_group_1, n_group_2, data_group_1, data_group_2, type_of_cor) {
    diff_p <- array(0, dim = c(m, p, p))
    pb <- txtProgressBar(min = 0, max = m, style = 3)
    for (t in 1 : m) {
        data_group_1_p <- matrix(0, n_group_1, p)
        for (i in 1 : p) {
            data_group_1_p[, i] <- data_group_1[sample(n_group_1), i]
        }
        data_group_2_p <- matrix(0, n_group_2, p)
        for (i in 1 : p) {
            data_group_2_p[, i] <- data_group_2[sample(n_group_2), i]
        }

    if (is.null(type_of_cor)) {
        cor_group_1_p <- cor(data_group_1_p, method = "pearson")
        cor_group_2_p <- cor(data_group_2_p, method = "pearson")
    } else {
        cor_group_1_p <- cor(data_group_1_p, method = "spearman")
        cor_group_2_p <- cor(data_group_2_p, method = "spearman")
    }
        diff_p[t, , ] <- cor_group_2_p - cor_group_1_p

        # update progress bar
        setTxtProgressBar(pb, t)
    }
    close(pb)
    return(diff_p)
}


#' @title Permutations to build differential network based on partial correlation analysis
#' @description A permutation test that randomly permutes the sample labels in distinct
#'     biological groups for each biomolecule. The difference in paired partial correlation
#'     is considered statistically significant if it falls into the 2.5% tails on either end of the 
#'     empirical distribution curve.  
#' @param m This is the number of permutations.
#' @param p This is the number of biomarker candidates.
#' @param n_group_1 This is the number of subjects in group 1.
#' @param n_group_2 This is the number of subjects in group 2.
#' @param data_group_1 This is a \eqn{n*p} matrix containing group 1 data.
#' @param data_group_2 This is a \eqn{n*p} matrix containing group 2 data.
#' @param rho_group_1_opt This is an optimal tuning parameter to obtain a sparse differential 
#'     network for group 1.
#' @param rho_group_2_opt This is an optimal tuning parameter to obtain a sparse differential
#'     network for group 2.
#' @return A multi-dimensional matrix that contains the permutation result.
#' @importFrom utils txtProgressBar setTxtProgressBar
#' @importFrom glasso glasso

permutation_pc <- function(m, p, n_group_1, n_group_2, data_group_1, data_group_2, rho_group_1_opt, 
                           rho_group_2_opt) {
    diff_p <- array(0, dim = c(m, p, p))
    pb <- txtProgressBar(min = 0, max = m, style = 3)
    for(t in 1 : m) {
        data_group_1_p <- matrix(0, n_group_1, p)
        for(i in 1 : p) {
            data_group_1_p[, i] <- data_group_1[sample(n_group_1), i]
        }
        data_group_2_p <- matrix(0, n_group_2, p)
        for(i in 1 : p) {
            data_group_2_p[, i] <- data_group_2[sample(n_group_2), i]
        }
        per_group_1 <- glasso(var(data_group_1_p), rho = rho_group_1_opt)
        per_group_2 <- glasso(var(data_group_2_p), rho = rho_group_2_opt)
        pc_group_1_p <- compute_par(per_group_1$wi)
        pc_group_2_p <- compute_par(per_group_2$wi)
        diff_p[t, , ] <- pc_group_2_p - pc_group_1_p
        # update progress bar
        setTxtProgressBar(pb, t)
    }
    close(pb)
    return(diff_p)
}


#' @title Calculate the positive and negative thresholds based on the permutation result
#' @description This function calculates the positive and negative thresholds based on the 
#'     permutation result.
#' @param thres_left This is the threshold representing 2.5 percent of the left tail of the 
#'     empirical distributuion curve.
#' @param thres_right This is the threshold representing 2.5 percent of the right tail of the 
#'     empirical distributuion curve.
#' @param p This is the number of biomarker candidates.
#' @param diff_p This is the permutation result from either permutation_cor or permutation_pc.
#' @return A list of positive and negative thresholds.
#' @importFrom stats quantile

permutation_thres <- function(thres_left, thres_right, p, diff_p) {
    significant_thres_p <- matrix(0, p, p)
    significant_thres_n <- matrix(0, p, p)
    for (i in 1 : (p-1)) {
        for (j in (i + 1) : p) {
            significant_thres_n[i, j] <- quantile(diff_p[, i, j], probs = thres_left)
            significant_thres_n[j, i] <- significant_thres_n[i, j]
            significant_thres_p[i, j] <- quantile(diff_p[, i, j], probs = thres_right)
            significant_thres_p[j, i] <- significant_thres_p[i, j]
        }
    }
    significant_thres <- list("positive" = significant_thres_p, "negative" = significant_thres_n)
    return(significant_thres)
}


#' @title Calculate the differential network score
#' @description This function calculates differential network score by using the binary link and 
#'     z-scores.
#' @param binary_link This is the binary correlation matrix with 1 indicating positive correlation 
#'     and -1 indicating negative correlation for each biomolecular pair.
#' @param z_score This is converted from the given or calculated p-value.
#' @return An activity score associated with each biomarker candidate.

compute_dns <- function(binary_link, z_score) {
    # get adjacent matrix
    diff_d <- abs(binary_link)
    # set diagonal elements to 1
    diag(diff_d) <- 1
    # compute differential network score for each row
    dns <- apply(diff_d, 1, function(x, y = z_score) sum(y[which(x == 1)]))
    return(dns)
}


#' @title Obtain p-values using logistic regression
#' @description This function calculates p-values using logistic regression in cases that p-values 
#'     are not provided.
#' @param x This is a data frame consists of data from group 1 and group 2.
#' @param class_label This is a binary array indicating 0 for group 1 and 1 for group 2.
#' @param Met_name This is an array of IDs.
#' @return p-values
#' @importFrom stats glm

pvalue_logit <- function(x, class_label, Met_name) {
    data_tp <- as.data.frame(t(x))    # n*p
    class_label_tp <- as.data.frame(t(class_label))
    pvalue <- c()
    # attach metabolites ID and class label in the data set
    X_df <- cbind(data_tp, class_label_tp)
    colnames(X_df)[1:(ncol(X_df)-1)] <- Met_name
    colnames(X_df)[ncol(X_df)] <- "Class"
    for (i in 1:(ncol(X_df)-1)) {
      X_df_tempt <- X_df[,c(i, ncol(X_df))]
      glm.fit <- glm(Class ~. , family = "binomial", data = X_df_tempt)
      pvalue_tempt <- summary(glm.fit)$coefficients[,4][2]
      pvalue <- c(pvalue, pvalue_tempt)
    }
    pvalue_df <- data.frame("ID" = Met_name, "p.value" = pvalue)
    return(pvalue_df)
}


#' @title Create log likelihood error
#' @description This function calculates the log likelihood error. 
#' @param data This is a matrix.
#' @param theta This is a precision matrix.
#' @return log likelihood error 

loglik_ave <- function(data, theta){
    loglik <- c()
    loglik <- log(det(theta)) - sum(diag(var(data) %*% theta))
    return(-loglik)
}


#' @title Draw error curve
#' @description This function draws error curve using cross-validation.
#' @param data This is a matrix.
#' @param n_fold This parameter specifies the n number in n-fold cross_validation.
#' @param rho This is the regularization parameter values to be evalueated in terms their errors.
#' @return A list of errors and their corresponding \eqn{log(rho)}.
#' @importFrom glasso glasso

choose_rho <- function(data, n_fold, rho) {
  # randomly shuffle the data
  Data <- data[sample(nrow(data)), ]
  # create n_fold equally size folds
  folds <- cut(seq(1, nrow(Data)), breaks = n_fold, labels = FALSE)
  # tune parameters
  d <- ncol(Data)

  loglik <- lapply(seq_along(rho), function(i) {
    vapply(seq_len(n_fold), function(j) {
      # segement your data by fold using the which() function
      testIndexes <- which(folds == j, arr.ind = TRUE)
      testData <- Data[testIndexes, ]
      trainData <- Data[-testIndexes, ]
      # use test and train data partitions however you desire...
      cov <- var(trainData) # compute the covariance matrix
      pre <- glasso(cov, rho = rho[i])
      loglik_ave(testData, pre$wi)
    }, numeric(1))})

  loglik_cv <- vapply(loglik, mean, numeric(1))
  loglik_rho <- vapply(loglik, function(x) sd(x) / sqrt(n_fold), numeric(1))

  #plot(rho, loglik_cv, xlab = expression(lambda), ylab = "Error")
  #lines(rho, loglik_cv)
  error <- list("log.cv" = loglik_cv, "log.rho" = loglik_rho)
  return(error)
}        


#' @title Compute p-value for edges
#' @description This function computes p-value for edges based on permutation result. 
#' @param p This is the number of biomarker candidates.
#' @param diff This is the delta correlation or partial correlation matrix.
#' @param diff_p This is the permutation result from either permutation_cor or permutation_pc.
#' @param m This is the number of permutations.
#' @return p-value for edges.

compute_pvalue_edge <- function(p, diff, diff_p, m) {
  significant_thres <- matrix(0, p, p)
  for (i in 1 : (p-1)) {
    for (j in (i + 1) : p) {
      significant_thres[i, j] <- min(length(diff_p[,i,j][diff_p[,i,j]>=diff[i,j]]),
                                     length(diff_p[,i,j][diff_p[,i,j]<=diff[i,j]])) 
      significant_thres[j, i] <- significant_thres[i, j]
    }
  }
  pvalue_edge <- significant_thres/m
  # adjust for two sides
  pvalue_edge[pvalue_edge <= 0.5] <- 2 * pvalue_edge[pvalue_edge <= 0.5]
  
  diag(pvalue_edge) <- 1
  return(pvalue_edge)
}


#' @title Compute fdr p-value for edges
#' @description This function computes fdr p-value for edges to adjust for multiple testing.
#' @param p This is the number of biomarker candidates.
#' @param pvalue_edge This is p-value for edges from compute_pvalue_edge.
#' @return Adjusted p-value for edges by fdr.
#' @importFrom stats p.adjust

compute_pvalue_edge_fdr <- function(p, pvalue_edge) {
  pvalue_edge_vector <- vector()
  for (i in 1:(p-1)){
    for (j in (i+1):p){
      pvalue_edge_vector = append(pvalue_edge_vector, c(i,j,pvalue_edge[i,j]))
    }
  }
  pvalue_edge_vector <- matrix(pvalue_edge_vector, ncol = 3, byrow = T)
  pvalue_edge_vector_fdr <- p.adjust(pvalue_edge_vector[,3], method = "fdr", n = length(pvalue_edge_vector[,3]))
  pvalue_edge_fdr <- matrix(0, p, p)
  num <- 1
  for (i in 1 : (p-1)) {
    for (j in (i + 1) : p) {
      pvalue_edge_fdr[i, j] <- pvalue_edge_vector_fdr[num]
      pvalue_edge_fdr[j, i] <- pvalue_edge_fdr[i, j]
      num <- num + 1
    }
  }
  diag(pvalue_edge_fdr) = 1
  return(pvalue_edge_fdr)
}


#' @title Compute edge weights
#' @description This function computes edge weights based on p-value for edges with directions.
#' @param pvalue_edge_fdr This is the p-value for edges possibly after multiple testing correction.
#' @param binary_link This is the binary edge connection.
#' @return Edge weights.
#' @importFrom stats qnorm

compute_edge_weights <- function(pvalue_edge_fdr, binary_link) {
  zscore_edge_fdr <- abs(qnorm(1 - pvalue_edge_fdr/2))
  # 1.5 is a predefined factor to cap zero-pvalue connection
  inf_cap <- 1.5 * max(zscore_edge_fdr[is.finite(zscore_edge_fdr)]) 
  zscore_edge_fdr[is.infinite(zscore_edge_fdr)] <- inf_cap
  weight_link <- zscore_edge_fdr * binary_link
  return(weight_link)
}