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


###################################################
### chunk number 2: 
###################################################
library("affy")


###################################################
### chunk number 3: readaffy eval=FALSE
###################################################
## library("affy")
## Data <- ReadAffy()


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


###################################################
### chunk number 5: 
###################################################
slotNames(spikein95)


###################################################
### chunk number 6: 
###################################################
pd <- data.frame(population=c(1,1,1,2,2,2),
replicate=c(1,2,3,1,2,3))
rownames(pd) <- sampleNames(spikein95)
vl <- list(population="1 is control, 2 is treatment",
           replicate="arbitrary numbering") 
phenoData(spikein95) <- new("phenoData",pData=pd,varLabels=vl)


###################################################
### chunk number 7: 
###################################################
library("ALLMLL")
data(MLL.B)
Data <- MLL.B[,c(2,1,3:5,14,6,13)]
sampleNames(Data) <- letters[1:8]
rm(MLL.B)
gc()


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


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


###################################################
### chunk number 10: 
###################################################
pData(eset)


###################################################
### chunk number 11: 
###################################################
Index1 <- which(eset$population==1)
Index2 <- which(eset$population==2)


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


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


###################################################
### chunk number 14: sumAbs
###################################################
sum(abs(d)>1)

###################################################
### chunk number 15: maplot
###################################################
par(cex.main=0.8)
plot(a,d,ylim=c(-1,1),main="A) MA-plot",pch=".")

###################################################
### chunk number 16: rowttests
###################################################
library("genefilter")
tt <- rowttests(e, factor(eset$population))


###################################################
### chunk number 17: volcano
###################################################
lod <- -log10(tt$p.value)
plot(d,lod,cex=.25,main="B) Volcano plot for $t$-test")
abline(h=2)


###################################################
### chunk number 18: volcano2
###################################################
lod <- -log10(tt$p.value)
plot(d,lod,cex=.25,main="B) Volcano plot for $t$-test")
abline(h=2)


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


###################################################
### chunk number 20: volcano3
###################################################
lod <- -log10(ebayes$p.value[,2])
mtstat<- ebayes$t[,2]
o1 <- order(abs(d),decreasing=TRUE)[1:25]
o2 <- order(abs(mtstat),decreasing=TRUE)[1:25]
o <- union(o1,o2)
plot(d[-o],lod[-o],cex=.25,xlim=c(-1,1),ylim=c(0,4),
main="D) Volcano plot for moderated $t$-test")
points(d[o1],lod[o1],pch=18,col="blue")
points(d[o2],lod[o2],pch=1,col="red")


###################################################
### chunk number 21: 
###################################################
sum(tt$p.value<=0.01)


###################################################
### chunk number 22: 
###################################################
data(spikein95)
spikedin <- colnames(pData(spikein95))
spikedIndex <- match(spikedin,geneNames(eset))


###################################################
### chunk number 23: 
###################################################
d.rank <- sort(rank(-abs(d))[spikedIndex])
t.rank <- sort(rank(-abs(tt$statistic))[spikedIndex])
mt.rank <- sort(rank(-abs(mtstat))[spikedIndex])
ranks <- cbind(mt.rank,d.rank,t.rank)
rownames(ranks) <- NULL
ranks


###################################################
### chunk number 24: 
###################################################
tab <- topTable(ebayes,coef=2,adjust="fdr",n=10)


###################################################
### chunk number 25: 
###################################################
tab[1:5,]


###################################################
### chunk number 26: 
###################################################
genenames <- as.character(tab$ID)


###################################################
### chunk number 27: 
###################################################
library("hgu95av2")
hgu95av2()
gn <- mget(genenames, hgu95av2GENENAME)
sapply(gn, function(x) if(!is.na(x)) strwrap(x) else x)


###################################################
### chunk number 28: 
###################################################
ls(hgu95av2GENENAME)[1:10]


###################################################
### chunk number 29: 
###################################################
library("GO")
library("annotate")
gos <- mget(genenames, hgu95av2GO)
gos[[1]][1:2]
goids <- sapply(gos, names)
goids[[1]]
goterms <- sapply(goids, function(x) mget(as.character(x), GOTERM))
sapply(goterms[[1]], Term)


###################################################
### chunk number 30: loadlibxml
###################################################
library("XML")
absts<-pm.getabst(genenames,"hgu95av2")


###################################################
### chunk number 31: absts12
###################################################
absts[[1]][[4]]


###################################################
### chunk number 32: pmtit2
###################################################
## wh 13.01.05: commented out pm.titles since it behaves weird. 
## Revert back to it when it is fixed?
##titl <- pm.titles(absts[2])
titl <- sapply(absts[[2]], articleTitle)
strwrap(titl, simplify=FALSE)


###################################################
### chunk number 33: 
###################################################
pro.res <- sapply(absts, function(x) pm.abstGrep("[Pp]rotein", x))
pro.res[[2]]


###################################################
### chunk number 34: 
###################################################
pmAbst2HTML(absts[[2]],filename="pm.html")


###################################################
### chunk number 35: 
###################################################
ll <- getLL(genenames,"hgu95av2")
ug <- unlist(lookUp(genenames, "hgu95av2", "UNIGENE"))
sym <- getSYMBOL(genenames,"hgu95av2")


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


###################################################
### chunk number 37: 
###################################################
library("annaffy")
atab <- aafTableAnn( genenames, "hgu95av2", aaf.handler() )
saveHTML(atab, file="report2.html")


###################################################
### chunk number 38: 
###################################################
library("GOstats")


###################################################
### chunk number 39: 
###################################################
index <- sample(1:12000, 100)
gn <- ls(hgu95av2GENENAME)[index]
lls <- unique(getLL(gn, "hgu95av2"))
hyp <- GOHyperG(lls, "hgu95av2", what = "MF")
index2 <- hyp$pvalues < 0.05
wh <- mget(names(hyp$pvalues)[index2], GOTERM)
terms <- sapply(wh, Term)
strwrap(terms[1:10], simplify = FALSE)
hyp$pvalues[index2][1:10]