R/DIscBIO-generic-ClustDiffGenes.R

In DIscBIO: A User-Friendly Pipeline for Biomarker Discovery in Single-Cell Transcriptomics

#' @title ClustDiffGenes
#' @description Creates a table of cluster differences
#' @param object \code{DISCBIO} class object.
#' @param K A numeric value of the number of clusters.
#' @param pValue A numeric value of the p-value. Default is 0.05.
#' @param fdr A numeric value of the false discovery rate. Default is 0.01.
#' @param export A logical vector that allows writing the final gene list in
#'   excel file. Default is TRUE.
#' @param quiet if `TRUE`, suppresses intermediate text output
#' @param filename_up Name of the exported "up" file (if `export=TRUE`)
#' @param filename_down Name of the exported "down" file (if `export=TRUE`)
#' @param filename_binom Name of the exported binomial file
#' @param filename_sigdeg Name of the exported sigDEG file
#' @importFrom stats pbinom median
#' @import org.Hs.eg.db
#' @rdname ClustDiffGenes
#' @return A list containing two tables.
#' @export
#' @examples
#' sc <- DISCBIO(valuesG1msTest)
#' sc <- Clustexp(sc, cln = 3, quiet = TRUE)
#' cdiff <- ClustDiffGenes(sc, K = 3, fdr = .3, export = FALSE)
#' str(cdiff)
#' cdiff[[2]]
setGeneric(
  "ClustDiffGenes",
  function(
      object, K, pValue = 0.05, fdr = .01, export = FALSE, quiet = FALSE,
      filename_up = "Up-DEG-cluster",
      filename_down = "Down-DEG-cluster",
      filename_binom = "binomial-DEGsTable",
      filename_sigdeg = "binomial-sigDEG") {
    standardGeneric("ClustDiffGenes")
  }
)
#' @export
#' @rdname ClustDiffGenes
setMethod(
  "ClustDiffGenes",
  signature = "DISCBIO",
  definition = function(
      object, K, pValue, fdr, export, quiet, filename_up, filename_down,
      filename_binom, filename_sigdeg) {
    # ======================================================================
    # Validating
    # ======================================================================
    ran_k <- length(object@kmeans$kpart) > 0
    ran_m <- length(object@MBclusters) > 0
    if (!is.numeric(fdr)) {
      stop("fdr has to be a number between 0 and 1")
    } else if (fdr < 0 | fdr > 1) {
      stop("fdr has to be a number between 0 and 1")
    }
    if (!is.numeric(pValue)) {
      stop("pValue has to be a number between 0 and 1")
    } else if (pValue < 0 | pValue > 1) {
      stop("pValue has to be a number between 0 and 1")
    }
    cdiff <- list()
    x <- object@ndata
    y <- object@expdata[, names(object@ndata)]
    if (ran_k) {
      part <- object@kmeans$kpart
    } else if (ran_m) {
      part <- object@MBclusters$clusterid
    } else {
      stop("Run clustering before running this function")
    }
    # ======================================================================
    # Operating
    # ======================================================================
    for (i in 1:max(part)) {
      if (sum(part == i) == 0) {
        next
      }
      m <- apply(x, 1, mean)
      n <-
        if (sum(part == i) > 1) {
          apply(x[, part == i], 1, mean)
        } else {
          x[, part == i]
        }
      no <-
        if (sum(part == i) > 1) {
          median(apply(y[, part == i], 2, sum)) /
            median(apply(x[, part == i], 2, sum))
        } else {
          sum(y[, part == i]) / sum(x[, part == i])
        }
      m <- m * no
      n <- n * no
      pv <- binompval(m / sum(m), sum(n), n)
      d <-
        data.frame(
          mean.all = m,
          mean.cl = n,
          fc = n / m,
          pv = pv
        )[order(pv, decreasing = FALSE), ]
      cdiff[[i]] <- d[d$pv < pValue, ]
    }
    DEGsE <- vector()
    DEGsS <- vector()
    DEGsTable <- data.frame()

    for (n in 1:K) {
      if (length(cdiff[[n]][, 1]) == 0) {
        next
      }

      if (length(cdiff[[n]][, 1]) > 0) {
        p.adj <- p.adjust(cdiff[[n]][, 4], method = "bonferroni")
        out <- cbind(cdiff[[n]], p.adj)
        out <- subset(out, out[, 5] < fdr)
        if (length(out[, 1]) > 0) {
          Regulation <- vector()
          for (i in seq_len(length(out[, 1]))) {
            if (out[i, 1] > out[i, 2]) {
              Regulation[i] <- "Down"
            } else {
              Regulation[i] <- "Up"
            }
          }
          out <- cbind(out, Regulation)
          if (quiet) {
            suppressMessages(
              geneList <-
                AnnotationDbi::select(
                  org.Hs.eg.db,
                  keys = keys(org.Hs.eg.db),
                  columns = c("SYMBOL", "ENSEMBL")
                )
            )
            GL <- c(1, "MTRNR2", "ENSG00000210082")
            GL1 <- c(1, "MTRNR1", "ENSG00000211459")
            geneList <- rbind(geneList, GL, GL1)
          } else {
            geneList <-
              AnnotationDbi::select(
                org.Hs.eg.db,
                keys = keys(org.Hs.eg.db),
                columns = c("SYMBOL", "ENSEMBL")
              )
            GL <- c(1, "MTRNR2", "ENSG00000210082")
            GL1 <- c(1, "MTRNR1", "ENSG00000211459")
            geneList <- rbind(geneList, GL, GL1)
          }
          genes <- rownames(out)
          gene_list <- geneList[, 3]
          idx_genes <- is.element(gene_list, genes)
          genes2 <- geneList[idx_genes, ]
          Final <- cbind(genes, out)

          Final <-
            merge(
              Final,
              genes2,
              by.x = "genes",
              by.y = "ENSEMBL",
              all.x = TRUE
            )
          Final <- Final[!duplicated(Final[, 1]), ]
          Final[is.na(Final[, 9]), c(1, 9)] <-
            Final[is.na(Final[, 9]), 1]
          rownames(Final) <- Final[, 1]
          Final[, 1] <- Final[, 9]
          Final <- Final[, -9]
          DEGsS <- c(DEGsS, Final[, 1])
          DEGsE <-
            c(DEGsE, as.character(rownames(Final)))
          Up <- subset(Final, Final[, 7] == "Up")
          cols_to_keep <- c(
            "Regulation", "genes", "pv", "mean.all", "mean.cl", "fc", "p.adj"
          )
          Up <- Up[, cols_to_keep]
          Up[, 3] <- rownames(Up)
          Up[, 6] <- log2(Up[, 6])
          Up[, 1] <- Up[, 2]
          colnames(Up) <- c(
            "Genes",
            "genes",
            "E.genes",
            "mean.all",
            "mean.cl",
            "log2.fc",
            "p.adj"
          )
          if (export) {
            write.csv(
              Up,
              file = paste0(filename_up, n, ".csv")
            )
          }

          Down <- subset(Final, Final[, 7] == "Down")
          cols_to_keep <- c(
            "Regulation", "genes", "pv", "mean.all", "mean.cl", "fc", "p.adj"
          )
          Down <- Down[, cols_to_keep]
          Down[, 3] <- rownames(Down)
          Down[, 6] <- log2(Down[, 6])
          Down[, 1] <- Down[, 2]
          colnames(Down) <- c(
            "Genes",
            "genes",
            "E.genes",
            "mean.all",
            "mean.cl",
            "log2.fc",
            "p.adj"
          )
          if (export) {
            write.csv(
              Down,
              file = paste0(filename_down, n, ".csv")
            )
          }

          sigDEG <- cbind(DEGsE, DEGsS)
          if (export) {
            write.csv(
              sigDEG,
              file = paste0(filename_sigdeg, ".csv")
            )
          }

          DEGsTable[n, 1] <- paste0("Cluster ", n)
          DEGsTable[n, 2] <- "Remaining Clusters"
          DEGsTable[n, 3] <- length(Up[, 1])
          DEGsTable[n, 4] <- paste0(filename_up, n, ".csv")
          DEGsTable[n, 5] <- length(Down[, 1])
          DEGsTable[n, 6] <- paste0(filename_down, n, ".csv")
        }
      }
    }
    if (length(DEGsTable) > 0) {
      colnames(DEGsTable) <- c(
        "Target Cluster", "VS", "Gene number", "File name",
        "Gene number", "File name"
      )
      if (export) {
        write.csv(DEGsTable, file = paste0(filename_binom, ".csv"))
      }
      return(list(sigDEG, DEGsTable))
    } else {
      print(paste("There are no DEGs with fdr =", fdr))
    }
  }
)