R/7_quick_double_enrich.R

In tinyarray: Expression Data Analysis and Visualization

##' quick_enrich
##'
##' Perform enrichment analysis for a gene list; human is the default species.
##'
##' @param genes a gene symbol or Entrez ID vector
##' @param kkgo_file optional RDS cache filename for KEGG and GO results;
##'   when `NULL`, the function recomputes on every call and does not write a cache
##' @param destdir destination directory used to save the cache file when
##'   `kkgo_file` is provided
##' @inheritParams trans_ensembl_exp
##' @return enrichment results and dotplots
##' @author Xiaojie Sun
##' @importFrom ggplot2 facet_grid
##' @export
##' @examples
##' \dontrun{
##' head(genes)
##' g = quick_enrich(genes,destdir = tempdir())
##' names(g)
##' g$kk
##' g$go
##' }
##' @seealso
##' \code{\link{double_enrich}}

quick_enrich <- function(genes,
                         kkgo_file = NULL,
                         destdir = getwd(),
                         species = "human"){
  if(any(is.na(suppressWarnings(as.numeric(genes))))){
    org = .tinyarray_orgdb_for_species(species)
    .tinyarray_require_clusterProfiler()
    s2e <- clusterProfiler::bitr(genes, fromType = "SYMBOL",
                toType = "ENTREZID",
                OrgDb = org$orgdb)
    s2e <- .tinyarray_dedup_by_key(s2e, "SYMBOL")
    genes = s2e$ENTREZID
  }
  compute_quick_enrich <- function(genes, species) {
    org = .tinyarray_orgdb_for_species(species)
    orgdb = org$orgdb
    .tinyarray_require_clusterProfiler()
    kk <- tryCatch(
      clusterProfiler::enrichKEGG(
        gene = genes,
        organism = org$kegg,
        pvalueCutoff = 0.05
      ),
      error = function(e) {
        message("KEGG enrichment failed: ", conditionMessage(e))
        NULL
      }
    )
    if (!is.null(kk)) {
      kk <- .tinyarray_make_readable_enrich(kk, OrgDb = orgdb, keyType = "ENTREZID")
      if (nrow(kk@result) == 0L) {
        message("No KEGG pathway enriched")
        kk <- NULL
      }
    } else {
      message("No KEGG pathway enriched")
    }
    go <- clusterProfiler::enrichGO(genes,
                                    OrgDb = orgdb,
                                    ont="all",
                                    readable = FALSE)
    go <- .tinyarray_make_readable_enrich(go, OrgDb = orgdb, keyType = "ENTREZID")
    if (!is.null(go) && nrow(go@result) == 0L) {
      message("No GO terms enriched")
      go <- NULL
    }
    list(kk = kk, go = go)
  }
  if (is.null(kkgo_file) || !length(kkgo_file) || !nzchar(kkgo_file[1L])) {
    obj = compute_quick_enrich(genes, species)
  } else {
    f = file.path(destdir,kkgo_file)
    legacy_f = if (grepl("\\.rds$", f, ignore.case = TRUE)) sub("\\.rds$", ".Rdata", f, ignore.case = TRUE) else f
    if(file.exists(f)){
      obj = .tinyarray_read_cache(f)
    }else if(legacy_f != f && file.exists(legacy_f)){
      obj = .tinyarray_read_cache(legacy_f)
      saveRDS(obj, file = f)
    }else{
      obj = compute_quick_enrich(genes, species)
      saveRDS(obj, file = f)
    }
  }
  kk = obj$kk
  go = obj$go
  k1 = if (is.null(kk)) 0L else sum(kk@result$p.adjust <0.05, na.rm = TRUE)
  k2 = if (is.null(go)) 0L else sum(go@result$p.adjust <0.05, na.rm = TRUE)
  if (k1 == 0L) {
    kk.dot = "no pathway enriched"
  } else{
    kk.dot = clusterProfiler::dotplot(kk)
  }

  if (k2 == 0L) {
    go.dot = "no terms enriched"
  } else{
    go.dot = clusterProfiler::dotplot(go, split="ONTOLOGY",font.size =10,showCategory = 5)+
      facet_grid(ONTOLOGY~., scales="free")
  }
  result = list(kk = kk,go = go,kk.dot = kk.dot,go.dot = go.dot)
  return(result)
}

##' draw enrichment bar plots for both up and down genes
##'
##' draw enrichment bar plots for both up and down genes,for human only.
##'
##' @param deg a data.frame contains at least two columns:"ENTREZID" and "change"
##' @param n how many terms will you perform for up and down genes respectively
##' @param color color for bar plot
##' @inheritParams quick_enrich
##' @return a list with kegg and go bar plot according to up and down genes enrichment result.
##' @author Xiaojie Sun
##' @importFrom stringr str_to_lower
##' @importFrom stringr str_wrap
##' @importFrom dplyr mutate
##' @importFrom dplyr arrange
##' @importFrom ggplot2 ggplot
##' @importFrom ggplot2 geom_bar
##' @importFrom ggplot2 coord_flip
##' @importFrom ggplot2 theme_light
##' @importFrom ggplot2 ylim
##' @importFrom ggplot2 scale_x_discrete
##' @importFrom ggplot2 scale_y_continuous
##' @importFrom ggplot2 theme
##' @export
##' @examples
##' \dontrun{
##' deg = data.frame(
##'   ENTREZID = genes[1:20],
##'   change = rep(c("up", "down"), 10)
##' )
##' double_enrich(deg)
##' }
##' @seealso
##' \code{\link{quick_enrich}}

double_enrich <- function(deg,n = 10,color = c("#2874C5", "#f87669"),species = "human"){
  deg$change = tolower(deg$change)
  up = quick_enrich(deg$ENTREZID[deg$change=="up"],"up.rds",destdir = tempdir(),species = species)
  down = quick_enrich(deg$ENTREZID[deg$change=="down"],"down.rds",destdir = tempdir(),species = species)
  ud_enrich = function(df){
    if (is.null(df) || !nrow(df)) {
      return(NULL)
    }
    padj = pmax(as.numeric(df$p.adjust), .Machine$double.xmin)
    df$pl = ifelse(df$change == "up",-log10(padj),log10(padj))
    df = df[order(df$change, df$pl), , drop = FALSE]
    df$Description = factor(df$Description,levels = unique(df$Description),ordered = TRUE)
    lm = if (any(is.finite(df$pl))) {
      c(floor(min(df$pl, na.rm = TRUE)), ceiling(max(df$pl, na.rm = TRUE)))
    } else {
      c(-1, 1)
    }
    ggplot(df, aes(x=Description, y= pl)) +
      geom_bar(stat='identity', aes(fill=change), width=.7)+
      scale_fill_manual(values = color)+
      coord_flip()+
      theme_light() +
      ylim(lm)+
      scale_x_discrete(labels=function(x) vapply(strwrap(x, width = 30), paste, character(1), collapse = "\n"))+
      theme(
        panel.border = element_blank()
      )
  }
  kk_ok = !is.null(up$kk) & !is.null(down$kk)
  go_ok = !is.null(up$go) & !is.null(down$go)
  if (kk_ok || go_ok) {
    result = list()
    take_top_terms = function(x) {
      if (is.null(x)) {
        return(NULL)
      }
      if (!nrow(x)) {
        return(x[0, , drop = FALSE])
      }
      x[seq_len(min(n, nrow(x))), , drop = FALSE]
    }
    if (kk_ok) {
      up$kk@result$change = "up"
      down$kk@result$change = "down"
      kk = rbind(take_top_terms(up$kk@result), take_top_terms(down$kk@result))
      result$kp = tryCatch(
        ud_enrich(kk),
        error = function(e) {
          warning("KEGG plot failed: ", conditionMessage(e))
          "KEGG plot unavailable"
        }
      )
    } else {
      warning("No KEGG pathway enriched; showing GO result only")
    }
    if (go_ok) {
      up$go@result$change = "up"
      down$go@result$change = "down"
      go = rbind(take_top_terms(up$go@result), take_top_terms(down$go@result))
      result$gp = tryCatch(
        ud_enrich(go),
        error = function(e) {
          warning("GO plot failed: ", conditionMessage(e))
          "GO plot unavailable"
        }
      )
    } else {
      warning("No GO term enriched; showing KEGG result only")
    }
  } else {
    warning("No KEGG or GO pathway was enriched; returning results from quick_enrich")
    result = list(up = up,
                  down = down)
  }
  class(result) <- c("tinyarray_double_enrich", class(result))
  return(result)

}

##' @export
print.tinyarray_double_enrich <- function(x, ...) {
  cat("<tinyarray_double_enrich>\n")
  if (length(x)) {
    cat("components:", paste(names(x), collapse = ", "), "\n")
  }
  invisible(x)
}

utils::globalVariables(c("change","pl","Description"))