R/mts_CrisprTS.R

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

# MultiSEp CRISPR function: Tissue specific
mts_CrisprTS <- function(resultList=resultList, crisprMatrix=crisprMatrix, fcVal=fcVal, pVal=pVal, cores=cores, tissueMatrix=tissueMatrix, allDisRes=allDisRes)
{
  # Error catching - from input
  if(missing(cores))
  {
    cores <- 1
    message("1 core selected")
  } else {
    message(paste(cores, " cores selected", sep=""))
  }
  
  if(missing(fcVal))
  {
    fcVal <- -0.1
    message("Fold Change Threshold not entered, using -0.1")
  } else {
    message(paste("Fold Change Threshold: ", fcVal, sep=""))
  }
  
  if(missing(pVal))
  {
    pVal <- 0.1
    message("P-Value Threshold not entered, using 0.1")
  } else {
    message(paste("P Value Threshold: ", pVal, sep=""))
  }
  
  if(missing(crisprMatrix))
  {
    stop("No CRISPR Matrix Provided")
  } else {
    if (!is.data.frame(crisprMatrix)) {
      stop("crisprMatrix must be a data frame")
      return(NULL)
    }
  }
  
  if(missing(resultList))
  {
    stop("No MultiSEp results list provided - please run mts_clusterAvg to generate")
  } else {
    if (!inherits(resultList, "list")) {
      stop("resultList must be a list of data frames")
      return(NULL)
    }
  }
  
  if(missing(tissueMatrix))
  {
    stop("No tissue matrix provided")
  } 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
  
  
  if(missing(allDisRes))
  {
    stop("No multisep-crispr result table provided - please run mts_Crispr to generate")
  } else {
    if (!is.data.frame(allDisRes)) {
      stop("allDisRes must be a data frame")
      return(NULL)
    }
    dbTab1 <- allDisRes
  }
  
  listGroupD <- pbmclapply(seq_len(nrow(dbTab1)), function(i) {
    if (nrow(dbTab1)==0) {
      stop("No results found by mts_CrisprTS() with the current threshold values")
      return(NULL)
    }
    if(dbTab1[i,3] == 5)
    {
      mGene <- dbTab1[i,1]
      cGene <- dbTab1[i,2]
      mGeneClass <- resultList[[mGene]][[2]]
      colnames(mGeneClass) <- c("cell_line", "mRNA_Expression", "mode")
      mGeneClass <- merge(mGeneClass, tissueMatrix)
      unD <- unique(mGeneClass[,4])
      disSubTab <- as.data.frame(do.call("rbind",lapply(1:length(unD), function(j)
      {
        mGeneClass2 <- mGeneClass[which(mGeneClass$tissue %in% unD[j]),]
        cl1 <- mGeneClass2[which(mGeneClass2[,3] == 1),1]
        cl2 <- mGeneClass2[which(mGeneClass2[,3] == 2),1]
        cl3 <- mGeneClass2[which(mGeneClass2[,3] == 3),1]
        cl4 <- mGeneClass2[which(mGeneClass2[,3] == 4),1]
        cl5 <- mGeneClass2[which(mGeneClass2[,3] == 5),1]
        if(length(cl1) > 2 & length(cl2) > 2)
        {
          if(length(which(!is.na(as.numeric(crisprMatrix[cGene,cl1])))) > 2 & length(which(!is.na(as.numeric(crisprMatrix[cGene,cl2])))) > 2)
          {
            fc1 <- mean(as.numeric(crisprMatrix[cGene,cl1]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
            pv1 <- tryCatch(t.test(as.numeric(crisprMatrix[cGene,cl1]), as.numeric(crisprMatrix[cGene,cl2]))$p.value, error = function(e) NA_real_)
            m1 <- mean(as.numeric(crisprMatrix[cGene,cl1]), na.rm=T)
            m2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
          } else {
            fc1 <- mean(as.numeric(crisprMatrix[cGene,cl1]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
            pv1 <- NA
            m1 <- mean(as.numeric(crisprMatrix[cGene,cl1]), na.rm=T)
            m2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
          }
        }
        else if(length(cl1) == 1 | length(cl2) == 1)
        {
          fc1 <- mean(as.numeric(crisprMatrix[cGene,cl1]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
          pv1 <- NA
          m1 <- mean(as.numeric(crisprMatrix[cGene,cl1]), na.rm=T)
          m2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
        }
        else
        {
          fc1 <- NA
          pv1 <- NA
          m1 <- NA
          m2 <- NA
        }
        # Mode 3 vs. 2
        if(length(cl2) > 2 & length(cl3) > 2)
        {
          if(length(which(!is.na(as.numeric(crisprMatrix[cGene,cl2])))) > 2 & length(which(!is.na(as.numeric(crisprMatrix[cGene,cl3])))) > 2)
          {
            fc2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T)
            pv2 <- tryCatch(t.test(as.numeric(crisprMatrix[cGene,cl2]), as.numeric(crisprMatrix[cGene,cl3]))$p.value, error = function(e) NA_real_)
            m2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
            m3 <- mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T)
          } else {
            fc2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T)
            pv2 <- NA
            m2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
            m3 <- mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T)
          }
        }
        else if(length(cl2) == 1 | length(cl3) == 1)
        {
          fc2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T)
          pv2 <- NA
          m2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
          m3 <- mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T)
        }
        else
        {
          fc2 <- NA
          pv2 <- NA
          m2 <- NA
          m3 <- NA
        }
        # Mode 4 vs. 3
        if(length(cl3) > 2 & length(cl4) > 2)
        {
          if(length(which(!is.na(as.numeric(crisprMatrix[cGene,cl3])))) > 2 & length(which(!is.na(as.numeric(crisprMatrix[cGene,cl4])))) > 2)
          {
            fc3 <- mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl4]), na.rm=T)
            pv3 <- tryCatch(t.test(as.numeric(crisprMatrix[cGene,cl3]), as.numeric(crisprMatrix[cGene,cl4]))$p.value, error = function(e) NA_real_)
            m3 <- mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T)
            m4 <- mean(as.numeric(crisprMatrix[cGene,cl4]), na.rm=T)
          } else {
            fc3 <- mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl4]), na.rm=T)
            pv3 <- NA
            m3 <- mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T)
            m4 <- mean(as.numeric(crisprMatrix[cGene,cl4]), na.rm=T)
          }
        }
        else if(length(cl3) == 1 | length(cl4) == 1)
        {
          fc3 <- mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl4]), na.rm=T)
          pv3 <- NA
          m3 <- mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T)
          m4 <- mean(as.numeric(crisprMatrix[cGene,cl4]), na.rm=T)
        }
        else
        {
          fc3 <- NA
          pv3 <- NA
          m3 <- NA
          m4 <- NA
        }
        # Mode 5 vs. 4
        if(length(cl4) > 2 & length(cl5) > 2)
        {
          if(length(which(!is.na(as.numeric(crisprMatrix[cGene,cl4])))) > 2 & length(which(!is.na(as.numeric(crisprMatrix[cGene,cl5])))) > 2)
          {
            fc4 <- mean(as.numeric(crisprMatrix[cGene,cl4]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl5]), na.rm=T)
            pv4 <- tryCatch(t.test(as.numeric(crisprMatrix[cGene,cl4]), as.numeric(crisprMatrix[cGene,cl5]))$p.value, error = function(e) NA_real_)
            m4 <- mean(as.numeric(crisprMatrix[cGene,cl4]), na.rm=T)
            m5 <- mean(as.numeric(crisprMatrix[cGene,cl5]), na.rm=T)
          } else {
            fc4 <- mean(as.numeric(crisprMatrix[cGene,cl4]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl5]), na.rm=T)
            pv4 <- NA
            m4 <- mean(as.numeric(crisprMatrix[cGene,cl4]), na.rm=T)
            m5 <- mean(as.numeric(crisprMatrix[cGene,cl5]), na.rm=T)
          }
        }
        else if(length(cl4) == 1 | length(cl5) == 1)
        {
          fc4 <- mean(as.numeric(crisprMatrix[cGene,cl4]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl5]), na.rm=T)
          pv4 <- NA
          m4 <- mean(as.numeric(crisprMatrix[cGene,cl4]), na.rm=T)
          m5 <- mean(as.numeric(crisprMatrix[cGene,cl5]), na.rm=T)
        }
        else
        {
          fc4 <- NA
          pv4 <- NA
          m4 <- NA
          m5 <- NA
        }
        return(c(mGene, cGene, unD[j], fc1, fc2, fc3,fc4, pv1, pv2, pv3, pv4, length(cl1), length(cl2), length(cl3), length(cl4), length(cl5),
                 m1, m2, m3, m4, m5, 5))
      })), stringsAsFactors=FALSE)
      for(k in 4:dim(disSubTab)[2]){disSubTab[,k] <- as.numeric(disSubTab[,k])}
      colnames(disSubTab) <- c("gene1", "DEMETER_Gene", "tissue", "fc1", "fc2", "fc3", "fc4", "pv1", "pv2", "pv3", "pv4", "len1", "len2", "len3", "len4", "len5", "mean1", "mean2", "mean3", "mean4", "mean5", "clusters")
      return(disSubTab)
      # listTest[[i]] <- disSubTab  ## IO commented out and uncommented the above line
    }
    else if(dbTab1[i,3] == 4)
    {
      mGene <- dbTab1[i,1]
      cGene <- dbTab1[i,2]
      mGeneClass <- resultList[[mGene]][[2]]
      colnames(mGeneClass) <- c("cell_line", "mRNA_Expression", "mode")
      mGeneClass <- merge(mGeneClass, tissueMatrix)
      unD <- unique(mGeneClass[,4])
      disSubTab <- as.data.frame(do.call("rbind",lapply(1:length(unD), function(j)
      {
        mGeneClass2 <- mGeneClass[which(mGeneClass$tissue %in% unD[j]),]
        cl1 <- mGeneClass2[which(mGeneClass2[,3] == 1),1]
        cl2 <- mGeneClass2[which(mGeneClass2[,3] == 2),1]
        cl3 <- mGeneClass2[which(mGeneClass2[,3] == 3),1]
        cl4 <- mGeneClass2[which(mGeneClass2[,3] == 4),1]
        if(length(cl1) > 2 & length(cl2) > 2)
        {
          if(length(which(!is.na(as.numeric(crisprMatrix[cGene,cl1])))) > 2 & length(which(!is.na(as.numeric(crisprMatrix[cGene,cl2])))) > 2)
          {
            fc1 <- mean(as.numeric(crisprMatrix[cGene,cl1]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
            pv1 <- tryCatch(t.test(as.numeric(crisprMatrix[cGene,cl1]), as.numeric(crisprMatrix[cGene,cl2]))$p.value, error = function(e) NA_real_)
            m1 <- mean(as.numeric(crisprMatrix[cGene,cl1]), na.rm=T)
            m2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
          } else {
            fc1 <- mean(as.numeric(crisprMatrix[cGene,cl1]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
            pv1 <- NA
            m1 <- mean(as.numeric(crisprMatrix[cGene,cl1]), na.rm=T)
            m2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
          }
        }
        else if(length(cl1) == 1 | length(cl2) == 1)
        {
          fc1 <- mean(as.numeric(crisprMatrix[cGene,cl1]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
          pv1 <- NA
          m1 <- mean(as.numeric(crisprMatrix[cGene,cl1]), na.rm=T)
          m2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
        }
        else
        {
          fc1 <- NA
          pv1 <- NA
          m1 <- NA
          m2 <- NA
        }
        # Mode 3 vs. 2
        if(length(cl2) > 2 & length(cl3) > 2)
        {
          if(length(which(!is.na(as.numeric(crisprMatrix[cGene,cl2])))) > 2 & length(which(!is.na(as.numeric(crisprMatrix[cGene,cl3])))) > 2)
          {
            fc2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T)
            pv2 <- tryCatch(t.test(as.numeric(crisprMatrix[cGene,cl2]), as.numeric(crisprMatrix[cGene,cl3]))$p.value, error = function(e) NA_real_)
            m2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
            m3 <- mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T)
          } else {
            fc2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T)
            pv2 <- NA
            m2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
            m3 <- mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T)
          }
        }
        else if(length(cl2) == 1 | length(cl3) == 1)
        {
          fc2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T)
          pv2 <- NA
          m2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
          m3 <- mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T)
        }
        else
        {
          fc2 <- NA
          pv2 <- NA
          m2 <- NA
          m3 <- NA
        }
        # Mode 4 vs. 3
        if(length(cl3) > 2 & length(cl4) > 2)
        {
          if(length(which(!is.na(as.numeric(crisprMatrix[cGene,cl3])))) > 2 & length(which(!is.na(as.numeric(crisprMatrix[cGene,cl4])))) > 2)
          {
            fc3 <- mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl4]), na.rm=T)
            pv3 <- tryCatch(t.test(as.numeric(crisprMatrix[cGene,cl3]), as.numeric(crisprMatrix[cGene,cl4]))$p.value, error = function(e) NA_real_)
            m3 <- mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T)
            m4 <- mean(as.numeric(crisprMatrix[cGene,cl4]), na.rm=T)
          } else {
            fc3 <- mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl4]), na.rm=T)
            pv3 <- NA
            m3 <- mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T)
            m4 <- mean(as.numeric(crisprMatrix[cGene,cl4]), na.rm=T)
          }
        }
        else if(length(cl3) == 1 | length(cl4) == 1)
        {
          fc3 <- mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl4]), na.rm=T)
          pv3 <- NA
          m3 <- mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T)
          m4 <- mean(as.numeric(crisprMatrix[cGene,cl4]), na.rm=T)
        }
        else
        {
          fc3 <- NA
          pv3 <- NA
          m3 <- NA
          m4 <- NA
        }
        return(c(mGene, cGene, unD[j], fc1, fc2, fc3, NA, pv1, pv2, pv3, NA, length(cl1), length(cl2), length(cl3), length(cl4), NA,
                 m1, m2, m3, m4, NA, 4))
      })), stringsAsFactors=FALSE)
      for(k in 4:dim(disSubTab)[2]){disSubTab[,k] <- as.numeric(disSubTab[,k])}
      colnames(disSubTab) <- c("gene1", "DEMETER_Gene", "tissue", "fc1", "fc2", "fc3", "fc4", "pv1", "pv2", "pv3", "pv4", "len1", "len2", "len3", "len4", "len5", "mean1", "mean2", "mean3", "mean4", "mean5", "clusters")
      return(disSubTab)
      # listTest[[i]] <- disSubTab  ## IO commented out and uncommented the above line
    }
    else if(dbTab1[i,3] == 3)
    {
      mGene <- dbTab1[i,1]
      cGene <- dbTab1[i,2]
      mGeneClass <- resultList[[mGene]][[2]]
      colnames(mGeneClass) <- c("cell_line", "mRNA_Expression", "mode")
      mGeneClass <- merge(mGeneClass, tissueMatrix)
      unD <- unique(mGeneClass[,4])
      disSubTab <- as.data.frame(do.call("rbind",lapply(1:length(unD), function(j)
      {
        mGeneClass2 <- mGeneClass[which(mGeneClass$tissue %in% unD[j]),]
        cl1 <- mGeneClass2[which(mGeneClass2[,3] == 1),1]
        cl2 <- mGeneClass2[which(mGeneClass2[,3] == 2),1]
        cl3 <- mGeneClass2[which(mGeneClass2[,3] == 3),1]
        if(length(cl1) > 2 & length(cl2) > 2)
        {
          if(length(which(!is.na(as.numeric(crisprMatrix[cGene,cl1])))) > 2 & length(which(!is.na(as.numeric(crisprMatrix[cGene,cl2])))) > 2)
          {
            fc1 <- mean(as.numeric(crisprMatrix[cGene,cl1]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
            pv1 <- tryCatch(t.test(as.numeric(crisprMatrix[cGene,cl1]), as.numeric(crisprMatrix[cGene,cl2]))$p.value, error = function(e) NA_real_)
            m1 <- mean(as.numeric(crisprMatrix[cGene,cl1]), na.rm=T)
            m2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
          } else {
            fc1 <- mean(as.numeric(crisprMatrix[cGene,cl1]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
            pv1 <- NA
            m1 <- mean(as.numeric(crisprMatrix[cGene,cl1]), na.rm=T)
            m2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
          }
        }
        else if(length(cl1) == 1 | length(cl2) == 1)
        {
          fc1 <- mean(as.numeric(crisprMatrix[cGene,cl1]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
          pv1 <- NA
          m1 <- mean(as.numeric(crisprMatrix[cGene,cl1]), na.rm=T)
          m2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
        }
        else
        {
          fc1 <- NA
          pv1 <- NA
          m1 <- NA
          m2 <- NA
        }
        # Mode 3 vs. 2
        if(length(cl2) > 2 & length(cl3) > 2)
        {
          if(length(which(!is.na(as.numeric(crisprMatrix[cGene,cl2])))) > 2 & length(which(!is.na(as.numeric(crisprMatrix[cGene,cl3])))) > 2)
          {
            fc2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T)
            pv2 <- tryCatch(t.test(as.numeric(crisprMatrix[cGene,cl2]), as.numeric(crisprMatrix[cGene,cl3]))$p.value, error = function(e) NA_real_)
            m2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
            m3 <- mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T)
          } else {
            fc2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T)
            pv2 <- NA
            m2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
            m3 <- mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T)
          }
        }
        else if(length(cl2) == 1 | length(cl3) == 1)
        {
          fc2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T)
          pv2 <- NA
          m2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
          m3 <- mean(as.numeric(crisprMatrix[cGene,cl3]), na.rm=T)
        }
        else
        {
          fc2 <- NA
          pv2 <- NA
          m2 <- NA
          m3 <- NA
        }
        return(c(mGene, cGene, unD[j], fc1, fc2, 0, 0, pv1, pv2, 0, 0, length(cl1), length(cl2), length(cl3), 0, 0, m1, m2, m3, 0, 0, 3))
      })), stringsAsFactors=FALSE)
      for(k in 4:dim(disSubTab)[2]){disSubTab[,k] <- as.numeric(disSubTab[,k])}
      colnames(disSubTab) <- c("gene1", "DEMETER_Gene", "tissue", "fc1", "fc2", "fc3", "fc4", "pv1", "pv2", "pv3", "pv4", "len1", "len2", "len3", "len4", "len5", "mean1", "mean2", "mean3", "mean4", "mean5", "clusters")
      return(disSubTab)
      # listTest[[i]] <- disSubTab  ## IO commented out and uncommented the above line

    }
    else if(dbTab1[i,3] == 2)
    {
      mGene <- dbTab1[i,1]
      cGene <- dbTab1[i,2]
      mGeneClass <- resultList[[mGene]][[2]]
      colnames(mGeneClass) <- c("cell_line", "mRNA_Expression", "mode")
      mGeneClass <- merge(mGeneClass, tissueMatrix)
      unD <- unique(mGeneClass[,4])
      disSubTab <- as.data.frame(do.call("rbind",lapply(1:length(unD), function(j)
      {
        mGeneClass2 <- mGeneClass[which(mGeneClass$tissue %in% unD[j]),]
        cl1 <- mGeneClass2[which(mGeneClass2[,3] == 1),1]
        cl2 <- mGeneClass2[which(mGeneClass2[,3] == 2),1]
        if(length(cl1) > 2 & length(cl2) > 2)
        {
          if(length(which(!is.na(as.numeric(crisprMatrix[cGene,cl1])))) > 2 & length(which(!is.na(as.numeric(crisprMatrix[cGene,cl2])))) > 2)
          {
            fc1 <- mean(as.numeric(crisprMatrix[cGene,cl1]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
            pv1 <- tryCatch(t.test(as.numeric(crisprMatrix[cGene,cl1]), as.numeric(crisprMatrix[cGene,cl2]))$p.value, error = function(e) NA_real_)
            m1 <- mean(as.numeric(crisprMatrix[cGene,cl1]), na.rm=T)
            m2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
          } else {
            fc1 <- mean(as.numeric(crisprMatrix[cGene,cl1]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
            pv1 <- NA
            m1 <- mean(as.numeric(crisprMatrix[cGene,cl1]), na.rm=T)
            m2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
          }
        }
        else if(length(cl1) == 1 | length(cl2) == 1)
        {
          fc1 <- mean(as.numeric(crisprMatrix[cGene,cl1]), na.rm=T) - mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
          pv1 <- NA
          m1 <- mean(as.numeric(crisprMatrix[cGene,cl1]), na.rm=T)
          m2 <- mean(as.numeric(crisprMatrix[cGene,cl2]), na.rm=T)
        }
        else
        {
          fc1 <- NA
          pv1 <- NA
          m1 <- NA
          m2 <- NA
        }
        return(c(mGene, cGene, unD[j], fc1, 0, 0, 0, pv1, 0, 0, 0, length(cl1), length(cl2), 0, 0, 0, m1, m2, 0, 0, 0, 2))
      })), stringsAsFactors=FALSE)
      for(k in 4:dim(disSubTab)[2]){disSubTab[,k] <- as.numeric(disSubTab[,k])}
      colnames(disSubTab) <- c("gene1", "DEMETER_Gene", "tissue", "fc1", "fc2", "fc3", "fc4", "pv1", "pv2", "pv3", "pv4", "len1", "len2", "len3", "len4", "len5", "mean1", "mean2", "mean3", "mean4", "mean5", "clusters")
      return(disSubTab)
      # listTest[[i]] <- disSubTab   ## IO commented out and uncommented the above line

    }
    else
    {
      # Do Nothing
    }
  }, mc.cores=cores)
  
  # New processing results section (for speed)
  # n=2 results
  
  if(dim(listGroupD[[1]])[2] > 1)
  {
    for(i in 1:length(listGroupD))
    {
      if(unique(listGroupD[[i]]$clusters == 2))
      {
        listGroupD[[i]]$qv1 <- p.adjust(listGroupD[[i]]$pv1, method="BH")
        listGroupD[[i]]$qv2 <- NA
        listGroupD[[i]]$qv3 <- NA
        listGroupD[[i]]$qv4 <- NA
      }
      else if(unique(listGroupD[[i]]$clusters == 3))
      {
        listGroupD[[i]]$qv1 <- p.adjust(listGroupD[[i]]$pv1, method="BH")
        listGroupD[[i]]$qv2 <- p.adjust(listGroupD[[i]]$pv2, method="BH")
        listGroupD[[i]]$qv3 <- NA
        listGroupD[[i]]$qv4 <- NA
      }
      else if(unique(listGroupD[[i]]$clusters == 4))
      {
        listGroupD[[i]]$qv1 <- p.adjust(listGroupD[[i]]$pv1, method="BH")
        listGroupD[[i]]$qv2 <- p.adjust(listGroupD[[i]]$pv2, method="BH")
        listGroupD[[i]]$qv3 <- p.adjust(listGroupD[[i]]$pv3, method="BH")
        listGroupD[[i]]$qv4 <- NA
      }
      else if(unique(listGroupD[[i]]$clusters == 5))
      {
        listGroupD[[i]]$qv1 <- p.adjust(listGroupD[[i]]$pv1, method="BH")
        listGroupD[[i]]$qv2 <- p.adjust(listGroupD[[i]]$pv2, method="BH")
        listGroupD[[i]]$qv3 <- p.adjust(listGroupD[[i]]$pv3, method="BH")
        listGroupD[[i]]$qv4 <- p.adjust(listGroupD[[i]]$pv4, method="BH")
      }
    }
    listGroupTABd <- do.call("rbind", listGroupD)
    
    tabD <- listGroupTABd[,c(1,2,3,22,12,13,14,15,16,17,18,19,20,21,4,5,6,7,8,9,10,11,23,24,25,26)]
    colnames(tabD) <- c("mRNA_gene", "crispr_gene", "tissue", "num_modes", "nMode1", "nMode2", "nMode3", "nMode4",
                        "nMode5", "mMode1", "mMode2", "mMode3", "mMode4","mMode5", "shift2_1", "shift3_2",
                        "shift4_3", "shift5_4", "pvalue2_1", "pvalue3_2", "pvalue4_3", "pvalue5_4",
                        "qvalue2_1", "qvalue3_2", "qvalue4_3", "qvalue5_4")
    for(i in 4:dim(tabD)[2]) { tabD[,i] <- as.numeric(tabD[,i])}
  } else {
    tabD <- matrix(ncol=26, nrow=1)
    colnames(tabD) <- c("mRNA_gene", "crispr_gene", "tissue", "num_modes", "nMode1", "nMode2", "nMode3", "nMode4",
                        "nMode5", "mMode1", "mMode2", "mMode3", "mMode4","mMode5", "shift2_1", "shift3_2",
                        "shift4_3", "shift5_4", "pvalue2_1", "pvalue3_2", "pvalue4_3", "pvalue5_4",
                        "qvalue2_1", "qvalue3_2", "qvalue4_3", "qvalue5_4")
  }
  
  if(length(which(is.na(tabD[,1]))) > 0)
  {
    for(i in 4:dim(tabD)[2]) {
      tabD[,i] <- round(tabD[,i], 2)
    }
    tabD <- tabD[-which(is.na(tabD[,1])),]
    return(tabD)
  } else {
    return(tabD)
  }
}