R/mva_plotting.R

   1 # MvA plotting for microarray data
   2
   3
   4 # Copyright 2003 By Don Armstrong
   5
   6 # This program is free software; you can redistribute it and/or modify
   7 # it under the terms of the GNU General Public License as published by
   8 # the Free Software Foundation; either version 2 of the License, or
   9 # (at your option) any later version.
  10
  11 # This program is distributed in the hope that it will be useful, but
  12 # WITHOUT ANY WARRANTY; without even the implied warranty of
  13 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 # General Public License for more details.
  15
  16 # You should have received a copy of the GNU General Public License
  17 # along with this program; if not, write to the Free Software
  18 # Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  19 # 02111-1307, USA.
  20
  21 # $Id: mva_plotting.R,v 1.5 2006/06/10 04:14:45 don Exp $
  22
  23
  24 # mva.plot
  25 # Plots a pairwise M versus A plot with a loess fit line
  26 mva.plot.pair <- function(array,title=F,small.positive=0.1){
  27   size <- ncol(array)
  28   npair <- size*(size-1)/2
  29   height <- min(round((npair*3/4)^0.5),3)
  30   width <- min(ceiling(npair/height),3)
  31   ## print(c("width",width,"height",height))
  32   def.par <- par(no.readonly = TRUE)
  33   use.grid <- 0
  34   if (NROW(array[,1]) > 1e4) {
  35     use.grid <- 1
  36   }
  37   if (use.grid) {
  38     mva.layout <- grid.layout(width,height)
  39     pushViewport(viewport(layout=mva.layout))
  40   }
  41   else {
  42     par(mfcol=c(width,height),ps=9,cex.main=1)
  43   }
  44   # we could also use combn here instead
  45   num.plot <- 0
  46   for(a in 1:(size-1)){
  47     for(b in (a+1):size){
  48 #      par(mfg=c(a,b-1))
  49       num.plot <- num.plot + 1
  50       array[array[,a]<=0,a] <- small.positive
  51       array[array[,b]<=0,b] <- small.positive
  52       A <- 0.5*log2(array[,a]*array[,b])
  53       M <- log2(array[,a]/array[,b])
  54       temp.frame <- data.frame(cbind(A,M))
  55       temp.frame <- temp.frame[abs(temp.frame$M) < 7.8 & temp.frame$A < 15.8 & temp.frame$A > 0,]
  56       if (use.grid) {
  57         ## print(c("num.plot",num.plot))
  58         ## print(c("row",floor((num.plot -1)/ width) %% height))
  59         ## print(c("col",(num.plot-1) %% width))
  60         if (num.plot > 1 && floor((num.plot -1)/ width) %% height == 0 && (num.plot-1) %% width == 0) {
  61           popViewport(1)
  62           grid.text(title,x=unit(0.0,"npc")+unit(1,"lines"),y=unit(0.5,"npc"),just="center",rot=90)
  63           grid.newpage()
  64           pushViewport(viewport(layout=mva.layout))
  65         }
  66         ## print(floor((num.plot -1)/ width) %% height+1)
  67         ## print((num.plot - 1)%% width+1)
  68         pushViewport(viewport(layout.pos.row=(floor((num.plot -1)/ width) %% height)+1, layout.pos.col=((num.plot - 1)%% width)+1))
  69         hb <- hexbin(temp.frame,
  70                      ##xlab="A",ylab="M",
  71                      ##ybnds = c(-8,8),
  72                      ##xbnds = c(0,16)
  73                      );
  74         # P <- plot(hb,type="n",
  75         #           ## main=paste(colnames(array)[a],colnames(array)[b],sep=" vs. "),
  76         #           xlab="",ylab="",newpage=FALSE,legend=FALSE)
  77         hb.viewport <- hexViewport(hb,mar=unit(0.1,"npc"))
  78         pushHexport(hb.viewport)
  79         grid.hexagons(hb,style="colorscale",border=gray(.1),pen=gray(.6),minarea=.1,maxarea=1.5,colramp=BTY)
  80         grid.xaxis()
  81         grid.yaxis()
  82         temp.loess <- loess.smooth(A,M,family="gaussian",
  83                                 control=loess.control(surface="interpolate",statistics="approximate",
  84                                   trace.hat="approximate",iterations=2,cell=0.6))
  85         panel.loess(A,M,
  86                     col="red",lwd=2)
  87         popViewport(1)
  88         ## grid.rect()
  89         grid.text(paste(colnames(array)[a],colnames(array)[b],sep=" vs. "),y=unit(0.9,"npc"),x=unit(0.5,"npc"),just="center")
  90         popViewport(1)
  91 #        grid.lines()
  92       }
  93       else {
  94         plot(A,M,pch=".",main=paste(colnames(array)[a],colnames(array)[b],sep=" vs. "),xlab="A",ylab="M",ylim=c(-8,8),xlim=c(0,16))
  95         lines(loess.smooth(A,M,family="gaussian",
  96                            control=loess.control(surface="interpolate",statistics="approximate",
  97                              trace.hat="approximate",iterations=2,cell=0.6)),
  98               col="blue",lwd=2)
  99       }
 100       ## abline(h=0)
 101     }
 102   }
 103   if (use.grid) {
 104     popViewport(1)
 105     grid.text(title,x=unit(0.0,"npc")+unit(1,"lines"),y=unit(0.5,"npc"),just="center",rot=90)
 106     ## grid.newpage()
 107   }
 108   else {
 109     title(title,outer=TRUE)
 110     par(def.par)
 111   }
 112 }