Code examples for Chapter 14: 'Hypergeometric Testing used for Gene Set Enrichment Analysis'

## Chunk 1
`library("Biobase") library("ALL") library("hgu95av2.db") library("annotate") library("genefilter") library("GOstats") library("RColorBrewer") library("xtable") library("Rgraphviz")`

## Chunk 2
`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)`

(solution chunk) ## Chunk 3
`numSamp = length(ALL_bcrneg$mol.biol) table(ALL_bcrneg$mol.biol)`

(solution chunk) ## Chunk 4
`annotation(ALL_bcrneg) length(featureNames(ALL_bcrneg))`

## Chunk 5
`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`

(solution chunk) ## Chunk 6
`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 7
`hasSymbol = sapply(mget(featureNames(ALLfilt_bcrneg), envir=hgu95av2SYMBOL), function(x) !(length(x) == 1 && is.na(x))) ALLfilt_bcrneg = ALLfilt_bcrneg[hasSymbol, ]`

## Chunk 8
`affyUniverse = featureNames(ALLfilt_bcrneg) entrezUniverse = unlist(mget(affyUniverse, hgu95av2ENTREZID))`

(solution chunk) ## Chunk 9
`## an alternate universe based on the entire chip chipAffyUniverse = featureNames(ALLfilt_bcrneg) chipEntrezUniverse = mget(chipAffyUniverse, hgu95av2ENTREZID) chipEntrezUniverse = unique(unlist(chipEntrezUniverse))`

## Chunk 10
`ttestCutoff = 0.05 ttests = rowttests(ALLfilt_bcrneg, "mol.biol") smPV = ttests$p.value < ttestCutoff pvalFiltered = ALLfilt_bcrneg[smPV, ] selectedEntrezIds = unlist(mget(featureNames(pvalFiltered), hgu95av2ENTREZID))`

(solution chunk) ## Chunk 11
`sumpv = sum(smPV)`

## Chunk 12
`## ## 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 13
`hgCutoff = 0.001 params = new("GOHyperGParams", geneIds=selectedEntrezIds, universeGeneIds=entrezUniverse, annotation="hgu95av2.db", ontology="BP", pvalueCutoff=hgCutoff, conditional=FALSE, testDirection="over")`

## Chunk 14
`hgOver = hyperGTest(params)`

## Chunk 15
`hgOver`

(solution chunk) ## Chunk 16
`df = summary(hgOver) names(df)`

(solution chunk) ## Chunk 17
`df = summary(hgOver, pvalue=0.05, categorySize=350) nrow(df)`

(solution chunk) ## Chunk 18
`## ? HyperGResult-accessors`

## Chunk 19
`htmlReport(hgOver, file="ALL_hgo.html")`

(solution chunk) ## Chunk 20
`## browseURL("ALL_hgo.html")`

## Chunk 21
`sigSub = termGraphs(hgOver)`

(solution chunk) ## Chunk 22
`numG = length(sigSub) sizes = sapply(sigSub, numNodes) sizes`

## Chunk 23
`plotGOTermGraph(sigSub[[1]], hgOver, max.nchar=100)`

## Chunk 24
`paramsCond = params conditional(paramsCond) = TRUE`

## Chunk 25
`hgCond = hyperGTest(paramsCond)`

## Chunk 26
`hgCond`

(solution chunk) ## Chunk 27
`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")]`

(solution chunk) ## Chunk 28
`stdIds = sigCategories(hgOver) condIds = sigCategories(hgCond) setdiff(stdIds, condIds)`

## Chunk 29
`terms = 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="") }) df[terms, ]`

(solution chunk) ## Chunk 30
`params = new("ChrMapHyperGParams", conditional=FALSE, testDirection="over", universeGeneIds=entrezUniverse, geneIds=selectedEntrezIds, annotation="hgu95av2", pvalueCutoff=0.05) paramsCond = params conditional(paramsCond) = TRUE`

(solution chunk) ## Chunk 31
`hgans = hyperGTest(params) hgansCond = hyperGTest(paramsCond)`

(solution chunk) ## Chunk 32
`summary(hgans, categorySize=10)`

(solution chunk) ## Chunk 33
`kparams = new("KEGGHyperGParams", geneIds=selectedEntrezIds, universeGeneIds=entrezUniverse, annotation="hgu95av2", pvalueCutoff=0.05, testDirection="over") kans = hyperGTest(kparams)`

(solution chunk) ## Chunk 34
`summary(kans) kparamsUnder = kparams testDirection(kparamsUnder) = "under"`

(solution chunk) ## Chunk 35
`kansUnder = hyperGTest(kparamsUnder)`

(solution chunk) ## Chunk 36
`summary(kansUnder)`

(solution chunk) ## Chunk 37
`pparams = new("PFAMHyperGParams", geneIds=selectedEntrezIds, universeGeneIds=entrezUniverse, annotation="hgu95av2", pvalueCutoff=hgCutoff, testDirection="over") pans = hyperGTest(pparams)`

(solution chunk) ## Chunk 38
`summary(pans)`