R/volcanoPlotDE.R

In biomics-pasteur-fr/RNADiff: R Tools for RNA-Seq Projects at Institut Pasteur (Biomics Platform)

#' volcano-plot of DE features
#'
#' volcano-plot of DE features
#'
#' @param complete list of results from \code{exportComplete()} or \code{exportComplete.edgeR()}
#' @param pvalCutoff cut-off to apply on each adjusted p-value
#' @param out \code{TRUE} to export the figure
#' @param versionName versionName of the project
#' @param padjlim top limit of the plot (and to replace circles by triangles) (`NULL`` by default)
#' @param baseMeanThreshold to display only gene with a baseMean (average expression intensity) higher than `baseMeanThreshold`
#' @return A volcano-plot
#' @author Hugo Varet

# created Sept 22nd, 2019

volcanoPlotDE <- function(complete, pvalCutoff=0.05, out=TRUE, versionName=".", 
                          padjlim=NULL, baseMeanThreshold=0){
  if (out) pdf(file=paste0("figures/", versionName, "-volcanoPlot.pdf"))
  for (name in names(complete)){
    complete.name <- complete[[name]]
    complete.name$padj[which(complete.name$padj==0)] <- .Machine$double.xmin
    complete.name <- complete.name[which(!is.na(complete.name$padj)),]
    complete.name <- complete.name[which(complete.name$baseMean >= baseMeanThreshold),]
    complete.name$DE <- factor(ifelse(complete.name$padj <= pvalCutoff, "yes", "no"), levels=c("no", "yes"))
    if (is.null(padjlim)) padjlim.name <- quantile(complete.name$padj, probs=0.01, na.rm=TRUE) else padjlim.name <- padjlim
    complete.name$outfield <- factor(ifelse(complete.name$padj < padjlim.name, "top", "in"), levels=c("in", "top"))
    complete.name$padj[which(complete.name$padj < padjlim.name)] <- padjlim.name
    reverselog_trans <- function(base = exp(1)) {
      trans <- function(x) -log(x, base)
      inv <- function(x) base^(-x)
      trans_new(paste0("reverselog-", format(base)), trans, inv,
                log_breaks(base = base),
                domain = c(.Machine$double.xmin, Inf))
    }
    p <- ggplot(data=complete.name, 
                aes(x=.data$log2FoldChange, y=.data$padj, color=.data$DE, shape=.data$outfield)) +
      geom_point(show.legend=FALSE, alpha=0.5) +
      scale_y_continuous(trans = reverselog_trans(10),
                         breaks = trans_breaks("log10", function(x) 10^x),
                         labels = trans_format("log10", math_format(~10^.x))) +
      scale_colour_manual(values=c("no"="black", "yes"="red"), drop=FALSE) +
      scale_shape_manual(values=c("in"=16, "top"=17), drop=FALSE) +
      xlab(expression(log[2]~fold~change)) +
      ylab("Adjusted P-value") +
      ggtitle(paste0("Volcano plot - ", gsub("_", " ", name)))
    print(p)
  }
  if (out) dev.off()
}