###################################################
### chunk number 1: options
###################################################
options(width=50, digits=2)


###################################################
### chunk number 2: ALL
###################################################
library("ALL")
data("ALL")


###################################################
### chunk number 3: bcrabl
###################################################
bcell = grep("^B", as.character(ALL$BT)) # B-cell ALL
moltyp = which(as.character(ALL$mol.biol) 
    %in% c("NEG", "BCR/ABL"))
bcrneg = ALL[, intersect(bcell, moltyp)]
bcrneg$mol.biol = factor(bcrneg$mol.biol) # drop unused levels


###################################################
### chunk number 4: bcrabl-result
###################################################
dim(bcrneg)


###################################################
### chunk number 5: nsFilter
###################################################
library("genefilter")
bcrneg_filt = nsFilter(bcrneg, var.cutoff=0.5, require.GOBP=TRUE)$eset
dim(bcrneg_filt)


###################################################
### chunk number 6: rowtttest
###################################################
rtt <- rowttests(bcrneg_filt, "mol.biol")
rttPrb <- rtt$p.value
names(rttPrb) <- featureNames(bcrneg_filt)
tThresh <- rttPrb < 0.05


###################################################
### chunk number 7: GOselected
###################################################
ids <- featureNames(bcrneg_filt)
map <- hgu95av2ENTREZID
universe <- unlist(mget(ids, map))
selected <- unlist(mget(ids[tThresh], map))


###################################################
### chunk number 8: GOTERM
###################################################
library(GO.db)
GOTERM[["GO:0006468"]]


###################################################
### chunk number 9: tally eval=FALSE
###################################################
## res <- eapply(hgu95av2GO[names(universe)], function(elt) {
##     any(sapply(elt, "[[", "GOID") ==  "GO:0006468")
## })
## df <- data.frame(EntrezId=universe,
##                  Selected=universe %in% selected,
##                  InGO=unlist(res))
## xtabs(~Selected+InGO, df)


###################################################
### chunk number 10: HyperFormulation
###################################################
library(Category)
library(GOstats)
params = new("GOHyperGParams", geneIds=selected,
  universeGeneIds=universe, annotation=annotation(bcrneg_filt),
  ontology="BP", pvalueCutoff=0.001,
  conditional=FALSE, testDirection="over")
(overRepresented = hyperGTest(params))


###################################################
### chunk number 11: HyperInterpretation
###################################################
head(summary(overRepresented), n=3)


###################################################
### chunk number 12: HyperReport eval=FALSE
###################################################
## fl <- tempfile()
## htmlReport(overRepresented, file=fl)
## browseURL(fl)