R/mts_pairwiseExpressionPlot.R

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

mts_pairwiseExpressionPlot <- function(exprsMatrix, gene1, gene2)
{
  if(missing(exprsMatrix))
  {
    stop("No Expression Matrix Provided")
  } else {
    if (!is.data.frame(exprsMatrix)) {
      stop("exprsMatrix must be a data frame")
    }
    if(missing(gene1))
    {
      stop("No Gene1 provided")
    }
    
    if(missing(gene2))
    {
      stop("No Gene2 provided")
    }
    
    if(length(which(rownames(exprsMatrix) %in% gene1)) == 0 | length(which(rownames(exprsMatrix) %in% gene2)) == 0)
    {
      stop("Queried genes are not available in gene expression matrix")
    }
    # restore the caller's graphics parameters on exit 
    oldpar <- par("mar")         
    on.exit(par(mar = oldpar))   
    
    exprsResult <- tryCatch({
      # Find the best fitting number of clusters for each CCLE gene
      test.mog2 <- EM.findk(as.numeric(exprsMatrix[gene2,]),model.types="V", num.gaussians=2:5)
      test.mog <- EM.findk(as.numeric(exprsMatrix[gene1,]),model.types="V", num.gaussians=2:5)
      
      m2 <- mog.density(as.numeric(exprsMatrix[gene2,]), test.mog2)
      m1 <- mog.density(as.numeric(exprsMatrix[gene1,]), test.mog)
      # the plot 
      par(mar=c(5,5,5,5))
      plot(as.numeric(exprsMatrix[gene1,]), as.numeric(exprsMatrix[gene2,]), pch=21, cex=2,
           xlab=paste(gene1, " log2 gene expression"), ylab=paste(gene2, " log2 gene expression"), bg="lightblue", col="red",
           cex.axis=2, cex.lab=2)
      
      # gene 1 / m1 
      num_clusters_G1 <- length(unique(m1$membership))
      # iterate cluster 1,2,3,4 etc, then calculate min and max values, j + 1 for next cluster up.
      for (i in 1:(num_clusters_G1 - 1)) {
        j = i + 1
        # min and max of each cluster 
        max1 <- max(m1$x[m1$membership==i])
        min2 <- min(m1$x[m1$membership==j])
        
        # maximum expression of cluster(i) and minimum of clusters(i+1)
        ab_valG1 <- (min2 + max1) / 2
        # abline for the plot 
        abline(v=ab_valG1, lwd=2, col="blue") # x-axis, add abline for each ab_valG1 
      }
      
      
      # gene 2 / m2 
      num_clusters_G2 <- length(unique(m2$membership))
      # iterate cluster 1,2,3,4 etc, then calculate min and max values, j + 1 for next cluster up.
      for (i in 1:(num_clusters_G2 - 1)) {
        j <- i + 1
        # min and max of each cluster 
        max1 <- max(m2$x[m2$membership==i])
        min2 <- min(m2$x[m2$membership==j])
        
        # maximum expression of cluster(i) and minimum of clusters(i+1)
        ab_valG2 <- (min2 + max1) / 2 
        # abline for the plot 
        abline(h=ab_valG2, lwd=2, col="blue") # y-axis 
      }
      
    }, error = function(err) {
      warning("Cannot produce models for specified genes")
      NULL
    })
    return(invisible(exprsResult))
  }
}