R/plotEvolution.r

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

#' Plots of the expression of some interesting features
#'
#' Plots of the expression of some interesting features
#'
#' @param mat a \code{matrix} containing expression measures
#' @param features a character vector of interesting features
#' @param target the target of the project
#' @param varInt1 first variable of interest (displayed on the x-axis)
#' @param varInt2 optional second variable of interest (used to colour points)
#' @param colors colors of the plot
#' @param plot "perGene" to plot each gene profile, "aggregated" to average over the features or "both"
#' @param ylab label for the y-axis
#' @param out \code{TRUE} to export the figure
#' @param versionName versionName of the project
#' @return A plot of the interesting features
#' @author Hugo Varet

# created Jan 15th, 2016
# modified February 23rd, 2017

plotEvolution <- function(mat, features, target, varInt1, varInt2=NULL, colors=c("red","blue","black"), 
                          plot=c("perGene", "aggregated", "both"), ylab="Expression", out=TRUE, versionName="."){
  plot <- plot[1]
  varAbs <- target[,varInt1]
  if (!is.null(varInt2)) varOrd <- target[,varInt2]
  x <- as.numeric(as.factor(varAbs))
  col <- colors[as.numeric(as.factor(if(!is.null(varInt2)){varOrd}else{varAbs}))]
  
  if (plot %in% c("perGene", "both")){
    if (out) pdf(file=paste0("figures/", versionName, "-plotEvolution.pdf"))
    for (i in 1:length(features)){
      y <- mat[features[i],]
      plot(x, y, col=col, main=features[i], pch=16, xaxt="n", ylim=range(y), ylab=ylab, xlab="")
      axis(1, at=unique(x), labels=levels(varAbs))
      legend("topleft",legend=if(!is.null(varInt2)){levels(varOrd)}else{levels(varAbs)},
             col=colors[1:nlevels(if(!is.null(varInt2)){varOrd}else{varAbs})],pch=16)
      # add lines between points
      if (!is.null(varInt2)){
        ymean <- tapply(y, list(varAbs,varOrd), mean)
        for (i in 1:nlevels(varOrd)){
          lines(1:nlevels(varAbs), ymean[,levels(varOrd)[i]], col=colors[i], lwd=1)
        }
      } else{
        lines(1:nlevels(varAbs), tapply(y, list(x=varAbs), mean), lwd=1)
      }
    }
    if (out) dev.off()
  }
  
  if (plot %in% c("aggregated", "both")){
    if (out) pdf(file=paste0("figures/", versionName, "-plotEvolutionAggregated.pdf"))
    y <- colMeans(matrix(mat[features,], nrow=length(features)))
    plot(x, y, col=col, main=paste(length(features), "features aggregated"), 
         pch=16, xaxt="n", ylim=range(y), ylab=ylab, xlab="")
    axis(1, at=unique(x), labels=levels(varAbs))
    legend("topleft",legend=if(!is.null(varInt2)){levels(varOrd)}else{levels(varAbs)},
           col=colors[1:nlevels(if(!is.null(varInt2)){varOrd}else{varAbs})],pch=16)
    # add lines between points
    if (!is.null(varInt2)){
      ymean <- tapply(y, list(varAbs,varOrd), mean)
      for (i in 1:nlevels(varOrd)){
        lines(1:nlevels(varAbs), ymean[,levels(varOrd)[i]], col=colors[i], lwd=1)
      }
    } else{
      lines(1:nlevels(varAbs), tapply(y, list(x=varAbs), mean), lwd=1)
    }
    if (out) dev.off()
  }
}