extract_DEGs.R

#' @title Extract DEGs from a Limma fit.
#'
#' @description The function extracts the differentially expressed genes from an existing
#' limma fit and for a predefined analysis (aka a coefficient). It saves the data in a TSV file
#' and creates several plots in the process (Venn diagrams, MD plots).
#'
#' @param fit A fit as provided by a complete Limma analysis (for instance using
#' ArrayUtils::run_limma).
#' @param limma_coeffs A list of all limma coefficient names associated with that fit.
#' It will be used to name the output files, and to identify the results in the Venn
#' diagrams.
#' @param k The index of the current coefficient (from limma_coeffs).
#' @param output_data_dir A string representing the folder in which the output files
#' are going to be stored.
#' @param file_suffix A string used as a suffix for the file names. Default to "".
#' @param file_prefix A string used to prefix the file names. Default to "".
#' @param pval_adjust_method A string code indicating the multiple testing correction
#' method to use. Default to BH.
#' @param verbose A boolean representing whether the function should display log information. This
#' is FALSE by default.
#' @return NULL
extract_DEGs <- function(fit, limma_coeffs, k, output_data_dir,
                         file_suffix        = "",
                         file_prefix        = "",
                         pval_adjust_method = "BH",
                         verbose            = FALSE) {

  # We create the output file names.
  if (file_suffix != "") {
    file_suffix <- paste0(file_suffix, "_")
  }
  if (file_prefix != "") {
    file_prefix <- paste0( "_", file_prefix)
  }
  venn_fn  <- paste0(output_data_dir, file_suffix, "venn_",
                     limma_coeffs[k], file_prefix, ".png")
  md_fn    <- paste0(output_data_dir, file_suffix, "MD_",
                     limma_coeffs[k], file_prefix, ".png")
  table_fn <- paste0(output_data_dir, file_suffix, "toptable_",
                     limma_coeffs[k], file_prefix, ".tsv")

  # We extract the DEGs using toptable.
  table <- limma::topTable(fit,
                           coef          = paste0("coeff", k),
                           n             = Inf,
                           adjust.method = pval_adjust_method,
                           confint       = TRUE)

  # Another way to extract the differentially expressed genes (in parallel to topTable).
  # This is only used to then create the Venn diagrams.
  results <- limma::decideTests(fit,
                                coefficient   = paste0("coeff", k),
                                adjust.method = pval_adjust_method)

  # We log information.
  if (verbose == TRUE) {
    message(paste0("[", Sys.time(), "] Differentially expressed genes extracted."))
  }

  # We plot Venn diagrams that summarize the results and show to which
  # extend the analyses overlap.
  grDevices::png(venn_fn)
  limma::vennDiagram(results,
                     names        = limma_coeffs,
                     include      = c("up", "down"),
                     circle.col   = c("darkorange", "darkorange2", "darkorange4"),
                     counts.col   = c("blue", "green"),
                     lwd          = 3,
                     cex          = c(1.0, 0.1, 0.1),
                     show.include = TRUE)
  grDevices::dev.off()

  # We clean and log information.
  rm(venn_fn)
  if (verbose == TRUE) {
    message(paste0("[", Sys.time(), "] Venn diagram created."))
  }

  # We then plot the Mean vs Deviation of the comparisons and highlight the significant genes.
  grDevices::png(md_fn)
  limma::plotMD(fit,
                coef   = paste0("coeff", k),
                status = results[, k],
                values = c(1, -1),
                hl.col = c("red", "blue"),
                hl.pch = c(16, 16),
                hl.cex = c(0.9, 0.9))
  grDevices::dev.off()

  # We clean and log information.
  rm(md_fn, results)
  if (verbose == TRUE) {
    message(paste0("[", Sys.time(), "] MD plot created."))
  }

  # Save the top tables to TSV files.
  utils::write.table(table, file = table_fn, sep = "\t", quote = FALSE, col.names = NA)

  # We clean and log information.
  rm(table, table_fn)
  if (verbose == TRUE) {
    message(paste0("[", Sys.time(), "] Differentially expressed gene list saved."))
  }
}