################################################### ### chunk number 1: ################################################### #line 41 "vignettes/topGO/inst/doc/topGO.Rnw" options(width = 95) ## we better keep the data in data frames as strings options(stringsAsFactors = FALSE) ################################################### ### chunk number 2: ################################################### #line 68 "vignettes/topGO/inst/doc/topGO.Rnw" x <- topGO:::.algoComp[, -8] x <- x[, colSums(x) > 0] yesPic <- "\\includegraphics[width=3mm]{green_ckmark.png}" noPic <- "\\includegraphics[width=3mm]{red_ckmark.png}" x[x == 1] <- yesPic x[x == "0"] <- noPic ################################################### ### chunk number 3: ################################################### #line 79 "vignettes/topGO/inst/doc/topGO.Rnw" if(require(xtable)) print(xtable(x, align = c("l", rep("c", ncol(x)))), type = "latex", sanitize.text.function = function(x) return(x), floating = FALSE) ################################################### ### chunk number 4: eval=FALSE ################################################### ## #line 113 "vignettes/topGO/inst/doc/topGO.Rnw" ## source("http://bioconductor.org/biocLite.R") ## biocLite() ################################################### ### chunk number 5: eval=FALSE ################################################### ## #line 120 "vignettes/topGO/inst/doc/topGO.Rnw" ## source("http://bioconductor.org/biocLite.R") ## biocLite("topGO") ################################################### ### chunk number 6: ################################################### #line 164 "vignettes/topGO/inst/doc/topGO.Rnw" library(topGO) library(ALL) data(ALL) data(geneList) ################################################### ### chunk number 7: ################################################### #line 178 "vignettes/topGO/inst/doc/topGO.Rnw" affyLib <- paste(annotation(ALL), "db", sep = ".") library(package = affyLib, character.only = TRUE) ################################################### ### chunk number 8: ################################################### #line 190 "vignettes/topGO/inst/doc/topGO.Rnw" sum(topDiffGenes(geneList)) ################################################### ### chunk number 9: ################################################### #line 198 "vignettes/topGO/inst/doc/topGO.Rnw" sampleGOdata <- new("topGOdata", description = "Simple session", ontology = "BP", allGenes = geneList, geneSel = topDiffGenes, nodeSize = 10, annot = annFUN.db, affyLib = affyLib) ################################################### ### chunk number 10: ################################################### #line 216 "vignettes/topGO/inst/doc/topGO.Rnw" sampleGOdata ################################################### ### chunk number 11: ################################################### #line 238 "vignettes/topGO/inst/doc/topGO.Rnw" resultFisher <- runTest(sampleGOdata, algorithm = "classic", statistic = "fisher") ################################################### ### chunk number 12: ################################################### #line 244 "vignettes/topGO/inst/doc/topGO.Rnw" resultFisher ################################################### ### chunk number 13: ################################################### #line 251 "vignettes/topGO/inst/doc/topGO.Rnw" resultKS <- runTest(sampleGOdata, algorithm = "classic", statistic = "ks") resultKS.elim <- runTest(sampleGOdata, algorithm = "elim", statistic = "ks") ################################################### ### chunk number 14: ################################################### #line 271 "vignettes/topGO/inst/doc/topGO.Rnw" allRes <- GenTable(sampleGOdata, classicFisher = resultFisher, classicKS = resultKS, elimKS = resultKS.elim, orderBy = "elimKS", ranksOf = "classicFisher", topNodes = 10) ################################################### ### chunk number 15: ################################################### #line 284 "vignettes/topGO/inst/doc/topGO.Rnw" if(require(xtable)) print(xtable(apply(allRes, 2, as.character)), floating = FALSE) ################################################### ### chunk number 16: ################################################### #line 301 "vignettes/topGO/inst/doc/topGO.Rnw" colMap <- function(x) { .col <- rep(rev(heat.colors(length(unique(x)))), time = table(x)) return(.col[match(1:length(x), order(x))]) } ################################################### ### chunk number 17: eval=FALSE ################################################### ## #line 308 "vignettes/topGO/inst/doc/topGO.Rnw" ## pValue.classic <- score(resultKS) ## pValue.elim <- score(resultKS.elim)[names(pValue.classic)] ## ## gstat <- termStat(sampleGOdata, names(pValue.classic)) ## gSize <- gstat$Annotated / max(gstat$Annotated) * 4 ## gCol <- colMap(gstat$Significant) ## ## plot(pValue.classic, pValue.elim, xlab = "p-value classic", ylab = "p-value elim", ## pch = 19, cex = gSize, col = gCol) ################################################### ### chunk number 18: ################################################### #line 324 "vignettes/topGO/inst/doc/topGO.Rnw" pValue.classic <- score(resultKS) pValue.elim <- score(resultKS.elim)[names(pValue.classic)] gstat <- termStat(sampleGOdata, names(pValue.classic)) gSize <- gstat$Annotated / max(gstat$Annotated) * 4 gCol <- colMap(gstat$Significant) par(mfcol = c(1, 2), cex = 1) plot(pValue.classic, pValue.elim, xlab = "p-value classic", ylab = "p-value elim", pch = 19, cex = gSize, col = gCol) plot(pValue.classic, pValue.elim, log = "xy", xlab = "log(p-value) classic", ylab = "log(p-value) elim", pch = 19, cex = gSize, col = gCol) ################################################### ### chunk number 19: ################################################### #line 355 "vignettes/topGO/inst/doc/topGO.Rnw" sel.go <- names(pValue.classic)[pValue.elim < pValue.classic] cbind(termStat(sampleGOdata, sel.go), elim = pValue.elim[sel.go], classic = pValue.classic[sel.go]) ################################################### ### chunk number 20: eval=FALSE ################################################### ## #line 375 "vignettes/topGO/inst/doc/topGO.Rnw" ## showSigOfNodes(sampleGOdata, score(resultKS.elim), firstSigNodes = 5, useInfo = 'all') ################################################### ### chunk number 21: ################################################### #line 379 "vignettes/topGO/inst/doc/topGO.Rnw" printGraph(sampleGOdata, resultKS.elim, firstSigNodes = 5, fn.prefix = "tGO", useInfo = "all", pdfSW = TRUE) ################################################### ### chunk number 22: ################################################### #line 412 "vignettes/topGO/inst/doc/topGO.Rnw" library(topGO) library(ALL) data(ALL) ################################################### ### chunk number 23: ################################################### #line 423 "vignettes/topGO/inst/doc/topGO.Rnw" BPterms <- ls(GOBPTerm) head(BPterms) ################################################### ### chunk number 24: ################################################### #line 438 "vignettes/topGO/inst/doc/topGO.Rnw" library(genefilter) selProbes <- genefilter(ALL, filterfun(pOverA(0.20, log2(100)), function(x) (IQR(x) > 0.25))) eset <- ALL[selProbes, ] ################################################### ### chunk number 25: ################################################### #line 540 "vignettes/topGO/inst/doc/topGO.Rnw" geneID2GO <- readMappings(file = system.file("examples/geneid2go.map", package = "topGO")) str(head(geneID2GO)) ################################################### ### chunk number 26: ################################################### #line 564 "vignettes/topGO/inst/doc/topGO.Rnw" GO2geneID <- inverseList(geneID2GO) str(head(GO2geneID)) ################################################### ### chunk number 27: ################################################### #line 580 "vignettes/topGO/inst/doc/topGO.Rnw" geneNames <- names(geneID2GO) head(geneNames) ################################################### ### chunk number 28: ################################################### #line 588 "vignettes/topGO/inst/doc/topGO.Rnw" myInterestingGenes <- sample(geneNames, length(geneNames) / 10) geneList <- factor(as.integer(geneNames %in% myInterestingGenes)) names(geneList) <- geneNames str(geneList) ################################################### ### chunk number 29: ################################################### #line 606 "vignettes/topGO/inst/doc/topGO.Rnw" GOdata <- new("topGOdata", ontology = "MF", allGenes = geneList, annot = annFUN.gene2GO, gene2GO = geneID2GO) ################################################### ### chunk number 30: ################################################### #line 622 "vignettes/topGO/inst/doc/topGO.Rnw" GOdata ################################################### ### chunk number 31: ################################################### #line 646 "vignettes/topGO/inst/doc/topGO.Rnw" y <- as.integer(sapply(eset$BT, function(x) return(substr(x, 1, 1) == 'T'))) table(y) ################################################### ### chunk number 32: ################################################### #line 657 "vignettes/topGO/inst/doc/topGO.Rnw" geneList <- getPvalues(exprs(eset), classlabel = y, alternative = "greater") ################################################### ### chunk number 33: ################################################### #line 671 "vignettes/topGO/inst/doc/topGO.Rnw" topDiffGenes <- function(allScore) { return(allScore < 0.01) } x <- topDiffGenes(geneList) sum(x) ## the number of selected genes ################################################### ### chunk number 34: ################################################### #line 682 "vignettes/topGO/inst/doc/topGO.Rnw" GOdata <- new("topGOdata", description = "GO analysis of ALL data; B-cell vs T-cell", ontology = "BP", allGenes = geneList, geneSel = topDiffGenes, annot = annFUN.db, nodeSize = 5, affyLib = affyLib) ################################################### ### chunk number 35: ################################################### #line 722 "vignettes/topGO/inst/doc/topGO.Rnw" allProb <- featureNames(ALL) groupProb <- integer(length(allProb)) + 1 groupProb[allProb %in% genes(GOdata)] <- 0 groupProb[!selProbes] <- 2 groupProb <- factor(groupProb, labels = c("Used", "Not annotated", "Filtered")) tt <- table(groupProb) tt ################################################### ### chunk number 36: eval=FALSE ################################################### ## #line 752 "vignettes/topGO/inst/doc/topGO.Rnw" ## pValue <- getPvalues(exprs(ALL), classlabel = y, alternative = "greater") ## geneVar <- apply(exprs(ALL), 1, var) ## dd <- data.frame(x = geneVar[allProb], y = log10(pValue[allProb]), groups = groupProb) ## xyplot(y ~ x | groups, data = dd, groups = groups) ################################################### ### chunk number 37: ################################################### #line 760 "vignettes/topGO/inst/doc/topGO.Rnw" pValue <- getPvalues(exprs(ALL), classlabel = y, alternative = "greater") geneVar <- apply(exprs(ALL), 1, var) dd <- data.frame(x = geneVar[allProb], y = log10(pValue[allProb]), groups = groupProb) library(lattice) trellis.device(device = pdf, theme = col.whitebg(), file = "whichProbe.pdf", width = 9, height = 7) legendLab <- paste(names(table(groupProb)), " (#", table(groupProb), ")", sep = "") pP <- xyplot(y ~ x | groups, data = dd, groups = groups, xlab = "Variance", ylab = "Log of p-values", layout = c(2, 2), key = list(text = list(lab = legendLab), points = list(pch = 20, cex = 2, col = Rows(trellis.par.get("superpose.symbol"), 1:3)$col), size = 7, padding.text = 3, x = .65, y = .7, corner = c(0, 0), border = TRUE, cex = 1), panel = function(x, y, ...) { selY <- y <= -2 panel.xyplot(x[selY], y[selY], pch = 2, ...) panel.xyplot(x[!selY], y[!selY], pch = 20, ...) panel.abline(h = -2, lty = 2, col = "black") }) print(pP) dev.off() ################################################### ### chunk number 38: ################################################### #line 806 "vignettes/topGO/inst/doc/topGO.Rnw" description(GOdata) description(GOdata) <- paste(description(GOdata), "Object modified on:", format(Sys.time(), "%d %b %Y"), sep = " ") description(GOdata) ################################################### ### chunk number 39: ################################################### #line 815 "vignettes/topGO/inst/doc/topGO.Rnw" a <- genes(GOdata) ## obtain the list of genes head(a) numGenes(GOdata) ################################################### ### chunk number 40: ################################################### #line 825 "vignettes/topGO/inst/doc/topGO.Rnw" selGenes <- sample(a, 10) gs <- geneScore(GOdata, whichGenes = selGenes) print(gs) ################################################### ### chunk number 41: ################################################### #line 832 "vignettes/topGO/inst/doc/topGO.Rnw" gs <- geneScore(GOdata, whichGenes = selGenes, use.names = FALSE) print(gs) gs <- geneScore(GOdata, use.names = FALSE) str(gs) ################################################### ### chunk number 42: ################################################### #line 841 "vignettes/topGO/inst/doc/topGO.Rnw" sg <- sigGenes(GOdata) str(sg) numSigGenes(GOdata) ################################################### ### chunk number 43: ################################################### #line 851 "vignettes/topGO/inst/doc/topGO.Rnw" .geneList <- geneScore(GOdata, use.names = TRUE) GOdata ## more available genes GOdata <- updateGenes(GOdata, .geneList, topDiffGenes) GOdata ## the available genes are now the feasible genes ################################################### ### chunk number 44: ################################################### #line 862 "vignettes/topGO/inst/doc/topGO.Rnw" graph(GOdata) ## returns the GO graph ug <- usedGO(GOdata) head(ug) ################################################### ### chunk number 45: ################################################### #line 871 "vignettes/topGO/inst/doc/topGO.Rnw" sel.terms <- sample(usedGO(GOdata), 10) num.ann.genes <- countGenesInTerm(GOdata, sel.terms) ## the number of annotated genes num.ann.genes ann.genes <- genesInTerm(GOdata, sel.terms) ## get the annotations head(ann.genes) ################################################### ### chunk number 46: ################################################### #line 883 "vignettes/topGO/inst/doc/topGO.Rnw" ann.score <- scoresInTerm(GOdata, sel.terms) head(ann.score) ann.score <- scoresInTerm(GOdata, sel.terms, use.names = TRUE) head(ann.score) ################################################### ### chunk number 47: ################################################### #line 893 "vignettes/topGO/inst/doc/topGO.Rnw" termStat(GOdata, sel.terms) ################################################### ### chunk number 48: ################################################### #line 974 "vignettes/topGO/inst/doc/topGO.Rnw" goID <- "GO:0044255" gene.universe <- genes(GOdata) go.genes <- genesInTerm(GOdata, goID)[[1]] sig.genes <- sigGenes(GOdata) ################################################### ### chunk number 49: ################################################### #line 984 "vignettes/topGO/inst/doc/topGO.Rnw" my.group <- new("classicCount", testStatistic = GOFisherTest, name = "fisher", allMembers = gene.universe, groupMembers = go.genes, sigMembers = sig.genes) contTable(my.group) runTest(my.group) ################################################### ### chunk number 50: ################################################### #line 1011 "vignettes/topGO/inst/doc/topGO.Rnw" set.seed(123) elim.genes <- sample(go.genes, length(go.genes) / 4) elim.group <- new("elimCount", testStatistic = GOFisherTest, name = "fisher", allMembers = gene.universe, groupMembers = go.genes, sigMembers = sig.genes, elim = elim.genes) contTable(elim.group) runTest(elim.group) ################################################### ### chunk number 51: ################################################### #line 1051 "vignettes/topGO/inst/doc/topGO.Rnw" test.stat <- new("classicCount", testStatistic = GOFisherTest, name = "Fisher test") resultFisher <- getSigGroups(GOdata, test.stat) ################################################### ### chunk number 52: ################################################### #line 1059 "vignettes/topGO/inst/doc/topGO.Rnw" resultFisher ################################################### ### chunk number 53: ################################################### #line 1068 "vignettes/topGO/inst/doc/topGO.Rnw" test.stat <- new("classicScore", testStatistic = GOKSTest, name = "KS tests") resultKS <- getSigGroups(GOdata, test.stat) ################################################### ### chunk number 54: eval=FALSE ################################################### ## #line 1083 "vignettes/topGO/inst/doc/topGO.Rnw" ## test.stat <- new("elimScore", testStatistic = GOKSTest, name = "Fisher test", cutOff = 0.01) ## resultElim <- getSigGroups(GOdata, test.stat) ################################################### ### chunk number 55: ################################################### #line 1090 "vignettes/topGO/inst/doc/topGO.Rnw" test.stat <- new("weightCount", testStatistic = GOFisherTest, name = "Fisher test", sigRatio = "ratio") resultWeight <- getSigGroups(GOdata, test.stat) ################################################### ### chunk number 56: ################################################### #line 1150 "vignettes/topGO/inst/doc/topGO.Rnw" resultFis <- runTest(GOdata, algorithm = "classic", statistic = "fisher") ################################################### ### chunk number 57: eval=FALSE ################################################### ## #line 1159 "vignettes/topGO/inst/doc/topGO.Rnw" ## weight01.fisher <- runTest(GOdata, statistic = "fisher") ## weight01.t <- runTest(GOdata, algorithm = "weight01", statistic = "t") ## elim.ks <- runTest(GOdata, algorithm = "elim", statistic = "ks") ## ## weight.ks <- runTest(GOdata, algorithm = "weight", statistic = "ks") #will not work!!! ################################################### ### chunk number 58: ################################################### #line 1171 "vignettes/topGO/inst/doc/topGO.Rnw" whichTests() whichAlgorithms() ################################################### ### chunk number 59: ################################################### #line 1207 "vignettes/topGO/inst/doc/topGO.Rnw" pvalFis <- score(resultFis) head(pvalFis) hist(pvalFis, 50, xlab = "p-values") ################################################### ### chunk number 60: ################################################### #line 1217 "vignettes/topGO/inst/doc/topGO.Rnw" head(score(resultWeight)) ################################################### ### chunk number 61: ################################################### #line 1227 "vignettes/topGO/inst/doc/topGO.Rnw" pvalWeight <- score(resultWeight, whichGO = names(pvalFis)) head(pvalWeight) cor(pvalFis, pvalWeight) ################################################### ### chunk number 62: ################################################### #line 1237 "vignettes/topGO/inst/doc/topGO.Rnw" geneData(resultWeight) ################################################### ### chunk number 63: ################################################### #line 1250 "vignettes/topGO/inst/doc/topGO.Rnw" allRes <- GenTable(GOdata, classic = resultFis, KS = resultKS, weight = resultWeight, orderBy = "weight", ranksOf = "classic", topNodes = 20) ################################################### ### chunk number 64: ################################################### #line 1262 "vignettes/topGO/inst/doc/topGO.Rnw" if(require(xtable)) print(xtable(apply(allRes, 2, as.character)), floating = FALSE) ################################################### ### chunk number 65: ################################################### #line 1285 "vignettes/topGO/inst/doc/topGO.Rnw" goID <- allRes[1, "GO.ID"] print(showGroupDensity(GOdata, goID, ranks = TRUE)) ################################################### ### chunk number 66: ################################################### #line 1310 "vignettes/topGO/inst/doc/topGO.Rnw" goID <- allRes[10, "GO.ID"] gt <- printGenes(GOdata, whichTerms = goID, chip = affyLib, numChar = 40) ################################################### ### chunk number 67: ################################################### #line 1317 "vignettes/topGO/inst/doc/topGO.Rnw" if(require(xtable)) print(xtable(gt), floating = FALSE) ################################################### ### chunk number 68: eval=FALSE ################################################### ## #line 1343 "vignettes/topGO/inst/doc/topGO.Rnw" ## showSigOfNodes(GOdata, score(resultFis), firstSigNodes = 5, useInfo = 'all') ## showSigOfNodes(GOdata, score(resultWeight), firstSigNodes = 5, useInfo = 'def') ################################################### ### chunk number 69: ################################################### #line 1348 "vignettes/topGO/inst/doc/topGO.Rnw" printGraph(GOdata, resultFis, firstSigNodes = 5, fn.prefix = "tGO", useInfo = "all", pdfSW = TRUE) printGraph(GOdata, resultWeight, firstSigNodes = 5, fn.prefix = "tGO", useInfo = "def", pdfSW = TRUE) ################################################### ### chunk number 70: eval=FALSE ################################################### ## #line 1380 "vignettes/topGO/inst/doc/topGO.Rnw" ## printGraph(GOdata, resultWeight, firstSigNodes = 10, resultFis, fn.prefix = "tGO", useInfo = "def") ## printGraph(GOdata, resultElim, firstSigNodes = 15, resultFis, fn.prefix = "tGO", useInfo = "all") ################################################### ### chunk number 71: ################################################### #line 1391 "vignettes/topGO/inst/doc/topGO.Rnw" toLatex(sessionInfo())