R/mts_clusterAvg.R

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

mts_clusterAvg <- function(exprsMatrix, crisprMatrix, cores) {
  
  # Input error checking
  if (missing(cores)) {
    cores <- 1
    message("1 core selected")
  } else if (!is.numeric(cores) || cores < 1) {
    stop("cores should be >= 1")
  } else {
    message(cores, " cores selected")
  }
  
  if (missing(exprsMatrix)) {
    stop("No Expression Matrix Provided: A gene by sample matrix is required to invoke this function")
  }
  
  if (missing(crisprMatrix)) {
    stop("No CRISPR Matrix Provided")
  }
  
  # filter rows with complete 0 values & ensure overlap between expression and crispr matrices
  exprsMatrix <- exprsMatrix[, which(colnames(exprsMatrix) %in% colnames(crisprMatrix))]
  exprsMatrix <- exprsMatrix[rowSums(exprsMatrix == 0) != ncol(exprsMatrix), ]
  crisprMatrix <- crisprMatrix[, which(colnames(crisprMatrix) %in% colnames(exprsMatrix))]
  
  # Generate result list in full format
  resultList1 <- pbmclapply(1:dim(exprsMatrix)[1], function(x) {
    
    # Appropriate error catching
    result <- tryCatch({
      # Find the best fitting number of clusters for each CCLE gene
      test.mog <- EM.findk(as.numeric(exprsMatrix[x, ]), model.types = "V", num.gaussians = 2:5)
      
      # Create mixture model values 
      m1 <- mog.density(as.numeric(exprsMatrix[x, ]), test.mog)
      
      # Max-min boundaries
      qMap2 <- c()
      num_clusters_G2 <- length(unique(m1$membership))
      for (i in 1:(num_clusters_G2 - 1)) {
        j <- i + 1
        max1 <- max(m1$x[m1$membership == i])
        min2 <- min(m1$x[m1$membership == j])
        qMap3 <- (min2 + max1) / 2 
        qMap2 <- c(qMap2, qMap3)  
      }
      
      # Resolve to table
      tab1 <- as.data.frame(cbind(colnames(exprsMatrix), as.numeric(exprsMatrix[x, ])), stringsAsFactors = FALSE)
      tab1[, 2] <- as.numeric(tab1[, 2])
      fullRange1 <- c(min(tab1[, 2] - 0.1), qMap2, max(tab1[, 2]))
      tab1$category1 <- cut(tab1[, 2], breaks = fullRange1, labels = 1:(length(fullRange1) - 1))
      
      un1 <- unique(as.numeric(tab1$category1))
      un1 <- un1[order(un1)]
      depList <- list(numeric(0))
      for (i in 1:length(un1))
      {
        cl1 <- tab1[which(tab1$category1 == un1[i]), 1]
        depSUB <- crisprMatrix[, which(colnames(crisprMatrix) %in% cl1)]
        meanDEP <- apply(depSUB, 1, function(x) mean(as.numeric(x[1:dim(depSUB)[2]]), na.rm = TRUE))
        depList[[i]] <- meanDEP
      }
      depTAB <- do.call("cbind", depList)
      colnames(depTAB) <- paste("Cluster", 1:dim(depTAB)[2], sep = " ")
      colnames(tab1) <- c("Sample", "Values", "Cluster_Assignment")
      oL1 <- list(data.frame(0))
      oL1[[1]] <- depTAB
      oL1[[2]] <- tab1
      return(oL1)
    }, error = function(err) {
      # If unable to generate clusters, return empty matrix
      oL1 <- list(data.frame(0))
      oL1 <- matrix(nrow = 1, ncol = 1)
      return(oL1)
    })
    
  }, mc.cores = cores)
  
  names(resultList1) <- rownames(exprsMatrix)
  return(resultList1)
}