R/plot_gsea.R

In BIGslu/BIGpicture: Standard Plots for BIGslu

#' Plot GSEA normalized enrichment
#'
#' @param gsea Data frame output by SEARchways::BIGsea including pathway, FDR, and NES
#' @param fdr_cutoff Numeric. Maximum FDR to plot. Default is 0.2
#' @param fdr_colors Numeric vector. Cutoffs for color groups. Default is c(0.01, 0.05, 0.1, 0.2)
#' @param show_overlap Logical if should show overlap across all facets even if some missing (TRUE) or give each facet it's own axis labels (FALSE). Default is TRUE
#'
#' @param fdr.cutoff Deprecated form of fdr_cutoff
#' @param fdr.colors Deprecated form of fdr_colors
#' @param show.overlap Deprecated form of show_overlap
#'
#' @return ggplot2 object
#' @export
#'
#' @examples
#' library(SEARchways)
#' library(dplyr)
#' #Get fold change information from example model
#' genes.FC <- example.model$lmerel %>%
#'             filter(variable %in% c("virus", "asthma")) %>%
#'             select(variable, gene, estimate)
#' #Run GSEA
#' example_gsea <- SEARchways::BIGsea(gene_df = genes.FC, category = "H", ID = "ENSEMBL")
#'
#' #Plot
#' plot_gsea(example_gsea, fdr_cutoff = 0.5,
#'          fdr_colors = c(0.1, 0.5))

plot_gsea <- function(gsea, fdr_cutoff = 0.2,
                      fdr_colors = c(0.01, 0.05, 0.1, 0.2),
                      show_overlap = TRUE,
                      #Deprecated
                      fdr.cutoff = NULL, fdr.colors = NULL, show.overlap = NULL
                      ){

  FDR <-NES<-Significance<-pathway<-group<- NULL

  # Back compatibility
  if(!is.null(fdr.cutoff)){fdr_cutoff <- fdr.cutoff}
  if(!is.null(fdr.colors)){fdr_colors <- fdr.colors}
  if(!is.null(show.overlap)){show_overlap <- show.overlap}

  #### Format data ####
  dat.signif <- gsea %>%
    dplyr::filter(FDR < fdr_cutoff)
  #keep nonsignif overlap if requested
  if(show_overlap){
    dat.format <- gsea %>%
      dplyr::filter(pathway %in% dat.signif$pathway) %>%
      dplyr::mutate(pathway = gsub("_", " ", pathway)) %>%
      #Add 0 enrichment values
      tidyr::complete(group, pathway) %>%
      dplyr::mutate(NES = ifelse(is.na(NES), 0, NES),
                    FDR = ifelse(is.na(FDR), 1, FDR))
  } else{
    dat.format <- dat.signif %>%
      dplyr::mutate(pathway = gsub("_", " ", pathway))
  }

  if(nrow(dat.format) == 0){stop("No gene sets are significant. Please increase fdr_cutoff.")}

  fdr_colors.sort <- sort(fdr_colors)
  dat.format$Significance <- NA
  for(i in 1:length(fdr_colors.sort)){
    if(i==1){
      dat.format$Significance[dat.format$FDR < fdr_colors.sort[i]] <- paste("FDR <", fdr_colors.sort[i])
    } else{
      dat.format$Significance[dat.format$FDR < fdr_colors.sort[i] & dat.format$FDR >= fdr_colors.sort[i-1]] <- paste("FDR <", fdr_colors.sort[i])
    }
  }

  #Enrichment score limits
  plot.lim <- max(abs(dat.format$NES))+0.1

  #### Plot ####
  p1 <- dat.format %>%
    ggplot2::ggplot(ggplot2::aes(stats::reorder(pathway, NES), NES)) +
    ggplot2::geom_segment(ggplot2::aes(stats::reorder(pathway, NES),
                     xend=pathway, y=0, yend=NES)) +
    ggplot2::geom_point(size=3, ggplot2::aes(fill = Significance),
               shape=21, stroke=1) +
    ggplot2::geom_hline(yintercept = 0) +

    ggplot2::scale_fill_brewer(palette = "RdYlBu", na.value="grey70") +
    ggplot2::lims(y=c(-plot.lim,plot.lim)) +
    ggplot2::coord_flip() +
    ggplot2::labs(x="", y="Normalized Enrichment Score") +
    ggplot2::theme_bw()

  if(show_overlap & length(unique(dat.format$group)) > 1){
    p1.facet <- p1 + ggplot2::facet_grid( ~ group)
  } else if(length(unique(dat.format$group)) > 1){
    p1.facet <- p1 + ggplot2::facet_wrap( ~ group, scales="free")
  } else{
    p1.facet <- p1
  }
  return(p1.facet)
}