R/DIscBIO-generic-DEGanalysis2clust.R

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

#' @title Determining differentially expressed genes (DEGs) between two
#'   particular clusters.
#' @description This function defines DEGs between particular clusters generated
#'   by either K-means or model based clustering.
#' @param object \code{DISCBIO} class object.
#' @param Clustering Clustering has to be one of the following:
#'   ["K-means","MB"]. Default is "K-means"
#' @param K A numeric value of the number of clusters.
#' @param fdr A numeric value of the false discovery rate. Default is 0.05.
#' @param name A string vector showing the name to be used to save the resulted
#'   tables.
#' @param First A string vector showing the first target cluster.  Default is
#'   "CL1"
#' @param Second A string vector showing the second target cluster.  Default is
#'   "CL2"
#' @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 plot if `TRUE`, plots are generated
#' @param filename_deg Name of the exported DEG table
#' @param filename_sigdeg Name of the exported sigDEG table
#' @param ... additional parameters to be passed to samr()
#' @importFrom graphics title
#' @importFrom utils write.csv capture.output
#' @importFrom AnnotationDbi keys
#' @return A list containing two tables.
setGeneric(
  "DEGanalysis2clust",
  function(object, K, Clustering = "K-means", fdr = 0.05, name = "Name",
           First = "CL1", Second = "CL2", export = FALSE, quiet = FALSE,
           plot = TRUE, filename_deg = "DEGsTable", filename_sigdeg = "sigDEG",
           ...) {
    standardGeneric("DEGanalysis2clust")
  }
)

#' @export
#' @rdname DEGanalysis2clust
setMethod(
  "DEGanalysis2clust",
  signature = "DISCBIO",
  definition = function(
    object, K, Clustering, fdr, name, First, Second, export, quiet, plot,
    filename_deg, filename_sigdeg, ...
  ) {
    if (!(Clustering %in% c("K-means", "MB"))) {
      stop("Clustering has to be either K-means or MB")
    }
    gene_list <- object@FinalGeneList
    gene_names <- rownames(object@expdata)
    idx_genes <- is.element(gene_names, gene_list)
    gene_names2 <- gene_names[idx_genes]
    dataset <- object@expdata[gene_names2, ]
    Nam <- colnames(dataset)
    if (Clustering == "K-means") {
      Cluster_ID <- object@cpart
      if (length(object@cpart) < 1) {
        stop("run Clustexp before running DEGanalysis2clust")
      }
    }

    if (Clustering == "MB") {
      Cluster_ID <- object@MBclusters$clusterid
      if (length(object@MBclusters$clusterid) < 1) {
        stop("run ExprmclustMB before running DEGanalysis2clust")
      }
    }
    num <- seq_len(K)
    num1 <- paste("CL", num, sep = "")
    for (n in num) {
      Nam <- ifelse((Cluster_ID == n), num1[n], Nam)
    }
    colnames(dataset) <- Nam
    sg1 <- dataset[, which(colnames(dataset) == First)]
    sg2 <- dataset[, which(colnames(dataset) == Second)]
    sg <- cbind(sg1, sg2)

    sg3 <- factor(
      gsub(paste0("(", First, "|", Second, ").*"), "\\1", colnames(sg)),
      levels = c(paste0(First), paste0(Second))
    )
    sg3 <- sg3[!is.na(sg3)]

    colnames(sg) <- sg3
    len <- c(
      length(sg[, which(colnames(sg) == First)]),
      length(sg[, which(colnames(sg) == Second)])
    )
    y <- c(rep(1:2, len))
    L <- as.matrix(sg)
    gname <- rownames(sg)
    x <- L
    data <- list(x = x, y = y, geneid = gname)
    if (quiet) {
      suppressMessages(
        invisible(capture.output({
          samr.obj <- sammy(
            data,
            resp.type = "Two class unpaired",
            assay.type = "seq",
            testStatistic = "wilcoxon",
            random.seed = 15,
            ...
          )
          delta.table <- samr.compute.delta.table(samr.obj)
        }))
      )
    } else {
      samr.obj <- sammy(
        data,
        resp.type = "Two class unpaired",
        assay.type = "seq",
        testStatistic = "wilcoxon",
        random.seed = 15,
        ...
      )
      delta.table <- samr.compute.delta.table(samr.obj)
    }
    DEGsTable <- data.frame()
    DEGsE <- vector()
    DEGsS <- vector()
    wm <- which.min(delta.table[, 5])
    if (delta.table[wm, 5] <= fdr) {
      w <- which(delta.table[, 5] <= fdr)
      if (is.null(w)) stop("No suitable deltas. Try a lower FDR value.")
      delta <- delta.table[w[1], 1] - 0.001
      if (plot) {
        samr.plot(samr.obj, delta)
        title(paste("DEGs in the", Second, "in", First, "VS", Second))
      }
      siggenes.table <- samr.compute.siggenes.table(
        samr.obj, delta, data, delta.table
      )
      # ------------------------------------------------------------------
      # Reformat siggenes.table as data.frame
      # ------------------------------------------------------------------
      siggenes.table$genes.lo <- reformatSiggenes(siggenes.table$genes.lo)
      siggenes.table$genes.up <- reformatSiggenes(siggenes.table$genes.up)

      FDRl <- as.numeric(siggenes.table$genes.lo[, 8]) / 100
      FDRu <- as.numeric(siggenes.table$genes.up[, 8]) / 100

      siggenes.table$genes.lo[, 8] <- FDRl
      siggenes.table$genes.up[, 8] <- FDRu

      DEGsTable[1, 1] <- paste0(First, " VS ", Second)
      DEGsTable[1, 2] <- Second
      DEGsTable[1, 3] <- length(FDRu)
      DEGsTable[1, 4] <- paste0(
        "Up-regulated-", name, Second, "in", First, "VS", Second,
        ".csv"
      )
      DEGsTable[1, 5] <- length(FDRl)
      DEGsTable[1, 6] <- paste0(
        "Low-regulated-", name, Second, "in", First, "VS", Second,
        ".csv"
      )
      DEGsTable[2, 1] <- paste0(First, " VS ", Second)
      DEGsTable[2, 2] <- First
      DEGsTable[2, 3] <- length(FDRu)
      DEGsTable[2, 4] <- paste0(
        "Low-regulated-", name, First, "in", First, "VS", Second,
        ".csv"
      )
      DEGsTable[2, 5] <- length(FDRl)
      DEGsTable[2, 6] <- paste0(
        "Up-regulated-", name, First, "in", First, "VS", Second,
        ".csv"
      )
      FinalDEGsL <- data.frame()
      if (length(FDRl) > 0) {
        genes <- siggenes.table$genes.lo[, 3]
        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)
        }
        FinalDEGsL <- cbind(genes, siggenes.table$genes.lo)
        gene_list <- geneList[, 3]
        idx_genes <- is.element(gene_list, genes)
        genes2 <- geneList[idx_genes, ]
        if (!is.null(FinalDEGsL)) {
          FinalDEGsL <- merge(
            FinalDEGsL,
            genes2,
            by.x = "genes",
            by.y = "ENSEMBL",
            all.x = TRUE
          )
          FinalDEGsL[, 3] <- FinalDEGsL[, 11]
          FinalDEGsL <- FinalDEGsL[, c(-1, -10, -11)]
          FinalDEGsL <- FinalDEGsL[order(FinalDEGsL[, 8]), ]
          FinalDEGsL[is.na(FinalDEGsL[, 2]), c(2, 3)] <-
            FinalDEGsL[is.na(FinalDEGsL[, 2]), 3]
        }
        if (export) {
          message("The results of DEGs are saved in your directory")
          message(
            "Low-regulated genes in the ", Second, " in ",
            First, " VS ", Second, "\n"
          )
          write.csv(
            FinalDEGsL,
            file = paste0(
              "Low-regulated-", name, Second, "in", First,
              "VS", Second, ".csv"
            )
          )
          write.csv(
            FinalDEGsL,
            file = paste0(
              "Up-regulated-", name, First, "in", First, "VS",
              Second, ".csv"
            )
          )
        }
        DEGsS <- c(DEGsS, FinalDEGsL[, 2])
        DEGsE <- c(DEGsE, as.character(FinalDEGsL[, 3]))
      }
      FinalDEGsU <- data.frame()
      if (length(FDRu) > 0) {
        genes <- siggenes.table$genes.up[, 3]
        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")
          geneList <- rbind(geneList, GL)
        } else {
          geneList <- AnnotationDbi::select(
            org.Hs.eg.db,
            keys = keys(org.Hs.eg.db),
            columns = c("SYMBOL", "ENSEMBL")
          )
          GL <- c(1, "MTRNR2", "ENSG00000210082")
          geneList <- rbind(geneList, GL)
        }
        FinalDEGsU <- cbind(genes, siggenes.table$genes.up)
        gene_list <- geneList[, 3]
        idx_genes <- is.element(gene_list, genes)
        genes2 <- geneList[idx_genes, ]
        if (!is.null(FinalDEGsU)) {
          FinalDEGsU <- merge(
            FinalDEGsU,
            genes2,
            by.x = "genes",
            by.y = "ENSEMBL",
            all.x = TRUE
          )
          FinalDEGsU[, 3] <- FinalDEGsU[, 11]
          FinalDEGsU <- FinalDEGsU[, c(-1, -10, -11)]
          FinalDEGsU <- FinalDEGsU[order(FinalDEGsU[, 8]), ]
          FinalDEGsU[is.na(FinalDEGsU[, 2]), c(2, 3)] <-
            FinalDEGsU[is.na(FinalDEGsU[, 2]), 3]
        }
        if (export) {
          message("The results of DEGs are saved in your directory")
          message(
            "Up-regulated genes in the ", Second, " in ", First,
            " VS ", Second, "\n"
          )
          write.csv(
            FinalDEGsU,
            file = paste0(
              "Up-regulated-", name, Second, "in", First, "VS",
              Second, ".csv"
            )
          )
          write.csv(
            FinalDEGsU,
            file = paste0(
              "Low-regulated-", name, First, "in", First, "VS",
              Second, ".csv"
            )
          )
        }
        DEGsS <- c(DEGsS, FinalDEGsU[, 2])
        DEGsE <- c(DEGsE, as.character(FinalDEGsU[, 3]))
      }
    } else {
      DEGsTable[1, 1] <- paste0(First, " VS ", Second)
      DEGsTable[1, 2] <- Second
      DEGsTable[1, 3] <- NA
      DEGsTable[1, 4] <- NA
      DEGsTable[1, 5] <- NA
      DEGsTable[1, 6] <- NA

      DEGsTable[2, 1] <- paste0(First, " VS ", Second)
      DEGsTable[2, 2] <- First
      DEGsTable[2, 3] <- NA
      DEGsTable[2, 4] <- NA
      DEGsTable[2, 5] <- NA
      DEGsTable[2, 6] <- NA
    }
    colnames(DEGsTable) <- c(
      "Comparisons",
      "Target cluster",
      "Gene number",
      "File name",
      "Gene number",
      "File name"
    )
    if (!quiet) print(DEGsTable)
    sigDEG <- cbind(DEGsE, DEGsS)
    if (export) {
      write.csv(DEGsTable, file = paste0(filename_deg, ".csv"))
      write.csv(sigDEG, file = paste0(filename_sigdeg, ".csv"))
    }
    return(
      list(
        sigDEG = sigDEG,
        DEGsTable = DEGsTable,
        FinalDEGsL = FinalDEGsL,
        FinalDEGsU = FinalDEGsU
      )
    )
  }
)