###################################################
### chunk number 1: loadlib
###################################################
#library("ArrayExpress")
#AEset = ArrayExpress("E-MEXP-461")
#save(AEset,file="AEset.RData")


###################################################
### chunk number 2: loaddata
###################################################
load("AEset.RData")
class(AEset)


###################################################
### chunk number 3: exploringAEset
###################################################
AEset
colnames(pData(AEset))


###################################################
### chunk number 4: arrayQM1col
###################################################
fac = colnames(pData(AEset))[grep("Factor",colnames(pData(AEset)))]
library("arrayQualityMetrics")
arrayQualityMetrics(expressionset = AEset,
                    outdir = "MEXP461NoNormalisation",
                    force = FALSE,
                    do.logtransform = TRUE,
                    intgroup = fac[2])


###################################################
### chunk number 5: arrayQM1colnorm
###################################################
rAEset = rma(AEset)
arrayQualityMetrics(expressionset = rAEset,
                    outdir = "MEXP461rmaNormalisation",
                    force = FALSE,
                    intgroup = fac[2])


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


###################################################
### chunk number 7: datapath
###################################################
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: XYcoordinates
###################################################
CCl4_RGList$genes[95:105,]
CCl4_RGList$genes$X = CCl4_RGList$genes$Row
CCl4_RGList$genes$Y = CCl4_RGList$genes$Column


###################################################
### chunk number 10: GCcontent
###################################################
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: ProbeMapping
###################################################
fData$hasTarget = (regexpr("^NM", CCl4_RGList$genes$Name) > 0)


###################################################
### chunk number 12: CCl4-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_"))
CCl4 <- new("NChannelSet",
              assayData = assayData,
              featureData = featureData,
              phenoData = p)


###################################################
### chunk number 13: Covariate
###################################################
cond = paste(p$RIN.Cy3,p$RIN.Cy5,sep="/")
poor = grep(cond,pattern="2.5")
medium = grep(cond,pattern="^5/|/5")
good = grep(cond,pattern="9.7")
cov = rep(0, length = nrow(p))
cov[good] = "Good"
cov[medium] = "Medium"
cov[poor] = "Poor"
phenoData(CCl4)$RNAintegrity = cov


###################################################
### chunk number 14: vsn-NChannelSet
###################################################
library("vsn")
nCCl4 = justvsn(CCl4, subsample=2000)


###################################################
### chunk number 15: arrayQM2col
###################################################
arrayQualityMetrics(expressionset = nCCl4,
                    outdir = "CCl4",
                    force = FALSE,
                    intgroup = "RNAintegrity")


###################################################
### chunk number 16: pkgs
###################################################
toLatex(sessionInfo())