R/plotStats.R

In SwathXtend: SWATH extended library generation and statistical data analysis

###########################################################################
#' Plot the statistical summary of the base and external libraries
#' @param datBaseLib a data frame for base library
#' @param datExtLib a data frame for external/addon library
#' @param wb an openxlsx workbook object
#' @param sheet A name or index of a worksheet
#' @param ... arguments to pass onto cleanLib function
#' @return .png files of statistical plots of the two libraries 
#' @examples 
#' libfiles <- paste(system.file("files",package="SwathXtend"),c("Lib2.txt","Lib3.txt"),sep="/")
#' datBaseLib <- readLibFile(libfiles[1])
#' datExtLib <- readLibFile(libfiles[2])
#' plotStats(datBaseLib, datExtLib)
############################################################################ 

plotStats <- function(datBaseLib, datExtLib, wb=NULL, sheet=NULL, ...)
{
 
  if(!missing(datExtLib)){ ## 2 libs
    
    ## protein peptide numbers
    nums.pro <- sapply(list(datBaseLib$uniprot_id,datExtLib$uniprot_id),
                       FUN=function(x){length(unique(x))})
    names(nums.pro) <- c("Base", "Addon")
                       
    nums.pep <- sapply(list(datBaseLib$stripped_sequence,datExtLib$stripped_sequence),
                       FUN=function(x){length(unique(x))})
  
    names(nums.pep) <- c("Base", "Addon")    
  
    png("protein_peptide_numbers.png")
    layout(matrix(1:2,nrow=1))
    barplot(nums.pro, col=c("darkblue","red"),beside=TRUE,main="Proteins")
    barplot(nums.pep, col=c("darkblue","red"),beside=TRUE,main="Peptides")
    dev.off()
    
    if(!is.null(wb) & !is.null(sheet))
      insertImage(wb, sheet, "protein_peptide_numbers.png", width=6, height=5, startRow=1, startCol=1)
    
      ## Venn diagrams
  
    #grid.newpage()  
     
    if(nums.pro[1]>0 && nums.pro[2]>0){
      
      png("proteinVennDigram.png")
  
      
      draw.pairwise.venn(nums.pro[1], nums.pro[2], 
                         length(proteinOverlap(datBaseLib,datExtLib, ...)), 
                         category = c("Base Library Proteins", "External Library Proteins"), 
                         lty = rep("blank", 2), 
                         fill = c("light blue", "pink"), 
                         alpha = rep(0.5, 2), cat.pos = c(0,20), 
                         cat.dist = rep(0.025, 2))
      
      dev.off()
      
      if(!is.null(wb) & !is.null(sheet))
        insertImage(wb, sheet, "proteinVennDigram.png", width=5, height=5, startRow=1, startCol=15)
      
    }
      
    if(nums.pep[1]>0 && nums.pep[2]>0){
      
      png("peptideVennDigram.png")
      draw.pairwise.venn(nums.pep[1], nums.pep[2], 
                         length(unique(intersect(datBaseLib$stripped_sequence,
                                                 datExtLib$stripped_sequence))), 
                         category = c("Base Library Peptides", "External Library Peptides"), 
                         lty = rep("blank", 2), 
                         fill = c("light blue", "pink"), 
                         alpha = rep(0.5, 2), cat.pos = c(0,20), 
                         cat.dist = rep(0.025, 2))
      
      
      dev.off()
      
      if(!is.null(wb) & !is.null(sheet))
        insertImage(wb, sheet, "peptideVennDigram.png", width=5, height=5, startRow=1, startCol=25)
    }
 
    # Density plots

    png("densityPlots.png")
    layout(matrix(1:4,ncol=2))
    
    plot(density(datBaseLib$confidence), main="Base lib conf", xlab="Confidence")
    abline(v=0.99,col="red")
    
    plot(density(datExtLib$confidence), main="Addon lib conf", xlab="Confidence")
    abline(v=0.99, col="red")
    
    plot(density(log(datBaseLib$relative_intensity)), main="Base lib ion intensity", xlab="log ion intensity")
    abline(v=1.6, col="red") #log5
    
    plot(density(log(datExtLib$relative_intensity)), main="Addon lib ion intensity", xlab="log ion intensity")
    abline(v=1.6, col="red")
    
    
    dev.off()
    
    if(!is.null(wb) & !is.null(sheet))
      insertImage(wb, sheet, "densityPlots.png", width=6, height=6, startRow=30, startCol=1)
    
    
    
    # histgrams of peptides/protein and ions/peptide
    png("histgram_peptide_ion.png")
    layout(matrix(1:4,ncol=2))
    
    x<-aggregate(datBaseLib$uniprot_id, by=list(Protein=datBaseLib$uniprot_id, Peptide=datBaseLib$stripped_sequence), 
                 FUN=function(x){length(x)})
    y<-aggregate(x$Peptide, by=list(Protein=x$Protein), FUN=function(x){length(x)})
    
    
    hist(x[,3], xlab="Number of ions per peptide", main="Base lib")
    hist(y[,2], breaks="Sturges",xlab="Number of peptides per protein", main="Base lib")
    
    x<-aggregate(datExtLib$uniprot_id, by=list(Protein=datExtLib$uniprot_id, Peptide=datExtLib$stripped_sequence), 
                 FUN=function(x){length(x)})
    y<-aggregate(x$Peptide, by=list(Protein=x$Protein), FUN=function(x){length(x)})
    
    
    hist(x[,3], xlab="Number of ions per peptide", main="Addon lib")    
    
    hist(y[,2], breaks=5, xlab="Number of peptides per protein", main="Addon lib", freq=FALSE)
    
    dev.off()
    
    if(!is.null(wb) & !is.null(sheet))
      insertImage(wb, sheet, "histgram_peptide_ion.png", width=6, height=6, startRow=30, startCol=15)
    
 
  } else { # single lib
    nums.pro <- length(unique(datBaseLib$uniprot_id)) 
    
    nums.pep <- length(unique(datBaseLib$stripped_sequence))
    
    
    png("protein_peptide_numbers.png")
    par(mfrow=c(1,2))
    barplot(nums.pro,col="lightblue",main="Proteins",ylim = c(0,max(nums.pro)))
    barplot(nums.pep,col="lightseagreen",main="Peptides",ylim = c(0,max(nums.pep)))
    dev.off()
  
    if(!is.null(wb) & !is.null(sheet))
      insertImage(wb, sheet, "protein_peptide_numbers.png", width=6, height=5, startRow=1, startCol=1)

    png("densityPlots.png")
    layout(matrix(1:2,ncol=2))
    
    plot(density(datBaseLib$confidence), main="Base lib conf", xlab="Confidence")
    abline(v=0.99,col="red")
    
    plot(density(log(datBaseLib$relative_intensity)), main="Ion intensity", xlab="log ion intensity")
    abline(v=1.6, col="red") #log5
    
    dev.off()
    
    if(!is.null(wb) & !is.null(sheet))
      insertImage(wb, sheet, "densityPlots.png", width=6, height=6, startRow=30, startCol=1)
    
  }
}