R/pairwiseScatterPlots.R

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

#' Scatter plots for pairwise comparaisons of log counts
#'
#' Scatter plots for pairwise comparaisons of log counts
#'
#' @param counts raw counts
#' @param group vector of the condition from which each sample belongs
#' @param out \code{TRUE} to export the figure
#' @param maxCol maximum number of samples to produce the "big" plot with all the crossings
#' @param versionName versionName of the project
#' @return A pairwise scatter plot
#' @author Marie-Agnes Dillies and Hugo Varet

# created Feb 21th, 2012
# modified Sept 27th, 2012 (pdf output file)
# modified Oct 30th, 2012 (png)
# modified Jan 16th, 2013 (pdf)
# modified Oct 25th, 2013 (modifications for multiple factors)
# modified Oct 29th, 2013 (added maxCol argument)
# modified Feb 12th, 2014 (condition: 30000 rows)
# modified Mar 05th, 2014 (condition 30000 rows also in the dimensions of the output file)
# modified Mar 06th, 2014 (added abline() to the plots)
# modified Mar 19th, 2014 (fixed a bug in the dimensions of the graphs produced)
# modified Mar 21st, 2014 (removed outputfile argument)
# modified Aug 5th, 2014 (removed graphDir argument)
# modified August 26th, 2019 (ggplot2)

pairwiseScatterPlots <- function(counts, group, out=TRUE, maxCol=8, versionName="."){

  counts <- removeNul(counts)
  ncol <- ncol(counts)
  separated_graphs <- I(ncol>maxCol)
  if (out){
    pdf(file=paste0("figures/", versionName, "-pairwiseScatter.pdf"),
        width=ifelse(separated_graphs, 5.5, 5*ncol),
        height=ifelse(separated_graphs, 5.5, 4*ncol))
  }

  group <- data.frame(group=apply(group, 1, paste, collapse="-"))
  d <- data.frame(counts+1)
  
  if (separated_graphs){
    for (i in 1:(ncol-1)){
      for (j in (i+1):ncol){
        print(ggplot(data=cbind(d, z=1), aes_string(x=names(d)[i], y=names(d)[j], z="z")) +
                stat_summary_2d(fun=function(z) log(sum(z)), bins=60, show.legend=FALSE) +
                scale_x_continuous(trans = log10_trans(),
                                   breaks = trans_breaks("log10", function(x) 10^x),
                                   labels = trans_format("log10", math_format(~10^.x))) +
                scale_y_continuous(trans = log10_trans(),
                                   breaks = trans_breaks("log10", function(x) 10^x),
                                   labels = trans_format("log10", math_format(~10^.x))) +
                xlab(names(d)[i]) +
                ylab(names(d)[j]) +
                ggtitle(paste(versionName," - Pairwise scatter plot\n", colnames(counts)[i]," - ", colnames(counts)[j], sep="")) +
                geom_abline(slope=1, intercept=0, linetype="dashed", col="lightgrey"))
      }
    }
  } else{
    p <- list()
    layout_matrix = matrix(NA, ncol=ncol, nrow=ncol)
    k <- 1
    for (i in 1:ncol){
      for (j in 1:ncol){
        if (i==j) next
        p[[k]] <- ggplot(data=cbind(d, z=1), aes_string(x=names(d)[i], y=names(d)[j], z="z")) +
          stat_summary_2d(fun=function(z) log(sum(z)), bins=60, show.legend=FALSE) +
          scale_x_continuous(trans = log10_trans(),
                             breaks = trans_breaks("log10", function(x) 10^x),
                             labels = trans_format("log10", math_format(~10^.x))) +
          scale_y_continuous(trans = log10_trans(),
                             breaks = trans_breaks("log10", function(x) 10^x),
                             labels = trans_format("log10", math_format(~10^.x))) +
          xlab(names(d)[i]) +
          ylab(names(d)[j]) +
          geom_abline(slope=1, intercept=0, linetype="dashed", col="lightgrey")
        layout_matrix[j,i] <- k
        k <- k+1
      }
    }
    tmpfun <- function(...) grid.arrange(..., ncol=ncol, nrow=ncol, layout_matrix=layout_matrix,
                                         top=textGrob(paste(versionName," - Pairwise scatter plot"), x=0.01, just="left", gp=gpar(fontsize=20)))
    do.call(tmpfun, p)
  }
  
  if (out) dev.off()
}