# Eklund_script04.txt
# 2007-12-06
# Aron Charles Eklund
# http://www.cbs.dtu.dk/suppl/biascorr/
#
# Generate "Additional File 1".
# (intensity-dependent correlation plots for various
#  data sets and various normalization algorithms)

#  R 2.5.1

#  CM plots of various data sets and various normalization algorithms.

library(squash)
library(bias)
library(affy)
library(gcrma)

## various tumor data sets on Affy platform
ds <- c("bhattacharjee", "bild.brc", "gds360", "gse1456.u133a",
        "gse2603", "gse3494.u133a", "hess", "signoretti", "wang")         


## first compute correlation matrices, one data set at a time
## (making it easier to parallelize)

mbei <- function(x) expresso(x, normalize.method = "invariantset", 
               bg.correct = FALSE, pmcorrect.method = "pmonly",
               summary.method = "liwong") 

calcCM.meansort <- function(x) {
  ord <- order(rowMeans(x))
  calcCM( x[ord,] )
}  


############################ OLD VERSION
## for each data set, compute several CMs and save to a file.
severalCM.OLD <- function(dsName) {
  batchFileName <- paste('~/aron/chb/data/int/', dsName, '.batch.RData', sep = '')
  load(batchFileName)
  batchObjectName <- paste(dsName, '.batch', sep = '')
  cmObjectName <- paste(dsName, '.cm', sep = '')
  cmFileName <- paste(cmObjectName, '.RData', sep = '')
  batch <- updateObject(get(batchObjectName)) ## update old data
  cm.rma   <- calcCM.meansort(rma  (batch))
  cm.mas5  <- calcCM.meansort(mas5 (batch))
  cm.mbei  <- calcCM.meansort(mbei (batch))
  cm.gcrma <- calcCM.meansort(gcrma(batch))
  assign(cmObjectName, list(rma = cm.rma, mas5 = cm.mas5, 
                 mbei = cm.mbei, gcrma = cm.gcrma))
  save(list = cmObjectName, file = cmFileName)
}


## for each data set, compute several CMs and save to a file.
severalCM <- function(dsName) {
  batchFileName <- paste('~/aron/chb/data/int/', dsName, '.batch.RData', sep = '')
  load(batchFileName)
  batchObjectName <- paste(dsName, '.batch', sep = '')
  cmObjectName <- paste(dsName, '.cm', sep = '')
  cmFileName <- paste(cmObjectName, '.RData', sep = '')
  batch <- updateObject(get(batchObjectName)) ## update old data

  x.rma <- exprs(rma(batch))
  cm.rma   <- calcCM.meansort(x.rma)
  cm.rma.qn   <- calcCM.meansort(quantileNormalize(x.rma))

  x.mas5 <- exprs(mas5 (batch))
  cm.mas5  <- calcCM.meansort(x.mas5)
  cm.mas5.qn  <- calcCM.meansort(quantileNormalize (x.mas5))

  x.mbei <- exprs(mbei(batch))
  cm.mbei  <- calcCM.meansort(x.mbei)
  cm.mbei.qn  <- calcCM.meansort(quantileNormalize (x.mbei))

  x.gcrma <- exprs(gcrma(batch))
  cm.gcrma  <- calcCM.meansort(x.gcrma)
  cm.gcrma.qn  <- calcCM.meansort(quantileNormalize (x.gcrma))

  assign(cmObjectName, list(
    rma = cm.rma, rma.qn = cm.rma.qn, 
    mas5 = cm.mas5, mas5.qn = cm.mas5.qn, 
    mbei = cm.mbei, mbei.qn = cm.mbei.qn, 
    gcrma = cm.gcrma, gcrma.qn = cm.gcrma.qn ))

  save(list = cmObjectName, file = cmFileName)
}


for(i in 1:9) {
  severalCM(ds[i])
}


########### this could be the start a fresh session, or we could just continue

library(squash)

ds <- c("bhattacharjee", "bild.brc", "gds360", "gse1456.u133a",
        "gse2603", "gse3494.u133a", "hess", "signoretti", "wang")         

loadIntoList <- function(filenames) { 
  e <- new.env() 
  for (f in filenames) {
    load(f, envir = e)
  }
  rev(as.list(e))
} 

ds.cm <- paste(ds, '.cm.RData', sep = '')
all.cm <- loadIntoList(ds.cm)
names(all.cm) <- ds


DatasetLabels <- c('Bhattacharjee', 'Bild', 'Chang', 'Pawitan', 
  'Minn', 'Miller', 'Hess', 'Signoretti', 'Wang')

map <- colormap(c(-1,1), colFn = orangeblue, n = 10)

alg <- names(all.cm[[1]])
## [1] "rma"      "rma.qn"   "mas5"     "mas5.qn"  "mbei"     "mbei.qn"  "gcrma"   
## [8] "gcrma.qn"

algorithmLabels <- c('RMA', 'RMA + PSQN', 'MAS5', 'MAS5 + PSQN',
                   'MBEI', 'MBEI + PSQN', 'GCRMA', 'GCRMA + PSQN')

## plot one normalization algorithm per page.
pdf('AdditionalDataFile1.pdf', width = 8, height = 10)
par(mar = c(2, 2, 2, 2), xpd = NA, mfrow = c(5,4))
for (j in 1:4) {
  for (i in 1:9) {
    for (k in 0:1) {
      colorgram(all.cm[[i]][[alg[j + k]]], 
        map = map, key = FALSE,
        xlab = '', ylab = '', 
        main = paste(DatasetLabels[i], algorithmLabels[j + k]))
    }
  }
  plot(0, 0, type = 'n', axes = F, xlab='', ylab='') # placeholder
  opa <- par(mar = c(2, 6, 2, 8), yaxs='r', las = 1)
  drawKey(map, lab='median correlation', line = 4)
  par(opa)
}
dev.off()