R/stat.ranksum.R

In vdumeaux/mixtR: Matched Interaction across Tissues (MIxT) data analysis

#' Compute association between gene sets/modules across tissues
#'
#' Interactions between modules across tissues are identified using a random permutation approach
#' based on the correlation between ranksums of gene expression in gene sets/modules across tissues.
#'
#' @param mixt.ranksum output of sig.ranksum()
#' @param tissue1 name of the first tissue the test is performed for.
#' should be a valid name of mixt.ranksum
#' @param tissue2 name of second tissue the test is performed for.
#' should be a valid name of mixt.ranksum.
#' @param corType a character string indicating which correlation coefficient
#' is to be computed.
#' Default 'p' for pearson
#' @param nRuns number of permutations
#' @param randomSeed seed number for random number generation.
#' Default set as '1234'
#' @param mc.cores number of cores
#' @param verbose numerical. default > 0 show informational text on progress
#'
#' @export
#'
stat.ranksum <- function(mixt.ranksum,
                         tissue1,
                         tissue2,
                      corType = "p",
                      nRuns=10000, randomSeed = 12345,
                      mc.cores = 2,
                      verbose = 2)

{
  if (!tissue1 %in% names(mixt.ranksum))
    stop ("tissue1 should be a name of mixt.ranksum")
  if (!tissue2 %in% names(mixt.ranksum))
    stop ("tissue1 should be a name of mixt.ranksum")


  dat.ranksum <- lapply(mixt.ranksum, function (bs) {
    data.frame(lapply(bs, function(x) unlist(lapply(x, "[", "ranksum"))))
    })

  cohort.name <- names(mixt.ranksum[[1]][[1]])

  perm_cor_p <- lapply(cohort.name, function(cohort){
    dat.ranksum <- lapply(dat.ranksum, function (x) {x[grep(cohort, rownames(x)), ]})

    result = list()

    if (tissue1 != names(mixt.ranksum)[1]) {
      tissue.1 = tissue2
      tissue.2 = tissue1
    } else {
      tissue.1=tissue1
      tissue.2=tissue2
    }

    nSamples = length(dat.ranksum[[tissue.1]][[1]])

    seedSaved = FALSE

    if (!is.null(randomSeed))
      {
      if (exists(".Random.seed"))
      {
        seedSaved = TRUE;
        savedSeed = .Random.seed
      }
      set.seed(randomSeed)
    }

    mods <- parallel::mclapply(1:nRuns, function(i) {
    set.seed(randomSeed + 2*i + 1);
    if (verbose > 0) print(paste("...working on run", i, ".."));

    if (i > 1)
    {
      useTissue1Samples = sample(nSamples)
      useTissue2Samples = sample(nSamples)

    } else {
      useTissue1Samples = c(1:nSamples)
      useTissue2Samples = c(1:nSamples)}

      mods <- cor(dat.ranksum[[tissue.1]][useTissue1Samples, ],
                  dat.ranksum[[tissue.2]][useTissue2Samples, ], use=corType)

      rownames(mods) <- names(dat.ranksum[[tissue.1]])
      colnames(mods) <- names(dat.ranksum[[tissue.2]])

      return(mods)
  }, mc.cores = mc.cores)

    n.list <- NULL
    for (i in 2:length(mods)){
    n.list[[i-1]]<- abs(mods[[i]]) >= abs(mods[[1]])}

    ret <-  Reduce('+', n.list)/9999
    return(ret)})

    names(perm_cor_p) <- cohort.name
    return(perm_cor_p)

}