LabAnnotate.R
################################################### ### chunk number 1: graphisoff ################################################### graphics.off()
################################################### ### chunk number 2: loadlibs ################################################### library("Biobase") library("ALL") library("hgu95av2") library("annotate") data("ALL") ALL
################################################### ### chunk number 3: subset ################################################### ALLs1 = ALL[, ALL$mol.biol %in% c("ALL1/AF4", "BCR/ABL")]
################################################### ### chunk number 4: filter1 ################################################### library("genefilter") f1 <- function(x)(IQR(x)>0.8) ff <- filterfun(f1) fil <- genefilter(ALLs1, ff)
################################################### ### chunk number 5: filter2 ################################################### table(fil) ALLs2 <- ALLs1[fil,]
################################################### ### chunk number 6: rt ################################################### rt = rowttests(ALLs2, "mol.biol") names(rt) str(rt)
################################################### ### chunk number 7: figrtsol1 ################################################### hist(rt$statistic, breaks=100, col="skyblue")
################################################### ### chunk number 8: figrtsol2 ################################################### hist(rt$p.value, breaks=100, col="mistyrose")
################################################### ### chunk number 9: ALLs3 ################################################### sel <- order(rt$p.value)[1:100] ALLs3 = ALLs2[sel,]
################################################### ### chunk number 10: tabletable ################################################### EGs3 = unlist(mget(geneNames(ALLs3), hgu95av2LOCUSID)) table(table(EGs3))
################################################### ### chunk number 11: figheatmap ################################################### hm = heatmap(exprs(ALLs3), ColSide = ifelse(ALLs3$mol=="ALL1/AF4", "skyblue", "salmon"), col = topo.colors(255), keep.dendro=TRUE)
################################################### ### chunk number 12: figdend ################################################### plot(hm$Rowv, leaflab="none")
################################################### ### chunk number 13: cutdend ################################################### cdend = cut(hm$Rowv, h=25) dendlabs = lapply(cdend$lower, labels) listLen(dendlabs) dendlabs[[1]]
################################################### ### chunk number 14: annaffy ################################################### library("annaffy") anncols <- aaf.handler(chip="hgu95av2")[c(1:3, 8:9, 11:13)] anntable <- aafTableAnn(dendlabs[[1]], "hgu95av2", anncols) saveHTML(anntable, "example1.html", title="Cluster Number 1") browseURL("example1.html")
################################################### ### chunk number 15: multiplicities ################################################### lls <- unlist(as.list(hgu95av2LOCUSID)) lls[1:12] table(table(lls))
################################################### ### chunk number 16: somebutnotall ################################################### probeSetsPerGene = split(names(lls), lls) myFun = function(x) { all(fil[x]) || !any(fil[x])} isConsistent = sapply(probeSetsPerGene, myFun) table(isConsistent)
################################################### ### chunk number 17: extraparanoia ################################################### stopifnot(!isConsistent["7013"])
################################################### ### chunk number 18: 7013a ################################################### isConsistent["7013"] probeSetsPerGene["7013"] fil[probeSetsPerGene$"7013"]
################################################### ### chunk number 19: 7013b ################################################### mget(probeSetsPerGene$"7013", hgu95av2SYMBOL)
################################################### ### chunk number 20: figduplo1 ################################################### j <- probeSetsPerGene$"7013" plot(t(exprs(ALL)[j[1:2], ]), pch=16)
################################################### ### chunk number 21: figduplo2 ################################################### matplot(t(exprs(ALL)[j, order(ALL$mol.biol)]), type="l", lty=1, col=1:7, ylab="Intensity", xlab="Sample Index", main = "Intensities for probesets for EntrezGene 7013")
################################################### ### chunk number 22: pchr ################################################### chrList = as.list(hgu95av2CHR) table(listLen(chrList)) pchr = unlist(chrList) table(pchr) pchr[1:10]
################################################### ### chunk number 23: chr2 ################################################### stopifnot(all(listLen(chrList)==1))
################################################### ### chunk number 24: lls ################################################### lls[1:10]
################################################### ### chunk number 25: chr ################################################### pchr <- pchr[order(names(pchr))] lls <- lls[order(names(lls))] stopifnot(identical(names(lls), names(pchr))) chr <- cbind("EntrezGene ID" = lls, "Chromosome"=pchr) chr <- chr[!duplicated(chr[,1]), ] chr[1:5, ]
################################################### ### chunk number 26: selGenes ################################################### selGenes <- as.character(lls[geneNames(ALLs3)])
################################################### ### chunk number 27: selGenes ################################################### isdiff <- chr[, "EntrezGene ID"] %in% selGenes tab <- table(isdiff, chr[, "Chromosome"]) fisher.test(tab, simulate.p.value=TRUE) chisq.test(tab)
###################################################
### chunk number 28: figtab
###################################################
plot(tab[1,], tab[2,], pch=16, xlab=not diff. expressed
, ylab=diff. expressed
)
################################################### ### chunk number 29: loadGOstats ################################################### library("GOstats")
################################################### ### chunk number 30: mfhyper ################################################### mfhyper = GOHyperG(selGenes)
################################################### ### chunk number 31: figmfhyper ################################################### hist(mfhyper$pvalues, breaks=50, col="mistyrose")
################################################### ### chunk number 32: loadRgraphviz ################################################### library("Rgraphviz")
################################################### ### chunk number 33: setuplotting ################################################### gGO <- makeGOGraph(selGenes, "MF", TRUE)
p <- mfhyper$pvalues[nodes(gGO)]
nodeColor <- ifelse(p < 0.1, "tomato", "lightblue") nodeLabel <- character(numNodes(gGO)) names(nodeColor) <- names(nodeLabel) <- nodes(gGO)
nA <- list(label = nodeLabel, fillcolor = nodeColor)
################################################### ### chunk number 34: figgGO ################################################### plot(gGO, nodeAttrs=nA)
################################################### ### chunk number 35: ################################################### sessionInfo()