## require R devel (R 2.6.0)
source('http://www.bioconductor.org/biocLite.R')
biocLite('hgu133acdf', 'ALLMLL', 'arrayQualityMetrics', 'CCl4', 'limma', 'simpleaffy', 'affyPLM', 'preprocessCore', 'gcrma', 'matchprobes', 'affydata', 'RColorBrewer', 'vsn', 'affy', 'affyio', 'genefilter', 'survival', 'geneplotter', 'lattice', 'annotate', 'Biobase', 'fortunes')

#####################################################
### chunk number 1: load library
#####################################################
library("arrayQualityMetrics")

#####################################################
### chunk number 2: data import
#####################################################
library("ALLMLL")
data("MLL.A")
MLL.A
class(MLL.A)

#####################################################
### chunk number 3: report on AffyBatch
#####################################################
arrayQualityMetrics(expressionset = MLL.A,
                    outdir = "MLL",
                    force = FALSE,
                    do.logtransform = TRUE,
                    split.plots = 10)

#####################################################
### chunk number 4: AffyBatch normalization
#####################################################
rMLL.A = rma(MLL.A)
class(rMLL.A)
show(rMLL.A)

#####################################################
### chunk number 5: report on normalized AffyBatch
#####################################################
arrayQualityMetrics(expressionset = rMLL.A,
                    outdir = "rMLL",
                    force = FALSE,
                    do.logtransform = FALSE,
                    split.plots = 10)

#####################################################
### chunk number 6: load libraries
#####################################################
library("Biobase")
library("limma")
library("CCl4")
library("matchprobes")

#####################################################
### chunk number 7: path
#####################################################
datapath = system.file("extdata", package="CCl4")
dir(datapath)

#####################################################
### chunk number 8: RGList
#####################################################
p = read.AnnotatedDataFrame("samplesInfo.txt", path=datapath)
p
CCl4_RGList = read.maimages(files=sampleNames(p), 
   path = datapath, 
   source = "genepix", 
   columns = list(R = 'F635 Median', Rb = 'B635 Median', 
                  G = 'F532 Median', Gb = 'B532 Median'))

#####################################################
### chunk number 9: coordinates
#####################################################
 CCl4_RGList$genes[95:105,]

#####################################################
### chunk number 10: GC content
#####################################################
seq = read.AnnotatedDataFrame("013162_D_SequenceList_20060815.txt", 
path=datapath)
if(any(duplicated(featureNames(seq))))
  cat("IDs of the sequence file are not unique \n")
bc = basecontent(seq$Sequence)
GC = ((bc[,"C"]+bc[,"G"])/rowSums(bc))*100
mt = match(featureNames(seq), CCl4_RGList$genes$ID)
stopifnot(!any(is.na(mt)))
fData = cbind(CCl4_RGList$genes,GC=rep(as.numeric("NA"),
nrow(CCl4_RGList$genes)))
fData$GC[mt] = GC

#####################################################
### chunk number 11: reporter mapping
#####################################################
fData$hasTarget = (regexpr("^NM", CCl4_RGList$genes$Name) > 0)

#####################################################
### chunk number 12: NChannelSet
#####################################################
featureData = new("AnnotatedDataFrame", data = fData)

assayData = with(CCl4_RGList, assayDataNew(R=R, G=G, Rb=Rb, Gb=Gb))

varMetadata(p)$channel=factor(c("G", "R", "G", "R"), 
                    levels=c(ls(assayData), "_ALL"))

CCl4cfd <- new("NChannelSet",
              assayData = assayData,
              featureData = featureData,
              phenoData = p)

#####################################################
### chunk number 13: NChannelSet normalization
#####################################################
library("vsn")
nCCl4 = justvsn(CCl4cfd, subsample=2000)
save(nCCl4, file = "nCCl4.RData")

#####################################################
### chunk number 14: report on normalized NCHannelSet
#####################################################
arrayQualityMetrics(expressionset = nCCl4,
                    outdir = "CCl4")