################################################### ### 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]