R/mts_Mutation.R

In MultiSEp: Predict Synthetic Lethality and Other Gene Dependency Relationships from Multiomics Data

mts_Mutation <- function(resultList=resultList, mutMatrix=mutMatrix, tissueMatrix=tissueMatrix, cores=cores, pVal=pVal)
{
  # Error catching - from input
  if(missing(cores))
  {
    cores <- 1
    message("1 core selected")
  } else {
    message(paste(cores, " cores selected", sep=""))
  }

  if(missing(pVal))
  {
    pVal <- 0.01
    message("P-Value Threshold not entered, using 0.01")
  } else {
    message(paste("P Value Threshold: ", pVal, sep=""))
  }

  if(missing(mutMatrix))
  {
    stop("No Mutation Matrix Provided")
  } else {
    if (!is.data.frame(mutMatrix)) {
      stop("mutMatrix must be a data frame")
      return(NULL)
    }
  }

  if(missing(resultList))
  {
    stop("No MultiSEp results list provided - please run mts_mixModelCluster() to generate")
  } else {
    if (!inherits(resultList, "list")) {
      stop("resultList must be a list of data frames - please run mts_mixModelCluster() to generate")
      return(NULL)
    }
  }

  if(missing(tissueMatrix))
  {
    message("No tissue matrix provided - calculating across all tissues")
    tissueMatrix <- as.data.frame(cbind(colnames(mutMatrix), "All"), stringsAsFactors=FALSE)
    colnames(tissueMatrix) <- c("cell_line", "tissue")
  } else {
    if (!is.data.frame(tissueMatrix)) {
      stop("tissueMatrix must be a data frame")
      return(NULL)
    }
  }
  
  colnames(tissueMatrix) <- c("cell_line", "tissue") # for when tissueMatrix has differing colnames

  per5 <- round(5*(dim(mutMatrix)[2] / 100))
  mC <- apply(mutMatrix, 1, function(x) length(x [ x != "WT" ] ))
  mutMat2 <- mutMatrix[which(mC > per5),] # 5%

  mutMat3 <- as.data.frame(t(mutMat2), stringsAsFactors=FALSE)
  mutMat3$cell_line <- rownames(mutMat3)
  lDF3 <- pbmclapply(1:length(resultList), function(z)
  {
    if(dim(resultList[[z]])[2] < 3)
    {
      return(c(names(resultList)[z], "NA", "NA", NA))
    }
    else
    {
      s1 <- resultList[[z]]
      colnames(s1) <- c("cell_line", "mRNA_Expression", "expression_mode")
      m2 <- merge(s1, tissueMatrix)
      m3 <- merge(m2, mutMat3)
      rm(m2)
      unD <- unique(m3$tissue)
      lDF2 <- as.data.frame(do.call("rbind", pblapply(1:length(unD), function(y)
      {
        s2 <- m3[which(m3$tissue == unD[y]),]
        unE <- as.numeric(unique(s2$expression_mode))
        lDF1 <- as.data.frame(do.call("rbind", lapply(5:dim(s2)[2], function(x) {
          contTab <- matrix(nrow=2, ncol=length(unique(s2$expression_mode)))
          for(k in 1:length(unE))
          {
            contTab[2,k] <- length(which(s2[which(s2[,3] == unE[k]),x] != "WT"))
            contTab[1,k] <- length(which(s2[which(s2[,3] == unE[k]),x] == "WT"))
          }
          return(c(names(resultList)[z], colnames(s2)[x], unD[y], suppressWarnings(chisq.test(contTab))$p.value, paste(contTab[2,], collapse = ",")))
        })), stringsAsFactors=FALSE)
        if(length(which(lDF1[,4] == "NaN")) > 0)
        {
          lDF1 <- lDF1[-which(lDF1[,4] == "NaN"),]
        } else {
          lDF1 <- lDF1
        }
        lDF1[,4] <- as.numeric(lDF1[,4])
        if(length(which(lDF1[,4] < pVal)) > 0)
        {
          return(lDF1[which(lDF1[,4] < pVal),])
        }
        else
        {
          # Do nothing
        }
    })), stringsAsFactors=FALSE)
    return(lDF2)
    }
  }, mc.cores=cores)
  outDF <- do.call("rbind", lDF3)
  if(length(which(outDF[,2] == "NA")) == 0)
  {
    outDF[,4] <- as.numeric(outDF[,4])
  } else {
    outDF <- outDF[-which(outDF[,2] == "NA"),]
    outDF[,4] <- as.numeric(outDF[,4])
  }
  outDF <- outDF[order(outDF[,4]),]
  colnames(outDF) <- c("mRNA_gene", "mutation_gene", "tissue", "chiSqPvalue", "mutation_per_mode")
  
  # Remove single mode associations
  if(length(which(unlist(lapply(strsplit(outDF[,5], split=","), length)) == 1)) > 0)
  {
    outDF <- outDF[-which(unlist(lapply(strsplit(outDF[,5], split=","), length)) == 1),]
  } else {
    outDF <- outDF
  }
  return(outDF)
}