#' @title Determining differentially expressed genes (DEGs) between all

#'   individual clusters.

#' @description This function defines DEGs between all individual 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 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

#' @importFrom graphics title

#' @importFrom utils write.csv capture.output

#' @param ... additional parameters to be passed to samr()

#' @return A list containing two tables.

setGeneric(

  name = "DEGanalysis",

  def = function(object,

                 K,

                 Clustering = "K-means",

                 fdr = 0.05,

                 name = "Name",

                 export = FALSE,

                 quiet = FALSE,

                 plot = TRUE,

                 filename_deg = "DEGsTable",

                 filename_sigdeg = "sigDEG",

                 ...) {

    standardGeneric("DEGanalysis")

  }

)

#' @export

#' @rdname DEGanalysis

setMethod(

  f = "DEGanalysis",

  signature = "DISCBIO",

  definition = function(object,

                        K,

                        Clustering = "K-means",

                        fdr = 0.05,

                        name = "Name",

                        export = FALSE,

                        quiet = FALSE,

                        plot = TRUE,

                        filename_deg,

                        filename_sigdeg,

                        ...) {

    # Validation

    if (!(Clustering %in% c("K-means", "MB"))) {

      stop("Clustering has to be either K-means or MB")

    }

    if (Clustering == "K-means") {

      Cluster_ID <- object@cpart

      if (length(object@cpart) < 1) {

        stop("run Clustexp before running DEGanalysisM")

      }

    }

    if (Clustering == "MB") {

      Cluster_ID <- object@MBclusters$clusterid

      if (length(object@MBclusters$clusterid) < 1) {

        stop("run ExprmclustMB before running DEGanalysisM")

      }

    }

    # Defining initial objects

    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)

    num <- 1:K

    num1 <- paste("CL", num, sep = "")

    for (n in num) {

      Nam <- ifelse((Cluster_ID == n), num1[n], Nam)

    }

    colnames(dataset) <- Nam

    if (!quiet) {

      message("The dataset is ready for differential expression analysis")

    }

    num1 <- paste("CL", num, sep = "")

    clustName <- paste("Cl", num, sep = "")

    ClusterLength <- vector()

    for (i in num) {

      d1 <- clustName[i]

      d2 <- dataset[, which(names(dataset) == num1[i])]

      ClusterLength[i] <- length(d2)

      assign(d1, d2)

    }

    ccdf <- data.frame(clustName, ClusterLength)

    ccdff <- ccdf[order(ClusterLength), ]

    clustName <- ccdff[, 1]

    if (!quiet) {

      print(clustName)

    }

    first <- vector()

    second <- vector()

    if (K < 2) {

      stop("K has to be at least 2")

    } else if (K == 2) {

      first <- c(paste0(clustName[1]))

      second <- c(paste0(clustName[2]))

    } else if (K == 3) {

      first <- c(

        rep(paste0(clustName[1]), 2), rep(paste0(clustName[2]), 1)

      )

      second <- c(

        paste0(clustName[2]),

        paste0(clustName[3]),

        paste0(clustName[3])

      )

    } else if (K == 4) {

      first <- c(

        rep(paste0(clustName[1]), 3),

        rep(paste0(clustName[2]), 1),

        rep(paste0(clustName[4]), 1),

        rep(paste0(clustName[3]), 1)

      )

      second <- c(

        paste0(clustName[2]),

        paste0(clustName[3]),

        paste0(clustName[4]),

        paste0(clustName[3]),

        paste0(clustName[2]),

        paste0(clustName[4])

      )

    } else if (K == 5) {

      first <- c(

        rep(paste0(clustName[1]), 4),

        rep(paste0(clustName[2]), 3),

        rep(paste0(clustName[3]), 2),

        rep(paste0(clustName[5]), 1)

      )

      second <- c(

        paste0(clustName[2]),

        paste0(clustName[3]),

        paste0(clustName[4]),

        paste0(clustName[5]),

        paste0(clustName[3]),

        paste0(clustName[4]),

        paste0(clustName[5]),

        paste0(clustName[4]),

        paste0(clustName[5]),

        paste0(clustName[4])

      )

    } else if (K == 6) {

      first <- c(

        rep(paste0(clustName[1]), 3),

        rep(paste0(clustName[5]), 1),

        rep(paste0(clustName[1]), 1),

        rep(paste0(clustName[2]), 2),

        rep(paste0(clustName[5]), 1),

        rep(paste0(clustName[2]), 1),

        rep(paste0(clustName[3]), 1),

        rep(paste0(clustName[5]), 1),

        rep(paste0(clustName[3]), 1),

        rep(paste0(clustName[5]), 1),

        rep(paste0(clustName[4]), 1),

        rep(paste0(clustName[5]), 1)

      )

      second <- c(

        paste0(clustName[2]),

        paste0(clustName[3]),

        paste0(clustName[4]),

        paste0(clustName[1]),

        paste0(clustName[6]),

        paste0(clustName[3]),

        paste0(clustName[4]),

        paste0(clustName[2]),

        paste0(clustName[6]),

        paste0(clustName[4]),

        paste0(clustName[3]),

        paste0(clustName[6]),

        paste0(clustName[4]),

        paste0(clustName[6]),

        paste0(clustName[6])

      )

    }

    o <- 1:K

    oo <- o[-length(o)]

    com <- sum(oo)

    if (!quiet) message("Number of comparisons: ", com * 2, "\n")

    comNUM <- paste("comp", seq_len(com), sep = "")

    DEGsTable <- data.frame()

    DEGsE <- vector()

    DEGsS <- vector()

    for (i in seq_len(com)) {

      FinalDEGsL <- data.frame()

      FinalDEGsU <- data.frame()

      FDRl <- vector()

      FDRu <- vector()

      d1 <- comNUM[i]

      d2 <- cbind(get(first[i]), get(second[i]))

      assign(d1, d2)

      len <-

        c(length(get(first[i])[1, ]), length(get(second[i])[1, ]))

      y <- c(rep(1:2, len))

      L <- as.matrix(get(comNUM[i]))

      gname <- rownames(get(comNUM[i]))

      x <- L

      data <- list(x = x, y = y, geneid = gname)

      if (quiet) {

        invisible(capture.output(

          samr.obj <- sammy(

            data,

            resp.type = "Two class unpaired",

            assay.type = "seq",

            testStatistic = "wilcoxon",

            random.seed = 15,

            ...

          )

        ))

      } else {

        samr.obj <- sammy(

          data,

          resp.type = "Two class unpaired",

          assay.type = "seq",

          testStatistic = "wilcoxon",

          random.seed = 15,

          ...

        )

      }

      if (quiet) {

        invisible(capture.output(

          delta.table <- samr.compute.delta.table(samr.obj)

        ))

      } else {

        delta.table <- samr.compute.delta.table(samr.obj)

      }

      wm <- which.min(delta.table[, 5])

      if (delta.table[wm, 5] <= fdr) {

        w <- which(delta.table[, 5] <= fdr)

        delta <- delta.table[w[1], 1] - 0.001

        if (plot) {

          samr.plot(samr.obj, delta)

          title(paste0(

            "DEGs in the ", second[i], " in ", first[i], " VS ",

            second[i]

          ))

        }

        siggenes.table <- samr.compute.siggenes.table(

          samr.obj, delta, data, delta.table

        )

        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[i, 1] <-

          paste0(first[i], " VS ", second[i])

        DEGsTable[i, 2] <- second[i]

        DEGsTable[i, 3] <- length(FDRu)

        DEGsTable[i, 4] <- paste0(

          "Up-regulated-", name, second[i], "in", first[i], "VS",

          second[i], ".csv"

        )

        DEGsTable[i, 5] <- length(FDRl)

        DEGsTable[i, 6] <- paste0(

          "Low-regulated-", name, second[i], "in", first[i], "VS",

          second[i], ".csv"

        )

        DEGsTable[i + com, 1] <-

          paste0(first[i], " VS ", second[i])

        DEGsTable[i + com, 2] <- first[i]

        DEGsTable[i + com, 3] <- length(FDRu)

        DEGsTable[i + com, 4] <- paste0(

          "Low-regulated-", name, first[i], "in", first[i], "VS",

          second[i], ".csv"

        )

        DEGsTable[i + com, 5] <- length(FDRl)

        DEGsTable[i + com, 6] <- paste0(

          "Up-regulated-", name, first[i], "in", first[i], "VS",

          second[i], ".csv"

        )

        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")

            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)

          }

          FinalDEGsL <- cbind(genes, siggenes.table$genes.lo)

          gene_list <- geneList[, 3]

          idx_genes <- is.element(gene_list, genes)

          genes2 <- geneList[idx_genes, ]

          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 (!quiet) {

            message(

              "Low-regulated genes in the ", second[i],

              " in ", first[i], " VS ", second[i], "\n"

            )

          }

          if (export) {

            write.csv(

              FinalDEGsL,

              file = paste0(

                "Low-regulated-", name, second[i], "in",

                first[i], "VS", second[i], ".csv"

              )

            )

            write.csv(

              FinalDEGsL,

              file = paste0(

                "Up-regulated-", name, first[i], "in", first[i],

                "VS", second[i], ".csv"

              )

            )

          }

          DEGsS <- c(DEGsS, FinalDEGsL[, 2])

          DEGsE <- c(DEGsE, as.character(FinalDEGsL[, 3]))

        }

        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")

            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)

          }

          FinalDEGsU <- cbind(genes, siggenes.table$genes.up)

          gene_list <- geneList[, 3]

          idx_genes <- is.element(gene_list, genes)

          genes2 <- geneList[idx_genes, ]

          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 (!quiet) {

            message(

              "Up-regulated genes in the ", second[i], " in ",

              first[i], " VS ", second[i], "\n"

            )

          }

          if (export) {

            write.csv(

              FinalDEGsU,

              file = paste0(

                "Up-regulated-", name, second[i], "in",

                first[i], "VS", second[i], ".csv"

              )

            )

            write.csv(

              FinalDEGsU,

              file = paste0(

                "Low-regulated-", name, first[i], "in",

                first[i], "VS", second[i], ".csv"

              )

            )

          }

          DEGsS <- c(DEGsS, FinalDEGsU[, 2])

          DEGsE <- c(DEGsE, as.character(FinalDEGsU[, 3]))

        }

      } else {

        DEGsTable[i, 1] <- paste0(first[i], " VS ", second[i])

        DEGsTable[i, 2] <- second[i]

        DEGsTable[i, 3] <- NA

        DEGsTable[i, 4] <- NA

        DEGsTable[i, 5] <- NA

        DEGsTable[i, 6] <- NA

        DEGsTable[i + com, 1] <- paste0(first[i], " VS ", second[i])

        DEGsTable[i + com, 2] <- first[i]

        DEGsTable[i + com, 3] <- NA

        DEGsTable[i + com, 4] <- NA

        DEGsTable[i + com, 5] <- NA

        DEGsTable[i + com, 6] <- NA

      }

    }

    if (!quiet & export) {

      message("The results of DEGs are saved in your directory")

    }

    colnames(DEGsTable) <- c(

      "Comparisons", "Target cluster", "Gene number", "File name",

      "Gene number", "File name"

    )

    if (export) write.csv(DEGsTable, file = paste0(filename_deg, ".csv"))

    if (!quiet) print(DEGsTable)

    sigDEG <- cbind(DEGsE, DEGsS)

    if (export) write.csv(sigDEG, file = paste0(filename_sigdeg, ".csv"))

    return(list(sigDEG, DEGsTable))

  }

)