R/groupFisher.R

In YunhuiQi/TestPMD: Post Selection Inference on Results from Penalized Matrix Decomposition (SCCA only)

#' Fisher Exact Test for group enrichment analysis of significant loadings
#'
#' This function generates the heatmap of qvalues from Fisher Exact Test of group enrichment analysis of sinigicant loadings.
#' @param data Data matrix corresponding to the canonical loadings, samples are in rows, variables are in columns.
#' @param loadings A named vector containing the estimated loadings.
#' @param qvalues A vector containing qvalues for the canonical loadings.
#' @param feature_groups A vector indicating the group information of variables.
#' @param alpha Significance level, 0.05 by default.
#' @param draw Whether draw the heatmap, default is TRUE.
#' @import ggplot2
#' @import grDevices
#' @import RColorBrewer
#' @importFrom stats fisher.test
#' @importFrom stats p.adjust
#' @export
#' @return A vector containing the qvalues of fisher exact tests
#' @examples
#' library(TestPMD)
#' data("covid")
#' out <- PMA::CCA(standsdmu(covid$metabolite), standsdmu(covid$protein),
#' typex = "standard", typez="standard",
#' penaltyx = 0.9, penaltyz = 0.9, K = 3, standardize = FALSE, trace = FALSE)
#' pvalue <- CPTloading(X = covid$metabolite, Y = covid$protein, side = "X", K = 1, r = 10,
#' penalty = "Fixed", rho_x = 0.9, rho_y = 0.9, permutation_no = 100)
#' feature_groups <- c(rep("A", 61), rep("B", 50))
#' groupFisher(data = covid$metabolite, loadings = out$u[,1], qvalues = pvalue,
#' feature_groups = feature_groups, alpha = 0.1, draw = FALSE)


groupFisher <- function(data, loadings, qvalues, feature_groups, alpha, draw = T){
  if(!is.null(feature_groups) & ncol(data) != length(feature_groups)){stop("Feature groups do not match data")}
  if(ncol(data) != length(loadings)){stop("Loadings do not match data")}
  if(length(qvalues) != length(loadings)){stop("qvalues do not match loadings")}
  if(!is.numeric(alpha)){stop("alpha must be numeric")}
  if(alpha >= 1 | alpha <= 0){stop("alpha must be within (0, 1)")}
  log_qvalue <- NULL
  qvalue <- NULL
  fisher_pvalue=c()
  n_variables <- length(feature_groups)
  n_groups <- length(unique(feature_groups))
  for (group_ind in 1:n_groups) {
    group <- unique(feature_groups)[group_ind]
    t11 <- sum((qvalues < alpha) & (feature_groups == group))
    t12 <- sum((qvalues >= alpha) & (feature_groups == group))
    t21 <- sum((qvalues < alpha) & (feature_groups != group))
    t22 <- sum((qvalues >= alpha) & (feature_groups != group))
    counttable <- matrix(c(t11, t12, t21, t22), nrow = 2, byrow = T)
    fisher_pvalue[group_ind] <- stats::fisher.test(counttable, alternative = "greater")$p.value
  }
  fisher_qvalue <- stats::p.adjust(fisher_pvalue, method="fdr")
  names(fisher_qvalue) <- unique(feature_groups)
  n_cols <- length(unique(feature_groups))
  mycolors <- grDevices::colorRampPalette(RColorBrewer::brewer.pal(8, "Set3"))(n_cols)
  fisher_dat <- data.frame(log_qvalue = log(fisher_qvalue), feature_groups = unique(feature_groups), qvalue = fisher_qvalue)
  suppressWarnings(
    p_fisher <- ggplot2::ggplot(data = fisher_dat, ggplot2::aes(x = factor(1) , y = feature_groups, fill = log_qvalue),) +
    ggplot2::geom_tile(show.legend = T)+
    ggplot2::scale_fill_gradient(low = "red", high = "lightgreen", limits=c(-10,0))+
    ggplot2::theme(plot.title  = ggplot2::element_text(hjust = 0.5))+
    ggplot2::geom_text(ggplot2::aes(factor(1), feature_groups, label = round(qvalue, 4)), color = "black", size = 4) +
    ggplot2::labs(fill = "ln(q-value)")+
    ggplot2::labs(y = "Groups", x = "",title = "Feature groups enrichment analysis")+
    ggplot2::theme(plot.title = element_text(hjust = 0.5))+
    ggplot2::theme_linedraw()+
    ggplot2::theme(axis.text.y = ggplot2::element_text(angle = 0, hjust = 1, color = mycolors, face = "bold", size = 12))+
    ggplot2::theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
                   panel.background = element_blank(), axis.line = element_line(colour = "black")))
  print(fisher_qvalue)
  if(draw){return(p_fisher)}
  else{return(fisher_qvalue)}

}