# 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)
  }
}