###################################################
### chunk number 1: cwd
###################################################
cwd=getwd()


###################################################
### chunk number 2: libraries
###################################################
library(affy)
library(estrogen)
library(vsn)
library(genefilter)


###################################################
### chunk number 3: datadir
###################################################
datadir = system.file("extdata", package="estrogen")
datadir
dir(datadir)
setwd(datadir)


###################################################
### chunk number 4: loadpdata
###################################################
pd = read.AnnotatedDataFrame("estrogen.txt", header=TRUE, sep="", row.names=1)
pData(pd)


###################################################
### chunk number 5: a.read
###################################################
a = ReadAffy(filenames = rownames(pData(pd)), 
             phenoData = pd,
             verbose=TRUE)


###################################################
### chunk number 6: a.show
###################################################
a


###################################################
### chunk number 7: x.calc
###################################################
x <- expresso(a,
	     bg.correct       = FALSE,   ## bg correction is done by vsn
             normalize.method = "vsn",
             normalize.param  = list(subsample=1000), 
             pmcorrect.method = "pmonly", 
             summary.method   = "medianpolish")


###################################################
### chunk number 8: x.show
###################################################
x


###################################################
### chunk number 9: normmeth
###################################################
normalize.methods(a)
express.summary.stat.methods


###################################################
### chunk number 10: images1a eval=FALSE
###################################################
## image(a[, 1])


###################################################
### chunk number 11: images1b
###################################################
jpeg(file.path(cwd, "estrogen-image1.jpg"), quality=100)
image(a[, 1])
dev.off()


###################################################
### chunk number 12: images2a eval=FALSE
###################################################
## badc = ReadAffy("bad.cel")
## image(badc)


###################################################
### chunk number 13: images2b eval=FALSE
###################################################
## jpeg(file.path(cwd, "estrogen-image2.jpg"), quality=100)
## badc = ReadAffy("bad.cel")
## image(badc)
## dev.off()


###################################################
### chunk number 14: setwd
###################################################
setwd(cwd)


###################################################
### chunk number 15: hist
###################################################
hist(log2(intensity(a[, 4])), breaks=100, col="blue")


###################################################
### chunk number 16: boxplota eval=FALSE
###################################################
## boxplot(a,col="red")
## boxplot(data.frame(exprs(x)), col="blue")


###################################################
### chunk number 17: boxplotb
###################################################
par(mfrow=c(1,2))
colnames(exprs(a)) = colnames(exprs(x)) = paste(1:8)
boxplot(a,col="red")
boxplot(data.frame(exprs(x)), col="blue")
par(mfrow=c(1,1))


###################################################
### chunk number 18: classx
###################################################
class(x)
class(exprs(x))


###################################################
### chunk number 19: scp-plots eval=FALSE
###################################################
## jpeg(file.path(cwd, "estrogen-scp1.jpg"), quality=100)
## plot(exprs(a)[, 1:2], log="xy", pch=".", main="all")
## dev.off()
## 
## jpeg(file.path(cwd, "estrogen-scp2.jpg"), quality=100)
## plot(pm(a)[, 1:2], log="xy", pch=".", main="pm")
## dev.off()
## 
## jpeg(file.path(cwd, "estrogen-scp3.jpg"), quality=100)
## plot(mm(a)[, 1:2], log="xy", pch=".", main="mm")
## dev.off()
## 
## jpeg(file.path(cwd, "estrogen-scp4.jpg"), quality=100)
## plot(exprs(x)[, 1:2], pch=".", main="x")
## dev.off()


###################################################
### chunk number 20: scp eval=FALSE
###################################################
## plot(exprs(a)[,1:2], log="xy", pch=".", main="all")
## plot(pm(a)[, 1:2], log="xy", pch=".", main="pm")
## plot(mm(a)[, 1:2], log="xy", pch=".", main="mm")
## plot(exprs(x)[, 1:2], pch=".", main="x")


###################################################
### chunk number 21: heatmapa
###################################################
rsd <- rowSds(exprs(x))
sel <- order(rsd, decreasing=TRUE)[1:50] 


###################################################
### chunk number 22: heatmapb eval=FALSE
###################################################
## heatmap(exprs(x)[sel,], col=gentlecol(256))


###################################################
### chunk number 23: heatmapc eval=FALSE
###################################################
## jpeg(file.path(cwd, "estrogen-heatmap.jpg"), quality=100)
## heatmap(exprs(x)[sel,], col=gentlecol(256))
## dev.off()


###################################################
### chunk number 24: deflm
###################################################
lm.coef = function(y) 
  lm(y ~ estrogen * time.h)$coefficients
eff = esApply(x, 1, lm.coef)


###################################################
### chunk number 25: showlmres
###################################################
dim(eff)
rownames(eff)
affyids <- colnames(eff)


###################################################
### chunk number 26: gn
###################################################
library(hgu95av2) 
ls("package:hgu95av2")


###################################################
### chunk number 27: eff-show
###################################################
par(mfrow=c(1,1))
hist(eff[2,], breaks=100, col="blue", main="estrogen main effect")
setwd(cwd)


###################################################
### chunk number 28: eff
###################################################
lowest <- sort(eff[2,], decreasing=FALSE)[1:3]
mget(names(lowest), hgu95av2GENENAME)

highest <- sort(eff[2,], decreasing=TRUE)[1:3]
mget(names(highest), hgu95av2GENENAME)

hist(eff[4,], breaks=100, col="blue", main="estrogen:time interaction")
highia <- sort(eff[4,], decreasing=TRUE)[1:3]
mget(names(highia), hgu95av2GENENAME)