R/mts_mixModelCluster.R

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

mts_mixModelCluster <- function(dataMatrix, 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(dataMatrix)) {
    stop("No data matrix provided")
  }
  
  if (!is.data.frame(dataMatrix)) {
    stop("The dataMatrix must be a data frame")
  }
  # filter rows with complete 0 values
  exprsMatrixFiltered <- dataMatrix[rowSums(dataMatrix == 0) != ncol(dataMatrix), ]
  
  # Generate result list in minimal format
  resultList1 <- pbmclapply(1:dim(exprsMatrixFiltered)[1], function(x) {
    
    # Appropriate error catching
    result <- suppressWarnings(tryCatch({
      
      # Find the best fitting number of clusters for each CCLE gene
      test.mog <- EM.findk(as.numeric(exprsMatrixFiltered[x, ]), model.types = "V", num.gaussians = 2:5)
      
      # Find the mixture model values
      m1 <- mog.density(as.numeric(exprsMatrixFiltered[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(exprsMatrixFiltered), as.numeric(exprsMatrixFiltered[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))
      colnames(tab1) <- c("Sample", "Values", "Cluster_Assignment")
      # Generate object to catch all results for meta analysis
      return(tab1)
    }, error = function(err) {
      oL1 <- matrix(nrow=1, ncol=2)
      return(oL1)
    }))
    
  }, mc.cores = cores)
  
  resultList1_named <- setNames(resultList1, rownames(exprsMatrixFiltered))
  return(resultList1_named)
}