R/diagLogNormalityDisp.R

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

#' Check the log-normality of the features dispersion estimates
#'
#' Check the log-normality of the features dispersion estimates
#'
#' @param dds a \code{DESeqDataSet} object
#' @param out \code{TRUE} to export the figure
#' @param versionName versionName of the project
#' @return A hisogram of the dispersions with the log-normal density
#' @author Marie-Agnes Dillies and Hugo Varet

# created Oct 29th, 2013
# modified Mar 21st, 2014 (removed outputfile argument)
# modified August 26th, 2019 (ggplot2)

diagLogNormalityDisp <- function(dds, out = TRUE, versionName="."){
  if (out) pdf(file=paste0("figures/", versionName, "-diagLogNormalityDisp.pdf"))
  disp <- mcols(dds)$dispGeneEst
  disp <- disp[!is.na(disp)]
  disp <- disp[disp>1e-8]
  disp <- log(disp)
  mean.disp <- mean(disp,na.rm=TRUE)
  sd.disp <- sd(disp,na.rm=TRUE)
  d <- data.frame(disp)
  print(ggplot(data=d, aes(x=.data$disp)) +
          geom_histogram(bins=80, aes(y=.data$..density..)) +
          scale_y_continuous(expand=expansion(mult=c(0.01, 0.05))) +
          xlab("Feature dispersion estimate") +
          ylab("Density") +
          ggtitle(paste0(versionName, " - log-normality dispersion diagnostic")) +
          stat_function(fun = dnorm, args = list(mean = mean.disp, sd = sd.disp)))
  if (out) dev.off()
}