## ----Setup, echo=FALSE, results='hide', message=FALSE------------------------- library("Biobase") library("annotate") library("GOstats") library("xtable") library("multtest") library("Rgraphviz") library("hgu95av2.db") library("GO.db") library("genefilter") library("ALL") ## ----Example, results='hide', echo=FALSE-------------------------------------- ## subset of interest: 37+42 samples data(ALL) eset <- ALL[, intersect( grep("^B", as.character(ALL$BT)), which(as.character(ALL$mol) %in% c("BCR/ABL", "NEG", "ALL1/AF4")) )] MB <- factor(eset$mol.bio) ## intensities above 100 in at least 7 of the samples f1 <- kOverA(7, log2(100)) f2 <- function(x) { gps <- split(2^x, MB) mns <- sapply(gps, mean) if (max(mns) - min(mns) < 100) FALSE else TRUE } ff <- filterfun(f1, f2) selected <- genefilter(eset, ff) sum(selected) esetSub <- eset[selected, ] ## ----mtFstats, results='hide', echo=FALSE------------------------------------- pvCutOff <- 0.05 Fstats <- mt.maxT(exprs(esetSub), as.numeric(MB) - 1, B = 1000, test = "f" ) eSet <- esetSub[Fstats$index[Fstats$adjp < pvCutOff], ] gN <- featureNames(eSet) lls <- unlist(mget(gN, hgu95av2ENTREZID, ifnotfound = NA)) eS <- eSet[!duplicated(lls), ] gN <- featureNames(eS) lls <- unlist(mget(gN, hgu95av2ENTREZID, ifnotfound = NA)) lls <- as.character(lls[!is.na(lls)]) syms <- unlist(mget(gN, hgu95av2SYMBOL, ifnotfound = NA)) ## ----inducedGO, echo=FALSE, results='hide'------------------------------------ gMF <- makeGOGraph(lls, "MF", chip = "hgu95av2") gBP <- makeGOGraph(lls, "BP", chip = "hgu95av2") gCC <- makeGOGraph(lls, "CC", chip = "hgu95av2") nMF <- list() lbs <- rep("", length(nodes(gMF))) names(lbs) <- nodes(gMF) nMF$label <- lbs nBP <- list() lbs <- rep("", length(nodes(gBP))) names(lbs) <- nodes(gBP) nBP$label <- lbs nCC <- list() lbs <- rep("", length(nodes(gCC))) names(lbs) <- nodes(gCC) nCC$label <- lbs if (require("Rgraphviz")) { gMFlo <- agopen(gMF, "gMF") gBPlo <- agopen(gBP, "gBP") gCClo <- agopen(gCC, "gCC") } ## ----GOMF, fig=TRUE, echo=FALSE, warning=FALSE, fig.cap = " The induced MF GO graph for the selected genes."---- if (require("Rgraphviz")) { plot(gMFlo, nodeAttrs = nMF, main = "Molecular Function") } ## ----GOBP, fig=TRUE, echo=FALSE, warning=FALSE, fig.cap = " The induced BP GO graph for the selected genes."---- if (require("Rgraphviz")) { plot(gBPlo, nodeAttrs = nBP, main = "Biological Process") } ## ----GOCC, fig=TRUE, echo=FALSE, warning=FALSE, fig.cap = "The induced CC GO graph for the selected genes."---- if (require("Rgraphviz")) { plot(gCClo, nodeAttrs = nCC, main = "Cellular Component") } ## ----findhigestmeans---------------------------------------------------------- mns <- apply(exprs(eS), 1, function(x) sapply(split(x, MB), mean)) whismax <- apply(mns, 2, function(x) match(max(x), x)) maxNames <- names(mns[, 1])[whismax] names(maxNames) <- names(whismax) table(maxNames) ## ----genesToGO---------------------------------------------------------------- CCnodes <- nodes(gCC) nodes2affy <- mget(CCnodes, hgu95av2GO2ALLPROBES, ifnotfound = NA) cts <- sapply(nodes2affy, function(x) { wh <- names(whismax) %in% x table(maxNames[wh]) }) all1cts <- sapply(cts, function(x) x["ALL1/AF4"]) all1cts <- ifelse(is.na(all1cts), 0, all1cts) ## put nice names on names(all1cts) <- names(cts) bcrcts <- sapply(cts, function(x) x["BCR/ABL"]) bcrcts <- ifelse(is.na(bcrcts), 0, bcrcts) names(bcrcts) <- names(cts) negcts <- sapply(cts, function(x) x["NEG"]) negcts <- ifelse(is.na(negcts), 0, negcts) names(negcts) <- names(cts) ctmat <- cbind(all1cts, bcrcts, negcts) ## ----layoutandrender, eval=FALSE, warning=FALSE, results="hide"--------------- # if (require(Rgraphviz)) { # opar <- par(xpd = NA) # plotPieChart <- function(curPlot, counts, main) { # renderNode <- function(x) { # force(x) # y <- x * 100 + 1 # function(node, ur, attrs = list(), radConv = 1) { # nodeCenter <- getNodeCenter(node) # pieGlyph(y, # xpos = getX(nodeCenter), # ypos = getY(nodeCenter), # radius = getNodeRW(node), # col = c("blue", "green", "red") # ) # } # } # drawing <- vector(mode = "list", length = nrow(counts)) # for (i in 1:length(drawing)) { # drawing[[i]] <- renderNode(counts[i, ]) # } # if (missing(main)) { # main <- "Example Pie Chart Plot" # } # # plot(curPlot, drawNode = drawing, main = main) # # legend( # x = "bottomleft", legend = c("ALL1/AF4", "BCR/ABL", "NEG"), # fill = c("blue", "green", "red") # ) # } # plotPieChart(gCClo, ctmat) # par(opar) # } ## ----imageMap, message=FALSE, results='hide'---------------------------------- getNodeBB <- function(ingraph) { xypos <- getNodeXY(ingraph) hts <- getNodeHeight(ingraph) wdsR <- getNodeRW(ingraph) wdsL <- getNodeLW(ingraph) llx <- xypos$x - wdsL lly <- xypos$y - (hts / 2) urx <- xypos$x + wdsR ury <- xypos$y + (hts / 2) cbind(llx, lly, urx, ury) } CCbb <- getNodeBB(gCClo) ## now we need to adjust the graphBB to the plot region bbG <- boundBox(gCClo) llx <- getX(botLeft(bbG)) lly <- getY(botLeft(bbG)) urx <- getX(upRight(bbG)) ury <- getY(upRight(bbG)) ## FIXME: work in progress here - we need to do a bit of mapping ## to get the right user coords on the jpg file... if (interactive()) { jpeg("mygraph.jpg", quality = 100, width = urx - llx, height = ury - lly) opar <- par(plt = c(0, 1, 0, 1), xpd = NA) plotPieChart(gCClo, ctmat) par(opar) dev.off() } ## FIXME (wh 20.1.2005): better to use imageMap function for Ragraph objects ## in the package Rgraphviz instead. if (require("geneplotter")) { con <- openHtmlPage("example", "An Image Map") imageMap(CCbb, con, tags = list( TITLE = getNodeNames(gCClo), HREF = rep("", length(nodes(gCC))) ), imgname = "mygraph.jpg" ) closeHtmlPage(con) } ## ----multip------------------------------------------------------------------- igenes <- unlist(mget(gN[1:10], hgu95av2ENTREZID)) igenes <- igenes[!is.na(igenes)] igenes <- as.character(igenes) psets <- mget(igenes, revmap(hgu95av2ENTREZID)) psets <- sapply(psets, function(x) x[1]) GOs <- mget(psets, hgu95av2GO, ifnotfound = NA) gBP <- sapply(GOs, getOntology, "BP") ## how many terms sum(sapply(gBP, length)) ## drop those with no BP annotations hasBP <- sapply(gBP, function(x) length(x) > 0 && !is.na(x[1])) gBP <- gBP[hasBP] ggs <- lapply(igenes[hasBP], makeGOGraph, "BP", chip = "hgu95av2.db") ## you could also do ggx = lapply(gBP, GOGraph, GOBPPARENTS) ## and should get the same thing simatM1 <- matrix(1, nr = sum(hasBP), nc = sum(hasBP)) for (i in 1:sum(hasBP)) { for (j in 1:sum(hasBP)) { if (i == j) { next } else { simatM1[i, j] <- simUI(ggs[[i]], ggs[[j]]) } } } library("RBGL") simatM2 <- matrix(1, nr = sum(hasBP), nc = sum(hasBP)) for (i in 1:sum(hasBP)) { for (j in 1:sum(hasBP)) { if (i == j) { next } else { simatM2[i, j] <- simLP(ggs[[i]], ggs[[j]]) } } } ## ----leavesEx----------------------------------------------------------------- gCC.leaves <- leaves(gCC, "in") length(gCC.leaves) gCC.leaves[1:5]