inst/doc/GOstatsHyperG.R

In GOstats: Tools for manipulating GO and microarrays

### R code from vignette source 'GOstatsHyperG.Rnw'

###################################################
### code chunk number 1: Setup
###################################################
library("ALL")
library("hgu95av2.db")
library("GO.db")
library("annotate")
library("genefilter")
library("GOstats")
library("RColorBrewer")
library("xtable")
library("Rgraphviz")


###################################################
### code chunk number 2: 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)


###################################################
### code chunk number 3: bcrAblOrNegSubset
###################################################
data(ALL, package="ALL")
## For this data we can have ALL1/AF4 or BCR/ABL
subsetType <- "ALL1/AF4"

## Subset of interest: 37BRC/ABL + 42NEG = 79 samples
Bcell <- grep("^B", as.character(ALL$BT))
bcrAblOrNegIdx <- which(as.character(ALL$mol) %in% c("NEG", subsetType))

bcrAblOrNeg <- ALL[, intersect(Bcell, bcrAblOrNegIdx)]
bcrAblOrNeg$mol.biol = factor(bcrAblOrNeg$mol.biol)


###################################################
### code chunk number 4: nsFiltering-noEntrez
###################################################
## Remove genes that have no entrezGene id
entrezIds <- mget(featureNames(bcrAblOrNeg), envir=hgu95av2ENTREZID)
haveEntrezId <- names(entrezIds)[sapply(entrezIds, function(x) !is.na(x))]
numNoEntrezId <- length(featureNames(bcrAblOrNeg)) - length(haveEntrezId)
bcrAblOrNeg <- bcrAblOrNeg[haveEntrezId, ]


###################################################
### code chunk number 5: nsFiltering-noGO
###################################################
## Remove genes with no GO mapping
haveGo <- sapply(mget(featureNames(bcrAblOrNeg), hgu95av2GO),
                 function(x) {
                     if (length(x) == 1 && is.na(x)) 
                       FALSE 
                     else TRUE
                 })
numNoGO <- sum(!haveGo)
bcrAblOrNeg <- bcrAblOrNeg[haveGo, ]


###################################################
### code chunk number 6: nsFiltering-IQR
###################################################
## Non-specific filtering based on IQR
iqrCutoff <- 0.5
bcrAblOrNegIqr <- apply(exprs(bcrAblOrNeg), 1, IQR)
selected <- bcrAblOrNegIqr > iqrCutoff

## Drop those that are on the Y chromosome
## because there is an imbalance of men and women by group
chrN <- mget(featureNames(bcrAblOrNeg), envir=hgu95av2CHR)
onY <- sapply(chrN, function(x) any(x=="Y"))
onY[is.na(onY)] <- FALSE
selected <- selected & !onY

nsFiltered <- bcrAblOrNeg[selected, ]


###################################################
### code chunk number 7: nsFiltering-unique
###################################################
numNsWithDups <- length(featureNames(nsFiltered))
nsFilteredIqr <- bcrAblOrNegIqr[selected]
uniqGenes <- findLargest(featureNames(nsFiltered), nsFilteredIqr, 
                         "hgu95av2")
nsFiltered <- nsFiltered[uniqGenes, ]
numSelected <- length(featureNames(nsFiltered))

##set up some colors
BCRcols = ifelse(nsFiltered$mol == subsetType, "goldenrod", "skyblue")
cols = brewer.pal(10, "RdBu")


###################################################
### code chunk number 8: defineGeneUniverse
###################################################
## Define gene universe based on results of non-specific filtering
affyUniverse <- featureNames(nsFiltered)
entrezUniverse <- unlist(mget(affyUniverse, hgu95av2ENTREZID))
if (any(duplicated(entrezUniverse)))
  stop("error in gene universe: can't have duplicate Entrez Gene Ids")

## Also define an alternate universe based on the entire chip
chipAffyUniverse <- featureNames(bcrAblOrNeg)
chipEntrezUniverse <- mget(chipAffyUniverse, hgu95av2ENTREZID)
chipEntrezUniverse <- unique(unlist(chipEntrezUniverse)) 


###################################################
### code chunk number 9: parametric1
###################################################
ttestCutoff <- 0.05
ttests = rowttests(nsFiltered, "mol.biol")

smPV = ttests$p.value < ttestCutoff

pvalFiltered <- nsFiltered[smPV, ]
selectedEntrezIds <- unlist(mget(featureNames(pvalFiltered),
                                 hgu95av2ENTREZID))


###################################################
### code chunk number 10: withYourData1 (eval = FALSE)
###################################################
## entrezUniverse <- unlist(mget(featureNames(yourData), 
##                               hgu95av2ENTREZID))
## if (any(duplicated(entrezUniverse)))
##   stop("error in gene universe: can't have duplicate Entrez Gene Ids")
## 
## pvalFiltered <- yourData[hasSmallPvalue, ]
## selectedEntrezIds <- unlist(mget(featureNames(pvalFiltered),
##                                  hgu95av2ENTREZID))


###################################################
### code chunk number 11: standardHyperGeo
###################################################
hgCutoff <- 0.001
params <- new("GOHyperGParams",
              geneIds=selectedEntrezIds,
              universeGeneIds=entrezUniverse,
              annotation="hgu95av2.db",
              ontology="BP",
              pvalueCutoff=hgCutoff,
              conditional=FALSE,
              testDirection="over")



###################################################
### code chunk number 12: condHyperGeo
###################################################
paramsCond <- params
conditional(paramsCond) <- TRUE


###################################################
### code chunk number 13: standardHGTEST
###################################################
hgOver <- hyperGTest(params)
hgCondOver <- hyperGTest(paramsCond)



###################################################
### code chunk number 14: HGTestAns
###################################################
hgOver
hgCondOver


###################################################
### code chunk number 15: summaryEx
###################################################
df <- summary(hgOver)
names(df)                               # the columns
dim(summary(hgOver, pvalue=0.1))
dim(summary(hgOver, categorySize=10))


###################################################
### code chunk number 16: resultAccessors
###################################################
pvalues(hgOver)[1:3]

oddsRatios(hgOver)[1:3]

expectedCounts(hgOver)[1:3]

geneCounts(hgOver)[1:3]
universeCounts(hgOver)[1:3]

length(geneIds(hgOver))
length(geneIdUniverse(hgOver)[[3]])

## GOHyperGResult _only_
## (NB: edges go from parent to child)
goDag(hgOver)

geneMappedCount(hgOver)
universeMappedCount(hgOver)

conditional(hgOver)
testDirection(hgOver)
testName(hgOver)



###################################################
### code chunk number 17: htmlReportExample
###################################################
htmlReport(hgCondOver, file="ALL_hgco.html")


###################################################
### code chunk number 18: helperFunc
###################################################
sigCategories <- function(res, p) {
    if (missing(p))
      p <- pvalueCutoff(res)
    pv <- pvalues(res)
    goIds <- names(pv[pv < p])
    goIds
}


###################################################
### code chunk number 19: plotFuns
###################################################
coloredGoPlot <- function(ccMaxGraph, hgOver, hgCondOver) {
    nodeColors <- sapply(nodes(ccMaxGraph),
                         function(n) {
                             if (n %in% sigCategories(hgCondOver))
                               "dark red"
                             else if (n %in% sigCategories(hgOver))
                               "pink"
                             else 
                               "gray"
                         })
    nattr <- makeNodeAttrs(ccMaxGraph,
                           label=nodes(ccMaxGraph),
                           shape="ellipse",
                           fillcolor=nodeColors,
                           fixedsize=FALSE)
    plot(ccMaxGraph, nodeAttrs=nattr)
}

getMaxConnCompGraph <- function(hgOver, hgCondOver) {
    ##uGoDagRev <- ugraph(goDag(hgOver))
    sigNodes <- sigCategories(hgOver)
    ##adjNodes <- unlist(adj(uGoDagRev, sigNodes))
    adjNodes <- unlist(adj(goDag(hgOver), sigNodes))
    displayNodes <- unique(c(sigNodes, adjNodes))
    displayGraph <- subGraph(displayNodes, goDag(hgOver))
    cc <- connComp(displayGraph)
    ccSizes <- listLen(cc)
    ccMaxIdx <- which(ccSizes == max(ccSizes))
    ccMaxGraph <- subGraph(cc[[ccMaxIdx]], displayGraph)
    ccMaxGraph 
}



###################################################
### code chunk number 20: info
###################################################
toLatex(sessionInfo())