###################################################
### chunk number 1: queryAE
###################################################
library("ArrayExpress")
sets = queryAE(keywords = "pneumonia", species = "homo+sapiens")


###################################################
### chunk number 2: ArrayExpress-raw
###################################################
rawset = ArrayExpress("E-MEXP-1422")


###################################################
### chunk number 3: ArrayExpress-columnsneeded
###################################################
eset = try(ArrayExpress("E-DKFZ-1"))


###################################################
### chunk number 4: ArrayExpress-withcolumns
###################################################
eset = ArrayExpress("E-DKFZ-1",
  rawcol=list(R="Software Unknown:Foreground red",
    G="Software Unknown:Foreground green",
    Rb="Software Unknown:Background red" ,
    Gb="Software Unknown:Background green"))


###################################################
### chunk number 5: getAE-full
###################################################
mexp1422 = getAE("E-MEXP-1422", type = "full")


###################################################
### chunk number 6: magetab2bioc-full
###################################################
rawset= magetab2bioc(files = mexp1422)


###################################################
### chunk number 7: getcolproc
###################################################
cn = getcolproc(mexp1422)
show(cn)


###################################################
### chunk number 8: procset
###################################################
proset = procset(mexp1422, cn[2])


###################################################
### chunk number 9: import
###################################################
AEset = ArrayExpress("E-MEXP-1416")


###################################################
### chunk number 10: norm
###################################################
library("affy")
AEsetnorm = rma(AEset)


###################################################
### chunk number 11: fac
###################################################
fac = grep("Factor.Value",colnames(pData(AEset)), value=T)


###################################################
### chunk number 12: qanorm eval=FALSE
###################################################
## if (suppressWarnings(require("arrayQualityMetrics", quietly=TRUE))) {
##   qanorm = arrayQualityMetrics(AEsetnorm, 
##     outdir = "QAnorm", 
##     intgroup = fac, 
##     grouprep = TRUE)
## }


###################################################
### chunk number 13: limma
###################################################
library("limma")
facs =  pData(AEset)[,fac]
facs[facs[,2]=="female",2]="F"
facs[facs[,2]=="male",2]="M"
facs[facs[,1]=="Parkinson disease",1]="parkinson"
facs = paste(facs[,1],facs[,2], sep=".")
f = factor(facs)
design = model.matrix(~0+f)
colnames(design) = levels(f)
fit = lmFit(AEsetnorm, design)
cont.matrix = makeContrasts(normal.FvsM = normal.F-normal.M,             
    parkinson.FvsM = parkinson.F-parkinson.M,  
    Diff=(parkinson.F-parkinson.M)-(normal.F-normal.M),
    levels=design)
fit2 = contrasts.fit(fit, cont.matrix)
fit2 = eBayes(fit2)
res = topTable(fit2, coef = "parkinson.FvsM", adjust = "BH")