add functions

[don_microarray.git] / R / mva_plotting.R

diff --git a/R/mva_plotting.R b/R/mva_plotting.R
new file mode 100644 (file)
index 0000000..f40be20
--- /dev/null
+++ b/R/mva_plotting.R
@@ -0,0 +1,112 @@
+# MvA plotting for microarray data
+
+
+# Copyright 2003 By Don Armstrong
+
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+# General Public License for more details.
+
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
+# 02111-1307, USA.
+
+# $Id: mva_plotting.R,v 1.5 2006/06/10 04:14:45 don Exp $
+
+
+# mva.plot
+# Plots a pairwise M versus A plot with a loess fit line
+mva.plot.pair <- function(array,title=F,small.positive=0.1){
+  size <- ncol(array)
+  npair <- size*(size-1)/2
+  height <- min(round((npair*3/4)^0.5),3)
+  width <- min(ceiling(npair/height),3)
+  ## print(c("width",width,"height",height))
+  def.par <- par(no.readonly = TRUE)
+  use.grid <- 0
+  if (NROW(array[,1]) > 1e4) {
+    use.grid <- 1
+  }
+  if (use.grid) {
+    mva.layout <- grid.layout(width,height)
+    pushViewport(viewport(layout=mva.layout))
+  }
+  else {
+    par(mfcol=c(width,height),ps=9,cex.main=1)
+  }
+  # we could also use combn here instead
+  num.plot <- 0
+  for(a in 1:(size-1)){
+    for(b in (a+1):size){
+#      par(mfg=c(a,b-1))
+      num.plot <- num.plot + 1
+      array[array[,a]<=0,a] <- small.positive
+      array[array[,b]<=0,b] <- small.positive
+      A <- 0.5*log2(array[,a]*array[,b])
+      M <- log2(array[,a]/array[,b])
+      temp.frame <- data.frame(cbind(A,M))
+      temp.frame <- temp.frame[abs(temp.frame$M) < 7.8 & temp.frame$A < 15.8 & temp.frame$A > 0,]
+      if (use.grid) {
+        ## print(c("num.plot",num.plot))
+        ## print(c("row",floor((num.plot -1)/ width) %% height))
+        ## print(c("col",(num.plot-1) %% width))
+        if (num.plot > 1 && floor((num.plot -1)/ width) %% height == 0 && (num.plot-1) %% width == 0) {
+          popViewport(1)
+          grid.text(title,x=unit(0.0,"npc")+unit(1,"lines"),y=unit(0.5,"npc"),just="center",rot=90)
+          grid.newpage()
+          pushViewport(viewport(layout=mva.layout))
+        }
+        ## print(floor((num.plot -1)/ width) %% height+1)
+        ## print((num.plot - 1)%% width+1)
+        pushViewport(viewport(layout.pos.row=(floor((num.plot -1)/ width) %% height)+1, layout.pos.col=((num.plot - 1)%% width)+1))
+        hb <- hexbin(temp.frame,
+                     ##xlab="A",ylab="M",
+                     ##ybnds = c(-8,8),
+                     ##xbnds = c(0,16)
+                     );
+        # P <- plot(hb,type="n",
+        #           ## main=paste(colnames(array)[a],colnames(array)[b],sep=" vs. "),
+        #           xlab="",ylab="",newpage=FALSE,legend=FALSE)
+        hb.viewport <- hexViewport(hb,mar=unit(0.1,"npc"))
+        pushHexport(hb.viewport)
+        grid.hexagons(hb,style="colorscale",border=gray(.1),pen=gray(.6),minarea=.1,maxarea=1.5,colramp=BTY)
+        grid.xaxis()
+        grid.yaxis()
+        temp.loess <- loess.smooth(A,M,family="gaussian",
+                                control=loess.control(surface="interpolate",statistics="approximate",
+                                  trace.hat="approximate",iterations=2,cell=0.6))
+        panel.loess(A,M,
+                    col="red",lwd=2)
+        popViewport(1)
+        ## grid.rect()
+        grid.text(paste(colnames(array)[a],colnames(array)[b],sep=" vs. "),y=unit(0.9,"npc"),x=unit(0.5,"npc"),just="center")
+        popViewport(1)
+#        grid.lines()
+      }
+      else {
+        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))
+        lines(loess.smooth(A,M,family="gaussian",
+                           control=loess.control(surface="interpolate",statistics="approximate",
+                             trace.hat="approximate",iterations=2,cell=0.6)),
+              col="blue",lwd=2)
+      }
+      ## abline(h=0)
+    }
+  }
+  if (use.grid) {
+    popViewport(1)
+    grid.text(title,x=unit(0.0,"npc")+unit(1,"lines"),y=unit(0.5,"npc"),just="center",rot=90)
+    ## grid.newpage()
+  }
+  else {
+    title(title,outer=TRUE)
+    par(def.par)
+  }
+}