###################################################
### chunk number 1: 
###################################################
options(width = 70)


###################################################
### chunk number 2:  eval=FALSE
###################################################
## library("affycoretools")
## pd <- new("AnnotatedDataFrame", 
##           data = read.table("pdata.txt", header = TRUE, row.names = 1))
## eset <- affystart(groups = rep(1:4, each = 3), 
##                   groupnames = unique(paste(pData(pd)[,1], 
##                   pData(pd)[,2], sep = "-")), 
##                   phenoData = pd)


###################################################
### chunk number 3: 
###################################################
library(affycoretools)
load("abatch.Rdata")
load("exprSet.Rdata")


###################################################
### chunk number 4: 
###################################################
plotHist(dat)


###################################################
### chunk number 5: 
###################################################
plotDeg(dat)


###################################################
### chunk number 6: 
###################################################
pd <- new("AnnotatedDataFrame", 
          data = read.table("pdata.txt", header = TRUE, row.names = 1))
sampleNames(pd) <- sampleNames(eset)
plotPCA(eset, groups = rep(1:4, each = 3), 
        groupnames = unique(paste(pData(pd)[,1], pData(pd)[,2], sep = "-")))
phenoData(eset) <- pd


###################################################
### chunk number 7: 
###################################################
plotPCA(eset, screeplot = TRUE)


###################################################
### chunk number 8:  eval=FALSE
###################################################
## eset1 <- affystart(filenames = list.celfiles()[1:6], 
##                    plot = FALSE, pca = FALSE)
## eset2 <- affystart(filenames = list.celfiles()[7:12], 
##                    plot = FALSE, pca = FALSE)
## eset <- new("ExpressionSet", 
##             exprs = cbind(exprs(eset1), exprs(eset2)), 
##             phenoData = pd,
##             annotation = annotation(eset1))


###################################################
### chunk number 9: 
###################################################
library(genefilter)
f1 <- kOverA(3, 6)
filt <- filterfun(f1)
index <- genefilter(eset, filt)
eset <- eset[index,]


###################################################
### chunk number 10: 
###################################################
library(limma)
grps <- paste(pData(eset)[,1], 
              pData(eset)[,2], sep = ".")
design <- model.matrix(~ 0 + factor(grps))
colnames(design) <- levels(factor(grps))
ugrps <- unique(grps)
contrasts <- matrix(c(1, -1, 0, 0, 0, 0, 1, -1),
                    ncol = 2, dimnames = list(ugrps,
                              paste(ugrps[c(1,3)], ugrps[c(2,4)],
                                    sep = " - ")))
fit <- lmFit(eset, design)
fit2 <- contrasts.fit(fit, contrasts)
fit2 <- eBayes(fit2)


###################################################
### chunk number 11: 
###################################################
design
contrasts


###################################################
### chunk number 12: 
###################################################
rslt <- decideTests(fit2)
vc <- vennCounts2(rslt, method = "same")
vennDiagram(vc, cex = 0.8)


###################################################
### chunk number 13: 
###################################################
vennDiagram(vc, cex = 0.6)


###################################################
### chunk number 14:  eval=FALSE
###################################################
## vennSelect(eset, design, rslt, contrasts, fit2)


###################################################
### chunk number 15:  eval=FALSE
###################################################
## limma2annaffy(eset, fit2, design, 
##               contrasts, annotation(eset), 
##               pfilt = 0.05)


###################################################
### chunk number 16:  eval=FALSE
###################################################
## index1 <- vennSelect(x = rslt, indices.only = TRUE)[[3]]
## probids <- unique(getLL(featureNames(eset)[index1], 
##                         annotation(eset)))
## univ <- unique(getLL(featureNames(eset),
##                      annotation(eset)))
## params <- new("GOHyperGParams", geneIds = probids,
##               universeGeneIds = univ, 
##               annotation = annotation(eset),
##               conditional = TRUE, ontology = "MF")
## hyp <- hyperGTest(params)
## htmlReport(hyp, file = "GO MF terms.html",
##            categorySize = 10)