R/gseaView.R

In WubingZhang/ggView: Practical package for scientific plotting

#' GSEA plot that mimic the plot generated by broad institute's GSEA software
#'
#' @title gseaView
#' @param x gseaResult object
#' @param geneSetID gene set ID
#' @param title plot title
#' @param seg_height height of segments
#' @param base_size base font size
#' @return ggplot object.
#' @export
#' @import ggplot2
#' @author Wubing Zhang
#'
gseaView <- function(x, geneSetID, title = "", seg_height = 0.08,
                     base_size = 11) {

  if (length(geneSetID) == 1) {
    gsdata <- gsInfo(x, geneSetID)
  } else {
    gsdata <- do.call(rbind, lapply(geneSetID, gsInfo, object = x))
  }
  gsdata = gsdata[order(gsdata$pvalue), ]
  gsdata$Description = paste0(gsdata$Description, " (",
                              format(gsdata$pvalue, scientific = TRUE, digits = 3), ")")
  gsdata$Description = factor(gsdata$Description, levels = unique(gsdata$Description))
  i <- round(min(gsdata$runningScore)-0.02, 2)
  for (term in unique(gsdata$Description)) {
    idx <- which(gsdata$ymin != 0 & gsdata$Description == term)
    gsdata[idx, "ymin"] <- i-seg_height
    gsdata[idx, "ymax"] <- i
    i <- i-seg_height
  }
  p <- ggplot(gsdata, aes_(x = ~x, color= ~Description))
  p <- p + geom_line(aes_(y = ~runningScore), size=1)
  p <- p + geom_linerange(aes_(ymin=~ymin, ymax=~ymax))
  p <- p + labs(x = NULL, y = "Running Enrichment Score", title = title)
  p <- p + scale_y_continuous(breaks = seq(round(min(gsdata$runningScore),1), round(max(gsdata$runningScore),1), by = 0.2),
                              expand=c(0,0))
  p <- p + theme_classic(base_size = 14)
  p <- p + theme(legend.title = element_blank(),
          legend.background = element_rect(fill = "transparent"))
  p <- p + theme(plot.margin=margin(t=.2, r = .2, b=0, l=.2, unit="cm"))

  return(p)
}


#' extract gsea result of selected geneSet
#'
#'
#' @title gsInfo
#' @param object gseaResult object
#' @param geneSetID gene set ID
#' @return data.frame
#' @author Guangchuang Yu
gsInfo <- function(object, geneSetID) {
  geneList <- object@geneList

  if (is.numeric(geneSetID))
    geneSetID <- object@result[geneSetID, "ID"]

  geneSet <- object@geneSets[[geneSetID]]
  exponent <- object@params[["exponent"]]
  df <- gseaScores(geneList, geneSet, exponent, fortify=TRUE)
  df$ymin=0
  df$ymax=0
  pos <- df$position == 1
  h <- diff(range(df$runningScore))/20
  df$ymin[pos] <- -h
  df$ymax[pos] <- h
  df$geneList <- geneList

  df$Description <- object@result[geneSetID, "Description"]
  df$pvalue <- object@result[geneSetID, "pvalue"]
  df$p.adjust <- object@result[geneSetID, "p.adjust"]
  return(df)
}
gseaScores <- getFromNamespace("gseaScores", "DOSE")