R/MAplotDE.R

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

#' MA-plot of DE features
#'
#' MA-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
#' @return A MA-plot
#' @author Marie-Agnes Dillies and Hugo Varet

# created Feb 14th, 2012
# modified April 26th (edgeR compatibility)
# modified Oct 30th, 2012 (png)
# modified April 4th, 2013 (graphDir)
# modified Sept 20th, 2013 (adapted to DESeq2)
# modified Oct 16th, 2013 (added adjMethod argument)
# modified Oct 28th, 2013 (first argument: res instead of dds)
# modified Nov 8th, 2013 (res is now a list)
# modified Nov 15th, 2013 (first argument is now complete for compatibility with edgeR)
# modified Mar 21st, 2014 (removed outputfile argument)
# modified July 31th, 2014 (removed adjMethod argument)
# modified Aug 5th, 2014 (removed graphDir argument)
# modified August 26th, 2019 (ggplot2)

MAplotDE <- function(complete, pvalCutoff = 0.05, out = TRUE, versionName="."){
  if (out) pdf(file=paste0("figures/", versionName, "-MAplot.pdf"))
  for (name in names(complete)){
    complete.name <- complete[[name]]
    complete.name <- complete.name[complete.name$baseMean>0,]
    complete.name$padj <- ifelse(is.na(complete.name$padj), 1, complete.name$padj)
    py <- complete.name$log2FoldChange
    q <- quantile(abs(py[is.finite(py)]), probs=0.99)
    complete.name$log2FoldChange[which(py > q)] <- q
    complete.name$log2FoldChange[which(py < -q)] <- -q
    complete.name$outfield <- factor(ifelse(py > q, "top", ifelse(py < -q, "bottom", "in")), 
                                     levels=c("bottom", "in", "top"))
    complete.name$padj <- ifelse(is.na(complete.name$padj), 1, complete.name$padj)
    complete.name$DE <- factor(ifelse(complete.name$padj <= pvalCutoff, "yes", "no"), levels=c("no", "yes"))
    p <- ggplot(data=complete.name, 
                aes(x=.data$baseMean, y=.data$log2FoldChange, color=.data$DE, fill=.data$DE, shape=.data$outfield)) +
      scale_x_continuous(trans = log10_trans(),
                         breaks = trans_breaks("log10", function(x) 10^x),
                         labels = trans_format("log10", math_format(~10^.x))) +
      geom_point(show.legend=FALSE, alpha=0.5, size=0.8) +
      scale_colour_manual(values=c("no"="black", "yes"="red"), drop=FALSE) +
      scale_shape_manual(values=c("bottom"=25, "in"=21, "top"=24), drop=FALSE) +
      scale_fill_manual(values=c("no"="black", "yes"="red"), drop=FALSE) +
      scale_y_continuous(expand=expansion(mult=c(0.03, 0.03))) +
      xlab("Mean of normalized counts") +
      ylab(expression(log[2]~fold~change)) +
      ggtitle(paste0(versionName," - MA-plot\nComparison ", name))
    print(p)
  }
  if (out) dev.off()
}