MAplot and volcanoPlot

DESCRIPTION

@@ -2,7 +2,7 @@ Package: SARTools
 Type: Package
 Title: Statistical Analysis of RNA-Seq Tools
 Version: 1.6.7
-Date: 2018-12-03
+Date: 2019-05-07
 Author: Marie-Agnes Dillies and Hugo Varet
 Maintainer: Hugo Varet <[email protected]>
 Depends: R (>= 3.3.0), DESeq2 (>= 1.12.0), edgeR (>= 3.12.0), xtable

NEWS

@@ -1,6 +1,7 @@
 CHANGES IN VERSION 1.6.7
 ------------------------
         o modified the installation guidelines for the Bioconductor packages
+	o new parameters in the MAPlot() and volcanoPlot() functions thanks to Ken Field suggestion on GitHub
 
 CHANGES IN VERSION 1.6.6
 ------------------------

R/MAPlot.R

@@ -5,10 +5,11 @@
 #' @param complete A \code{list} of \code{data.frame} containing features results (from \code{exportResults.DESeq2()} or \code{exportResults.edgeR()})
 #' @param alpha cut-off to apply on each adjusted p-value
 #' @param outfile TRUE to export the figure in a png file
+#' @param log2FClim numeric vector containing both upper and lower y-axis limits for all the MA-plots produced (NULL by default to set them automatically)
 #' @return A file named MAPlot.png in the figures directory containing one MA-plot per comparison
 #' @author Marie-Agnes Dillies and Hugo Varet
 
-MAPlot <- function(complete, alpha=0.05, outfile=TRUE){
+MAPlot <- function(complete, alpha=0.05, outfile=TRUE, log2FClim=NULL){
   ncol <- ifelse(length(complete)<=4, ceiling(sqrt(length(complete))), 3)
   nrow <- ceiling(length(complete)/ncol)
   if (outfile) png(filename="figures/MAPlot.png", width=cairoSizeWrapper(1800*ncol), height=cairoSizeWrapper(1800*nrow), res=300)
@@ -18,7 +19,11 @@ MAPlot <- function(complete, alpha=0.05, outfile=TRUE){
         complete.name <- complete.name[complete.name$baseMean>0,]
         complete.name$padj <- ifelse(is.na(complete.name$padj),1,complete.name$padj)
         log2FC <- complete.name$log2FoldChange
-        ylim <- 1.1 * c(-1,1) * quantile(abs(log2FC[is.finite(log2FC)]), probs=0.99)
+        if (is.null(log2FClim)){
+          ylim <- 1.1 * c(-1,1) * quantile(abs(log2FC[is.finite(log2FC)]), probs=0.99)
+        } else{
+          ylim <- log2FClim
+        }
         plot(complete.name$baseMean, pmax(ylim[1], pmin(ylim[2], log2FC)),
 	         log = "x", cex=0.45, las = 1, ylim = ylim,
              col = ifelse(complete.name[,"padj"] < alpha, "red", "black"),

R/summarizeResults.DESeq2.r

@@ -8,11 +8,14 @@
 #' @param cooksCutoff outliers detection threshold (TRUE to let DESeq2 choosing it or FALSE to disable the outliers detection)
 #' @param alpha significance threshold to apply to the adjusted p-values
 #' @param col colors for the plots
+#' @param log2FClim numeric vector containing both upper and lower y-axis limits for all the MA-plots produced (NULL by default to set them automatically)
+#' @param padjlim numeric value between 0 and 1 for the adjusted p-value upper limits for all the volcano plots produced (NULL by default to set them automatically)
 #' @return A list containing: (i) a list of \code{data.frames} from \code{exportResults.DESeq2()}, (ii) the table summarizing the independent filtering procedure and (iii) a table summarizing the number of differentially expressed features
 #' @author Hugo Varet
 
 summarizeResults.DESeq2 <- function(out.DESeq2, group, independentFiltering=TRUE, cooksCutoff=TRUE,
-                                    alpha=0.05, col=c("lightblue","orange","MediumVioletRed","SpringGreen")){
+                                    alpha=0.05, col=c("lightblue","orange","MediumVioletRed","SpringGreen"),
+                                    log2FClim=NULL, padjlim=NULL){
   # create the figures/tables directory if does not exist
   if (!I("figures" %in% dir())) dir.create("figures", showWarnings=FALSE)
   if (!I("tables" %in% dir())) dir.create("tables", showWarnings=FALSE)
@@ -50,10 +53,10 @@ summarizeResults.DESeq2 <- function(out.DESeq2, group, independentFiltering=TRUE
   rawpHist(complete=complete)
 
   # MA-plots
-  MAPlot(complete=complete, alpha=alpha)
+  MAPlot(complete=complete, alpha=alpha, log2FClim=log2FClim)
 
   # Volcano plots
-  volcanoPlot(complete=complete, alpha=alpha)
+  volcanoPlot(complete=complete, alpha=alpha, padjlim=padjlim)
 
   return(list(complete=complete, tabIndepFiltering=tabIndepFiltering, nDiffTotal=nDiffTotal))
 }

R/summarizeResults.edgeR.r

@@ -7,11 +7,14 @@
 #' @param counts matrix of raw counts
 #' @param alpha significance threshold to apply to the adjusted p-values
 #' @param col colors for the plots
+#' @param log2FClim numeric vector containing both upper and lower y-axis limits for all the MA-plots produced (NULL by default to set them automatically)
+#' @param padjlim numeric value between 0 and 1 for the adjusted p-value upper limits for all the volcano plots produced (NULL by default to set them automatically)
 #' @return A list containing: (i) a list of \code{data.frames} from \code{exportResults.edgeR()} and (ii) a table summarizing the number of differentially expressed features
 #' @author Hugo Varet
 
 summarizeResults.edgeR <- function(out.edgeR, group, counts, alpha=0.05,
-                                   col=c("lightblue","orange","MediumVioletRed","SpringGreen")){  
+                                   col=c("lightblue","orange","MediumVioletRed","SpringGreen"),
+                                   log2FClim=NULL, padjlim=NULL){  
   # create the figures/tables directory if does not exist
   if (!I("figures" %in% dir())) dir.create("figures", showWarnings=FALSE)
   if (!I("tables" %in% dir())) dir.create("tables", showWarnings=FALSE)
@@ -34,10 +37,10 @@ summarizeResults.edgeR <- function(out.edgeR, group, counts, alpha=0.05,
   rawpHist(complete=complete)
 
   # MA-plots
-  MAPlot(complete=complete, alpha=alpha)
+  MAPlot(complete=complete, alpha=alpha, log2FClim=log2FClim)
 
   # Volcano plots
-  volcanoPlot(complete=complete, alpha=alpha)
+  volcanoPlot(complete=complete, alpha=alpha, padjlim=padjlim)
 
   return(list(complete=complete, nDiffTotal=nDiffTotal))
 }

R/volcanoPlot.r

@@ -5,10 +5,11 @@
 #' @param complete A \code{list} of \code{data.frame} containing features results (from \code{exportResults.DESeq2()} or \code{exportResults.edgeR()})
 #' @param alpha cut-off to apply on each adjusted p-value
 #' @param outfile TRUE to export the figure in a png file
+#' @param padjlim numeric value between 0 and 1 for the adjusted p-value upper limits for all the volcano plots produced (NULL by default to set them automatically)
 #' @return A file named volcanoPlot.png in the figures directory containing one volcano plot per comparison
 #' @author Hugo Varet
 
-volcanoPlot <- function(complete, alpha=0.05, outfile=TRUE){
+volcanoPlot <- function(complete, alpha=0.05, outfile=TRUE, padjlim=NULL){
   ncol <- ifelse(length(complete)<=4, ceiling(sqrt(length(complete))), 3)
   nrow <- ceiling(length(complete)/ncol)
   if (outfile) png(filename="figures/volcanoPlot.png", width=cairoSizeWrapper(1800*ncol), height=cairoSizeWrapper(1800*nrow), res=300)
@@ -17,7 +18,11 @@ volcanoPlot <- function(complete, alpha=0.05, outfile=TRUE){
       complete.name <- complete[[name]]
       complete.name$padj[which(complete.name$padj==0)] <- .Machine$double.xmin
       log10pval <- -log10(complete.name$padj)
-      ylim <- c(0,1) * quantile(log10pval, probs=0.99, na.rm=TRUE)
+      if (is.null(padjlim)){
+        ylim <- c(0,1) * quantile(log10pval, probs=0.99, na.rm=TRUE)
+      } else{
+        ylim <- c(0,1) * -log10(padjlim)
+      }
       plot(complete.name$log2FoldChange, pmin(ylim[2], log10pval), ylim=ylim, las=1, cex=0.45,
            xlab=expression(log[2]~fold~change), ylab=expression(-log[10]~adjusted~P~value),
            col=ifelse(complete.name$padj <= alpha, "red", "black"), pch=ifelse(log10pval >= ylim[2], 2, 20),