### R code from vignette source 'vignettes/estrogen/inst/doc/estrogen.Rnw' ################################################### ### code chunk number 1: cwd ################################################### cwd=getwd() ################################################### ### code chunk number 2: libraries ################################################### library(affy) library(estrogen) library(vsn) library(genefilter) ################################################### ### code chunk number 3: datadir ################################################### datadir = system.file("extdata", package="estrogen") datadir dir(datadir) setwd(datadir) ################################################### ### code chunk number 4: loadpdata ################################################### pd = read.AnnotatedDataFrame("estrogen.txt", header=TRUE, sep="", row.names=1) pData(pd) ################################################### ### code chunk number 5: a.read ################################################### a = ReadAffy(filenames = rownames(pData(pd)), phenoData = pd, verbose=TRUE) ################################################### ### code chunk number 6: a.show ################################################### a ################################################### ### code 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") ################################################### ### code chunk number 8: x.show ################################################### x ################################################### ### code chunk number 9: normmeth ################################################### normalize.methods(a) express.summary.stat.methods() ################################################### ### code chunk number 10: images1a (eval = FALSE) ################################################### ## image(a[, 1]) ################################################### ### code chunk number 11: images1b ################################################### jpeg(file.path(cwd, "estrogen-image1.jpg"), quality=100) image(a[, 1]) dev.off() ################################################### ### code chunk number 12: images2a (eval = FALSE) ################################################### ## badc = ReadAffy("bad.cel") ## image(badc) ################################################### ### code chunk number 13: images2b (eval = FALSE) ################################################### ## jpeg(file.path(cwd, "estrogen-image2.jpg"), quality=100) ## badc = ReadAffy("bad.cel") ## image(badc) ## dev.off() ################################################### ### code chunk number 14: setwd ################################################### setwd(cwd) ################################################### ### code chunk number 15: hist ################################################### hist(log2(intensity(a[, 4])), breaks=100, col="blue") ################################################### ### code chunk number 16: boxplota (eval = FALSE) ################################################### ## boxplot(a,col="red") ## boxplot(data.frame(exprs(x)), col="blue") ################################################### ### code 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)) ################################################### ### code chunk number 18: classx ################################################### class(x) class(exprs(x)) ################################################### ### code 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() ################################################### ### code 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") ################################################### ### code chunk number 21: heatmapa ################################################### rsd <- rowSds(exprs(x)) sel <- order(rsd, decreasing=TRUE)[1:50] ################################################### ### code chunk number 22: heatmapb (eval = FALSE) ################################################### ## heatmap(exprs(x)[sel,], col=gentlecol(256)) ################################################### ### code chunk number 23: heatmapc (eval = FALSE) ################################################### ## jpeg(file.path(cwd, "estrogen-heatmap.jpg"), quality=100) ## heatmap(exprs(x)[sel,], col=gentlecol(256)) ## dev.off() ################################################### ### code chunk number 24: deflm ################################################### lm.coef = function(y) lm(y ~ estrogen * time.h)$coefficients eff = esApply(x, 1, lm.coef) ################################################### ### code chunk number 25: showlmres ################################################### dim(eff) rownames(eff) affyids <- colnames(eff) ################################################### ### code chunk number 26: gn ################################################### library(hgu95av2.db) ls("package:hgu95av2.db") ################################################### ### code chunk number 27: eff-show ################################################### par(mfrow=c(1,1)) hist(eff[2,], breaks=100, col="blue", main="estrogen main effect") setwd(cwd) ################################################### ### code 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)