###################################################
### chunk number 1: loadaffy
###################################################
library("affy")


###################################################
### chunk number 2: ReadAffy eval=FALSE
###################################################
## MyData <- ReadAffy()


###################################################
### chunk number 3: ReadAffy2 eval=FALSE
###################################################
## ReadAffy(filenames=filenames.of.your.data)


###################################################
### chunk number 4: readdata
###################################################
library("SpikeInSubset")
data(spikein95)


###################################################
### chunk number 5: addPhenoData
###################################################
pd <- data.frame(population=factor(c("Pop1", "Pop1", "Pop1", 
      "Pop2", "Pop2", "Pop2")),
       replicate=c(1, 2, 3, 1, 2, 3))
rownames(pd) <- sampleNames(spikein95)
vl <- data.frame(labelDescription= c("Pop1 is control, Pop2 is treatment",
           "1, 2, 3, arbitrary numbering"),
            row.names = c("population", "replicate"))
phenoData(spikein95) <- new("AnnotatedDataFrame", data=pd, varMetadata=vl)


###################################################
### chunk number 6: qareport eval=FALSE
###################################################
## 
##  library("affyQCReport")
## ##we need to fix up the sample names
## ##until the affyQAReport function is fixed
##  s2 = spikein95
##  sampleNames(s2) = sub("_", ".", sampleNames(s2))
##  qarep = affyQAReport(s2)
## 


###################################################
### chunk number 7: rma
###################################################
eset <- rma(spikein95)


###################################################
### chunk number 8: justrma eval=FALSE
###################################################
## eset<- just.rma(filenames=list.celfiles())


###################################################
### chunk number 9: exprsDemo
###################################################
e <- exprs(eset)
dim(e)


###################################################
### chunk number 10: howmanyprobesets
###################################################
dim(e)
dim(exprs(spikein95))


###################################################
### chunk number 11: 
###################################################
pData(eset)
pData(spikein95)


###################################################
### chunk number 12: columnAccess
###################################################
Index1 <- which(eset$population == "Pop1")
Index2 <- which(eset$population == "Pop2")


###################################################
### chunk number 13: gcrma
###################################################
library("gcrma")
eset2 <- gcrma(spikein95)


###################################################
### chunk number 14: threestep
###################################################
library("affyPLM")
eset3 <- threestep(spikein95, background=FALSE,
                   summary.method="tukey.biweight")


###################################################
### chunk number 15: rowM
###################################################
d <- rowMeans(e[, Index2]) - rowMeans(e[, Index1])


###################################################
### chunk number 16: rowMeans
###################################################
a <- rowMeans(e)


###################################################
### chunk number 17: figdvsa
###################################################
par(mfrow=c(1,2))
plot(a, d, pch=".", ylim=c(-2,2))
abline(h=0,col="blue")
plot(rank(a), d, pch=".", ylim=c(-2,2))
abline(h=0,col="pink")


###################################################
### chunk number 18: rowttests
###################################################
library("genefilter") 
tt <- rowttests(eset, "population")


###################################################
### chunk number 19: eBayesEx
###################################################
library("limma") 
design <- model.matrix(~eset$population) 
fit <- lmFit(eset, design)
ebayes <- eBayes(fit) 


###################################################
### chunk number 20: volcano
###################################################
lod <- -log10(tt$p.value)

plot(d, lod, cex=.25, 
  main="Volcano plot for $t$-test")
abline(h=2)


###################################################
### chunk number 21: volcanoeb
###################################################
plot(d, -log10(ebayes$p.value[,2]), xlim=c(-1,1), cex=.25, 
main="Volcano plot for t-test")
abline(h=2)


###################################################
### chunk number 22: volcano2
###################################################
plot(d,lod,cex=.25,xlim=c(-1,1),ylim=range(lod),main="Volcano plot")
o1 <- order(abs(d),decreasing=TRUE)[1:25]
points(d[o1],lod[o1],pch=18,col="blue")


###################################################
### chunk number 23: 
###################################################
table(tt$p.value<=0.01)


###################################################
### chunk number 24: toptable
###################################################
tab <- topTable(ebayes,coef=2,adjust.method="BH",n=10)
genenames <- as.character(tab$ID)


###################################################
### chunk number 25: annotate1
###################################################
library("annotate")


###################################################
### chunk number 26: annotate2
###################################################
annotation(eset)


###################################################
### chunk number 27: annotate3
###################################################
library("hgu95a")


###################################################
### chunk number 28: entrezGeneAndSybol
###################################################
ll <- getLL(genenames, "hgu95a")
sym <- getSYMBOL(genenames, "hgu95a")


###################################################
### chunk number 29: output
###################################################
tab <- data.frame(sym, signif(tab[,-1], 3))
htmlpage(ll, filename="GeneList1.html",
         title="HTML report",
         othernames=tab,
         table.head=c("Locus ID",colnames(tab)), table.center=TRUE)


###################################################
### chunk number 30: 
###################################################
browseURL("GeneList1.html")


###################################################
### chunk number 31: 
###################################################
library("KEGG")
library("GO")
library("annaffy") 
atab <- aafTableAnn(genenames, "hgu95av2",
    aaf.handler()) 
saveHTML(atab, file="GeneList2.html") 


###################################################
### chunk number 32: annaffy2
###################################################
aaf.handler()
atab <- aafTableAnn(genenames, "hgu95av2", aaf.handler()[c(2,5,8,12)])
saveHTML(atab, file="GeneList3.html")