R/BCVPlot.R

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

#' BCV plot (for edgeR object)
#'
#' BCV plot (for edgeR object)
#'
#' @param dge a \code{DGEList} object
#' @param out \code{TRUE} to export the figure
#' @param versionName versionName of the project
#' @return A BCV plot
#' @author Marie-Agnes Dillies and Hugo Varet

# created Nov 14th, 2013
# modified Mar 21st, 2014 (removed outputfile argument)
# modified Aug 5th, 2014 (removed graphDir argument)

BCVPlot <- function(dge, out = TRUE, versionName="."){
  if (out) pdf(file=paste0("figures/", versionName, "-BCV.pdf"), width=9)
  A <- dge$AveLogCPM
  if (is.null(A)) A <- aveLogCPM(dge$counts, offset = getOffset(dge))
  disp <- getDispersion(dge)
  if (is.null(disp)) stop("No dispersions to plot")
  if (attr(disp, "type") == "common") disp <- rep(disp, length = length(A))
  d <- data.frame(A=A,
                  sqrtdisp=sqrt(disp),
                  sqrttagwise=sqrt(dge$tagwise.dispersion),
                  sqrttrended=sqrt(dge$trended.dispersion),
                  sqrtcommon=sqrt(dge$common.dispersion))
  print(ggplot() +
    scale_color_manual("", values=c("black", "blue", "red"), labels=c("Tagwise", "Trend", "Common")) +
    geom_point(data=d, mapping=aes(x=.data$A, y=.data$sqrttagwise, color="a"), size=0.5, alpha=0.5) +
    geom_line(data=d, mapping=aes(x=.data$A, y=.data$sqrttrended, color="b")) +
    geom_hline(data=d, aes(yintercept=.data$sqrtcommon, color="c")) +
    xlab("Average log CPM") +
    ylab("Biological coefficient of variation") +
    ggtitle(paste0(versionName, " - BCV plot")))
  if (out) dev.off()      
}