Skip to content.

bioconductor.org

Bioconductor is an open source and open development software project
for the analysis and comprehension of genomic data.

Sections

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()

News
2009-10-26

BioC 2.5, consisting of 352 packages and designed to work with R 2.10.z, was released today.

2009-01-07

R, the open source platform used by Bioconductor, featured in a series of articles in the New York Times.