################################################### ### chunk number 1: cwd ################################################### cwd=getwd() ################################################### ### chunk number 2: libraries ################################################### library(affy) library(estrogen) library(vsn) ################################################### ### chunk number 3: datadir ################################################### datadir = system.file("extdata", package="estrogen") datadir dir(datadir) setwd(datadir) ################################################### ### chunk number 4: loadpdata ################################################### pd = read.phenoData("estrogen.txt", header=TRUE, 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: get backup eval=FALSE ################################################### ## library(affy) ## library(estrogen) ## library(vsn) ## datadir = system.file("extdata", package="estrogen") ## setwd(datadir) ## load("workspace.RData") ################################################### ### chunk number 10: normmeth ################################################### normalize.methods(a) express.summary.stat.methods ################################################### ### chunk number 11: images1a eval=FALSE ################################################### ## image(a[, 1]) ################################################### ### chunk number 12: images1b ################################################### jpeg(file.path(cwd, "estrogen-image1.jpg"), quality=100) image(a[, 1]) dev.off() ################################################### ### chunk number 13: images2a eval=FALSE ################################################### ## badc = ReadAffy("bad.cel") ## image(badc) ################################################### ### chunk number 14: images2b ################################################### jpeg(file.path(cwd, "estrogen-image2.jpg"), quality=100) badc = ReadAffy("bad.cel") image(badc) dev.off() ################################################### ### chunk number 15: setwd ################################################### setwd(cwd) ################################################### ### chunk number 16: hist ################################################### hist(log2(intensity(a[, 4])), breaks=100, col="blue") ################################################### ### chunk number 17: boxplota eval=FALSE ################################################### ## boxplot(a,col="red") ## boxplot(data.frame(exprs(x)), col="blue") ################################################### ### chunk number 18: 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 19: classx ################################################### class(x) class(exprs(x)) ################################################### ### chunk number 20: scp-plots ################################################### 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 21: 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 22: heatmapa ################################################### rsd <- rowSds(exprs(x)) sel <- order(rsd, decreasing=TRUE)[1:50] ################################################### ### chunk number 23: heatmapb eval=FALSE ################################################### ## heatmap(exprs(x)[sel,], col=gentlecol(256)) ################################################### ### chunk number 24: heatmapc ################################################### jpeg(file.path(cwd, "estrogen-heatmap.jpg"), quality=100) heatmap(exprs(x)[sel,], col=gentlecol(256)) dev.off() ################################################### ### chunk number 25: deflm ################################################### lm.coef = function(y) lm(y ~ estrogen * time.h)$coefficients eff = esApply(x, 1, lm.coef) ################################################### ### chunk number 26: showlmres ################################################### dim(eff) rownames(eff) affyids <- colnames(eff) ################################################### ### chunk number 27: gn ################################################### library(hgu95av2) ls("package:hgu95av2") ################################################### ### chunk number 28: eff-show ################################################### par(mfrow=c(1,1)) hist(eff[2,], breaks=100, col="blue", main="estrogen main effect") setwd(cwd) ################################################### ### chunk number 29: 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)