R/mts_mixModelCluster_XPR.R

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

mts_mixModelCluster_XPR <- function(dataMatrix, 
                                    GeneXPRthresh = 3.321928, 
                                    NumSampleThresh = 20,
                                    cores) {
  
  # Input error checking
  if (missing(dataMatrix))        stop("No data matrix provided")
  if (!is.data.frame(dataMatrix)) stop("The dataMatrix must be a data frame")
  
  
  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(GeneXPRthresh))
  {
    GeneXPRthresh <- 3.321928
    message("Gene expression threshold not entered, using 3.321928")
  } else {
    message(paste("Gene expression threshold: ", GeneXPRthresh, sep=""))
  }
  
  
  if(missing(NumSampleThresh))
  {
    NumSampleThresh <- 20
    message("Number of samples threshold not entered, using 20 samples")
  } else {
    message(paste("Number of samples threshold: ", NumSampleThresh, sep=""))
  }
  
  # drop genes that are zero in every sample, then genes with NA rownames
  dataMatrix <- dataMatrix[rowSums(is.na(dataMatrix)) == 0, , drop = FALSE]
  dataMatrix = dataMatrix[rowSums(dataMatrix == 0) != ncol(dataMatrix), ]
  dataMatrix <- dataMatrix[!is.na(rownames(dataMatrix)), , drop = FALSE]
  
  # keep genes expressed above GeneXPRthresh in more than NumSampleThresh samples
  counts_over_thresh = rowSums(dataMatrix > GeneXPRthresh,
                               na.rm = TRUE)
  
  keepGenes = names(counts_over_thresh)[counts_over_thresh > NumSampleThresh]
  
  if (length(keepGenes) == 0) {
    stop(sprintf("No genes have an expression > %.6g in over %d samples",
                 GeneXPRthresh, NumSampleThresh))
  }
  
  exprsMatrixFiltered = dataMatrix[keepGenes,]
  
  
  # Generate result list in minimal format
  resultList1 <- pbmclapply(1:dim(exprsMatrixFiltered)[1], function(x) {
    
    # Appropriate error catching
    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
      tab1
    }, error = function(err) {
      matrix(nrow=1, ncol=2)
    }))
  }, mc.cores = cores)
  
  setNames(resultList1, rownames(exprsMatrixFiltered))
}