R/compare.sigs.R

In jhart99/gseasier: Run the GSEA java pacakage from inside R

# Copyright ---------------------------------------------------------------
# 2018 The Scripps Research Institute Author: Jonathan Ross Hart

# Author ------------------------------------------------------------------
# Jonathan Ross Hart(jonathan@jonathanrosshart.com)

# Description --------------------------------------------------------------

# A function to compare gene set overlap to determine if the genesets are driven
# by the same genes.


#' Compares sparse matrixes of gene signatures
#'
#' @param sig.matrix Signature sparse matrix
#' @param test.matrix Test signatures sparse matrix
#'
#' @return A data table of the comparisons of sig.matrix and test.matrix
#'   including odds ratios, p.values and FDR adjusted q.values
#' @examples
#' compare.sigs(sig.matrix, test.matrix)
compare.sigs <- function(sig.matrix, test.matrix) {
  # make the comparison
  system.time(results.vec <- sig.matrix %*% t(test.matrix))
  # get the other parts we need
  genes.per.test <- rowSums(test.matrix)
  genes.per.sig <- rowSums(sig.matrix)
  total.genes <- ncol(sig.matrix)
  gc()
  if ("pacakage:parallel" %in% search()) {
    system.time(p.values <- mcmapply(function(i, j) {
      phyper(results.vec[i, j], genes.per.sig[i],
             total.genes - genes.per.sig[i], genes.per.test[j], F)
    },
    i = rep(1:nrow(sig.matrix), times = length(genes.per.test)),
    j = rep(1:length(genes.per.test), each = nrow(sig.matrix)), mc.cores = 10))
  } else {
    system.time(p.values <- mapply(function(i, j) {
      phyper(results.vec[i, j], genes.per.sig[i],
             total.genes - genes.per.sig[i], genes.per.test[j], F)
    },
    i = rep(1:nrow(sig.matrix), times = length(genes.per.test)),
    j = rep(1:length(genes.per.test), each = nrow(sig.matrix))))
  }

  dim(p.values) <- c(nrow(sig.matrix), length(genes.per.test))
  rownames(p.values) <- row.names(sig.matrix)
  colnames(p.values) <- names(genes.per.test)
  odds.ratio <- total.genes * results.vec / genes.per.sig ^ 2

  output <- reshape2::melt(odds.ratio)
  output2 <- reshape2::melt(p.values)
  output$p <- output2$value
  colnames(output) <- c("sig", "test.sig", "odds.ratio", "p.value")
  output$q.value <- p.adjust(output$p.value, "BH")

  output
}