R/plotDispEstimates.R

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

#' Scatter plots representing dispersion estimates vs mean expression
#'
#' Scatter plots representing dispersion estimates vs mean expression
#'
#' @param dds a \code{DESeqDataSet} object
#' @param out \code{TRUE} to export the figure
#' @param versionName versionName of the project
#' @return A plot of the mean-dispersion relationship
#' @author Marie-Agnes Dillies and Hugo Varet

# created Feb 14th, 2012
# modified Oct 30th, 2012 (png)
# modified Jan 16th, 2013 (pdf)
# modified Sept 20th, 2013 (adapted for DESeq2)
# modified Mar 21st, 2014 (removed outputfile argument)
# modified Aug 5th, 2014 (removed graphDir argument)
# modified August 26th, 2019 (ggplot2)

plotDispEstimates <- function(dds, out = TRUE, versionName="."){
  if (out) pdf(file=paste0("figures/", versionName, "-dispersionEstimates.pdf"), width=9)
  d <- as.data.frame(mcols(dds)[,c("baseMean", "dispGeneEst", "dispFit", "dispersion")])
  d <- d[which(d$baseMean > 0),]
  d <- data.frame(baseMean=rep(d$baseMean, 3),
                  value=c(d$dispGeneEst, d$dispFit, d$dispersion),
                  variable=factor(rep(c("dispGeneEst", "dispFit", "dispersion"), each=nrow(d)),
                                  levels=c("dispGeneEst", "dispFit", "dispersion")))
  print(ggplot(d, aes(x=.data$baseMean, y=.data$value, colour=.data$variable)) + 
    geom_point(size=0.1) +
      scale_x_continuous(trans = log10_trans(),
                         breaks = trans_breaks("log10", function(x) 10^x),
                         labels = trans_format("log10", math_format())) +
      scale_y_continuous(trans = log10_trans(),
                         breaks = trans_breaks("log10", function(x) 10^x),
                         labels = trans_format("log10", math_format())) +
    ylab("Dispersion") + 
    xlab("Mean of normalized counts") +
    scale_colour_manual(
      values=c("Black", "#e41a1c", "#377eb8"), 
      breaks=c("dispGeneEst", "dispFit", "dispersion"), 
      labels=c("Estimate", "Fit", "Final"),
      name="") +
    guides(colour = guide_legend(override.aes = list(size=2))) +
    ggtitle(paste(versionName, "Dispersions", sep=" - ")))
  if (out) dev.off()
}