1 # MvA plotting for microarray data
4 # Copyright 2003 By Don Armstrong
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.
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.
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
21 # $Id: mva_plotting.R,v 1.5 2006/06/10 04:14:45 don Exp $
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){
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)
34 if (NROW(array[,1]) > 1e4) {
38 mva.layout <- grid.layout(width,height)
39 pushViewport(viewport(layout=mva.layout))
42 par(mfcol=c(width,height),ps=9,cex.main=1)
44 # we could also use combn here instead
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,]
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) {
62 grid.text(title,x=unit(0.0,"npc")+unit(1,"lines"),y=unit(0.5,"npc"),just="center",rot=90)
64 pushViewport(viewport(layout=mva.layout))
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,
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)
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))
89 grid.text(paste(colnames(array)[a],colnames(array)[b],sep=" vs. "),y=unit(0.9,"npc"),x=unit(0.5,"npc"),just="center")
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)),
105 grid.text(title,x=unit(0.0,"npc")+unit(1,"lines"),y=unit(0.5,"npc"),just="center",rot=90)
109 title(title,outer=TRUE)