################################################### ### 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) ## There's a problem with the name of the annotation data ## so we fix it here. Normally, this should not be needed. spikein95@cdfName <- "hgu95av2"

################################################### ### 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 <- list(population="Pop1 is control, Pop2 is treatment", replicate="1, 2, 3, arbitrary numbering") phenoData(spikein95) <- new("phenoData", pData=pd, varLabels=vl)

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

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

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

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

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

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

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

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

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

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

################################################### ### chunk number 16: 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 17: rowttests ################################################### library("genefilter") tt <- rowttests(eset, "population")

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

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

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

################################################### ### chunk number 20: 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 21: 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 22: ################################################### table(tt$p.value<=0.01)

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

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

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

################################################### ### chunk number 26: annotate3 ################################################### library("hgu95av2")

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

################################################### ### chunk number 28: 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 29: ################################################### browseURL("GeneList1.html")

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

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