################################################### ### chunk number 1: load0 ################################################### library("BiocCaseStudies") ################################################### ### chunk number 2: loadAllLibsFirst ################################################### library("Biobase") library("ALL") library("hgu95av2.db") library("annotate") library("genefilter") library("GOstats") library("RColorBrewer") library("xtable") library("Rgraphviz") ################################################### ### chunk number 3: setup1 ################################################### ## For this data we can have ALL1/AF4 or BCR/ABL ## This is redefined below to make the lab easier to follow, ## but we want it pre-defined so we can use this in Sexpr{}, etc. subsetType = "BCR/ABL" varCut = 0.5 ################################################### ### chunk number 4: Setup ################################################### library("Biobase") library("ALL") library("hgu95av2.db") library("annotate") library("genefilter") library("GOstats") library("RColorBrewer") library("xtable") library("Rgraphviz") ################################################### ### chunk number 5: universeSizeVsPvalue ################################################### hg_tester = function(size) { numFound = 10 numDrawn = 400 numAtCat = 40 numNotAtCat = size - numAtCat phyper(numFound-1, numAtCat, numNotAtCat, numDrawn, lower.tail=FALSE) } pv1000 = hg_tester(1000) pv5000 = hg_tester(5000) ################################################### ### chunk number 6: subset ################################################### data(ALL) bcell = grep("^B", as.character(ALL$BT)) types = c("NEG", "BCR/ABL") moltyp = which(as.character(ALL$mol.biol) %in% types) ALL_bcrneg = ALL[, intersect(bcell, moltyp)] ALL_bcrneg$mol.biol = factor(ALL_bcrneg$mol.biol) ################################################### ### chunk number 7: sol1 ################################################### numSamp = length(ALL_bcrneg$mol.biol) table(ALL_bcrneg$mol.biol) ################################################### ### chunk number 8: sol2 ################################################### annotation(ALL_bcrneg) length(featureNames(ALL_bcrneg)) ################################################### ### chunk number 9: nsfilter ################################################### varCut = 0.5 filt_bcrneg = nsFilter(ALL_bcrneg, require.entrez=TRUE, require.GOBP=TRUE, remove.dupEntrez=TRUE, var.func=IQR, var.cutoff=varCut, feature.exclude="^AFFX") names(filt_bcrneg) ALLfilt_bcrneg = filt_bcrneg$eset ################################################### ### chunk number 10: nsFiltering-Y ################################################### chrN = mget(featureNames(ALLfilt_bcrneg), envir=hgu95av2CHR) onY = sapply(chrN, function(x) any(x == "Y")) onY[is.na(onY)] = FALSE ALLfilt_bcrneg = ALLfilt_bcrneg[!onY, ] ################################################### ### chunk number 11: removeMissingSYMBOL ################################################### hasSymbol = sapply(mget(featureNames(ALLfilt_bcrneg), envir=hgu95av2SYMBOL), function(x) !(length(x) == 1 && is.na(x))) ALLfilt_bcrneg = ALLfilt_bcrneg[hasSymbol, ] ################################################### ### chunk number 12: defineGeneUniverse ################################################### affyUniverse = featureNames(ALLfilt_bcrneg) entrezUniverse = unlist(mget(affyUniverse, hgu95av2ENTREZID)) ################################################### ### chunk number 13: universeCheck ################################################### if (any(duplicated(entrezUniverse))) stop("error in gene universe: can't have duplicate Entrez Gene Ids") ################################################### ### chunk number 14: altUniv ################################################### ## an alternate universe based on the entire chip chipAffyUniverse = featureNames(ALLfilt_bcrneg) chipEntrezUniverse = mget(chipAffyUniverse, hgu95av2ENTREZID) chipEntrezUniverse = unique(unlist(chipEntrezUniverse)) ################################################### ### chunk number 15: parametric1 ################################################### ttestCutoff = 0.05 ttests = rowttests(ALLfilt_bcrneg, "mol.biol") smPV = ttests$p.value < ttestCutoff pvalFiltered = ALLfilt_bcrneg[smPV, ] selectedEntrezIds = unlist(mget(featureNames(pvalFiltered), hgu95av2ENTREZID)) ################################################### ### chunk number 16: sumPV ################################################### sumpv = sum(smPV) ################################################### ### chunk number 17: withYourData1 eval=FALSE ################################################### ## ## if you are following along with your own data... ## entrezUniverse = unlist(mget(featureNames(yourData), ## hgu95av2ENTREZID)) ## pvalFiltered = yourData[hasSmallPvalue, ] ## selectedEntrezIds = unlist(mget(featureNames(pvalFiltered), ## hgu95av2ENTREZID)) ################################################### ### chunk number 18: standardHyperGeo ################################################### hgCutoff = 0.001 params = new("GOHyperGParams", geneIds=selectedEntrezIds, universeGeneIds=entrezUniverse, annotation="hgu95av2.db", ontology="BP", pvalueCutoff=hgCutoff, conditional=FALSE, testDirection="over") ################################################### ### chunk number 19: standardHGTEST ################################################### hgOver = hyperGTest(params) ################################################### ### chunk number 20: HGTestAns ################################################### hgOver ################################################### ### chunk number 21: summaryEx ################################################### df = summary(hgOver, categorySize=100) if (nrow(df) < 6) stop("unexpected hyperGTest summary result, check categorySize?") df = df[1:6, ] df$Term = substr(df$Term, 0, 30) sizeOrd = order(df$Size, decreasing=TRUE) tab = xtable(df[sizeOrd, c("Count", "Size", "Term")], label="tab:hgOver", caption=paste("Top six overrepresented GO terms with ", "at least 100 gene annotations identified by", "\\Rfunction{hyperGTest}.")) print(tab) ################################################### ### chunk number 22: summaryNames ################################################### df = summary(hgOver) names(df) ################################################### ### chunk number 23: summary2 ################################################### df = summary(hgOver, pvalue=0.05, categorySize=350) nrow(df) ################################################### ### chunk number 24: achelp eval=FALSE ################################################### ## ? HyperGResult-accessors ################################################### ### chunk number 25: htmlReportExample ################################################### htmlReport(hgOver, file="ALL_hgo.html") ################################################### ### chunk number 26: htmlBrowse eval=FALSE ################################################### ## browseURL("ALL_hgo.html") ################################################### ### chunk number 27: termGraphs ################################################### sigSub = termGraphs(hgOver) ################################################### ### chunk number 28: sigSub ################################################### numG = length(sigSub) sizes = sapply(sigSub, numNodes) sizes ################################################### ### chunk number 29: cellComFig ################################################### plotGOTermGraph(sigSub[[1]], hgOver, max.nchar=100) ################################################### ### chunk number 30: condHyperGeo ################################################### paramsCond = params conditional(paramsCond) = TRUE ################################################### ### chunk number 31: condhgRun ################################################### hgCond = hyperGTest(paramsCond) ################################################### ### chunk number 32: condhgResult ################################################### hgCond ################################################### ### chunk number 33: condhgSummary ################################################### dfcond = summary(hgCond, categorySize=50) ## trim the term names for display purposes trimTerm = function(x) { if (nchar(x) <= 20) x else paste(substr(x, 1, 20), "...", sep="") } dfcond$Term = sapply(dfcond$Term, trimTerm) sizeOrd = order(dfcond$Size, decreasing=TRUE) dfcond[sizeOrd, c("Count", "Size", "Term")] ################################################### ### chunk number 34: compareResults ################################################### stdIds = sigCategories(hgOver) condIds = sigCategories(hgCond) setdiff(stdIds, condIds) ################################################### ### chunk number 35: cellcomcond ################################################### termIDs = nodes(sigSub[[1]]) df = summary(hgCond, pvalue=0.5)[ , c("Term", "Pvalue")] df$Pvalue = round(df$Pvalue, 3) df$Term = sapply(df$Term, function(x) { if (nchar(x) <= 20) x else paste(substr(x, 1, 20), "...", sep="") }) terms = mget(as.character(termIDs), GOTERM, ifnotfound=NA) terms = as.character(unlist(lapply(terms, function(x){x@Term}))) index = df$Term %in% terms df[index, ] ################################################### ### chunk number 36: basicParams ################################################### params = new("ChrMapHyperGParams", conditional=FALSE, testDirection="over", universeGeneIds=entrezUniverse, geneIds=selectedEntrezIds, annotation="hgu95av2", pvalueCutoff=0.05) paramsCond = params conditional(paramsCond) = TRUE ################################################### ### chunk number 37: basicTest ################################################### hgans = hyperGTest(params) hgansCond = hyperGTest(paramsCond) ################################################### ### chunk number 38: ChromResult ################################################### summary(hgans, categorySize=10) ################################################### ### chunk number 39: kegg1 ################################################### kparams = new("KEGGHyperGParams", geneIds=selectedEntrezIds, universeGeneIds=entrezUniverse, annotation="hgu95av2", pvalueCutoff=0.05, testDirection="over") kans = hyperGTest(kparams) ################################################### ### chunk number 40: kegg2 ################################################### summary(kans) kparamsUnder = kparams testDirection(kparamsUnder) = "under" ################################################### ### chunk number 41: keggUnder ################################################### kansUnder = hyperGTest(kparamsUnder) ################################################### ### chunk number 42: kegg3 ################################################### summary(kansUnder) ################################################### ### chunk number 43: pfam1 ################################################### pparams = new("PFAMHyperGParams", geneIds=selectedEntrezIds, universeGeneIds=entrezUniverse, annotation="hgu95av2", pvalueCutoff=hgCutoff, testDirection="over") pans = hyperGTest(pparams) ################################################### ### chunk number 44: pfam2 ################################################### summary(pans) ################################################### ### chunk number 45: sessionInfo ################################################### mySessionInfo()