R/QQplot.R

In DanteR: Proteomics data analysis tool

# Written by Ashoka D. Polpitiya
# for the Department of Energy (PNNL, Richland, WA)
# Copyright 2007, Battelle Memorial Institute
# E-mail: ashoka.polpitiya@pnl.gov
# Website: http://omics.pnl.gov/software
# -------------------------------------------------------------------------
#
# Licensed under the Apache License, Version 2.0; you may not use this file except
# in compliance with the License.  You may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0
#
# R Plotting functions used in DAnTE
# -------------------------------------------------------------------------


QQplot.dialog <- list(
  Data.dataframeItem = "T_Data", label="Data Selection",
    signal = c("default", "get.colnames", "Data.Columns"),
  Data.Columns.variableSelectorItem=character(0), label="Select Columns", 
  BREAK=T,
  reference.radiobuttonItem = c(value=1,2,3,4), item.labels=c("Normal", "Exponential", "Student", "Weibull"), label = "Reference Distribution",
    signal = c("default", function(item, exprate, degfree, wshape, wscale){
       sapply(c(exprate, degfree, wshape, wscale), gtkWidgetSetSensitive, FALSE)
       litems <- list(list(), list(exprate), list(degfree), list(wshape, wscale))
       sapply(litems[[get.value(item)]], gtkWidgetSetSensitive, TRUE)
      }, "exprate", "degfree", "wshape", "wscale"),
  exprate.numericItem=4, label = "Exponential: rate", indent=10,
  degfree.numericItem=4, label = "Student: degrees of freedom", indent=10,  
  wshape.numericItem=2, label = "Weibull: rate", indent=10,
  wscale.numericItem=1, label = "Weibull: scale", indent=10
)                    

QQplot <- function(Data, 
                      Data.Columns,
                      reference = "Normal",
                      wshape = 2,
                      wscale = 1,
                      degfree = 4,
                      exprate=1,
                      ncols = 2,
                      file="deleteme.png",
                      bkground="white",
                      colF="#ffc38a",
                      colB="#5FAE27",
                      colL="#FF0000",
                      ...)
{
    # Plot histograms, the distribution profile and the reference profile

    #png(filename=file,width=1152,height=864,pointsize=12,
    #        bg=bkground,res=600)
    ##require(Cairo)
    ##CairoPNG(filename=file,width=IMGwidth,height=IMGheight,pointsize=FNTsize,bg=bkground)
   data <- Data[,Data.Columns]
   ncols <- ceiling(sqrt(length(Data.Columns)))
    if (ncols == 0) ncols <- 1
    m <- ceiling((NCOL(data))/ncols)
    par(mfrow=c(m,ncols), cex=.6, mex=.6, oma=c(5, 2, 2, 2), mar=c(4,5,5,1))

#    tryCatch(
#    {
        for (i in 1:NCOL(data))
        {
          if (NCOL(data) == 1)
          {
            par(mfrow=c(1,1))
            xx <- data
            qtitle = "Normal Q-Q Plot"
          }
          else
          {
            xx <- data[,i]
            qtitle = colnames(data)[i]
          }

          #qqnorm(xx, col=colB, bg=colF, pch=21, main=qtitle)
          #qqline(xx, col = colL)
          qqplot.1(xx, reference = reference, shape = wshape,
                scale = wscale, def = degfree, exprate=exprate,
                colb=colB, bgrnd=colF,
                psize=21, main=qtitle)
          qqline.1(xx, reference = reference, colL = colL, shape = wshape,
                scale = wscale, def = degfree, exprate=exprate)
        }
#    },
#    interrupt = function(ex)
#    {
#      cat("An interrupt was detected.\n");
#      print(ex);
#    },
#    error = function(ex)
#    {
#      plot(c(1,1),type="n",axes=F,xlab="",ylab="")
#      text(1.5,1,paste("Error:", ex),cex=2)
#      cat("An error was detected.\n");
#      print(ex);
#    },
#    finally =
#    {
#      #cat("Releasing tempfile...");
#      par(mfrow=c(1,1),pch=1)
#      ##dev.off()
#      #cat("done\n");
#    }) # tryCatch()
  return(recordPlot())      
}

#-------------------------------------------------------------------
## Generic qqplots for any distribution
## (Normal, Exponential, Student and Weibull)
qqplot.1 <- function (data, position = 0.5, reference = "Normal",
                shape=2, scale=1, def = 4, exprate = 1,
                main='QQ plot', colb="#5FAE27", bgrnd="#ffc38a", psize=21, ...)
{
    data <- data[!is.na(data)]
    n <- length(data)
    plot.points <- ppoints(n, position)
    xpoints <- switch(reference,
      Normal = qnorm(plot.points, ...),
      Exponential = qexp(plot.points, rate=exprate, ...),
      Student = qt(plot.points, df=def, ...),
      Weibull = qweibull(plot.points, shape, scale, ...)
    )
    plot(xpoints, sort(data), col = colb, bg = bgrnd, pch = 21,
        xlab = paste("Theoretical Quantiles - ", reference),
        ylab = "Sample Quantiles", main = main)
}

#-------------------------------------------------------------------
## qqline.1 is based on qqline with additional, optional arguments.
qqline.1 <- function(x, reference = "Normal", exprate = 1,
              shape=2, scale=1, def = 4, colL='red', ...) {
  x <- x[!is.na(x)]
  plot.points <- c(0.25, 0.75)
  data.quartiles <- quantile(x, plot.points, na.rm = T)
  norm.quartiles <- switch(reference,
      Normal = qnorm(plot.points, ...),
      Exponential = qexp(plot.points, rate=exprate, ...),
      Student = qt(plot.points, df=def, ...),
      Weibull = qweibull(plot.points, shape, scale, ...)
  )
  b <- (data.quartiles[2] - data.quartiles[1]) /
    (norm.quartiles[2] - norm.quartiles[1])
  a <- data.quartiles[1] - norm.quartiles[1] * b
  abline(a, b, col=colL, ...)
  ans <- as.vector(c(a, b))
  names(ans) <- c("intercept", "slope")
  invisible(ans)
}