###################################################
### chunk: Setup
###################################################
library("Biobase")
library("ALL")
library("hgu95av2.db")
library("annotate")
library("genefilter")
library("GOstats")
library("RColorBrewer")
library("xtable")
library("Rgraphviz")


###################################################
### chunk: 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: sol1
###################################################
numSamp = length(ALL_bcrneg$mol.biol)
table(ALL_bcrneg$mol.biol)


###################################################
### chunk: sol2
###################################################
annotation(ALL_bcrneg)
length(featureNames(ALL_bcrneg))


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


###################################################
### chunk: defineGeneUniverse
###################################################
affyUniverse = featureNames(ALLfilt_bcrneg)
entrezUniverse = unlist(mget(affyUniverse, 
    hgu95av2ENTREZID))


###################################################
### chunk: altUniv
###################################################
## an alternate universe based on the entire chip
chipAffyUniverse = featureNames(ALLfilt_bcrneg)
chipEntrezUniverse = mget(chipAffyUniverse, hgu95av2ENTREZID)
chipEntrezUniverse = unique(unlist(chipEntrezUniverse))


###################################################
### chunk: 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: sumPV
###################################################
sumpv = sum(smPV)


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


###################################################
### chunk: standardHGTEST
###################################################
hgOver = hyperGTest(params)


###################################################
### chunk: HGTestAns
###################################################
hgOver


###################################################
### chunk: summaryNames
###################################################
df = summary(hgOver)
names(df)


###################################################
### chunk: summary2
###################################################
df = summary(hgOver, pvalue=0.05, categorySize=350)
nrow(df)


###################################################
### chunk: achelp eval=FALSE
###################################################
## ? HyperGResult-accessors


###################################################
### chunk: htmlReportExample
###################################################
htmlReport(hgOver, file="ALL_hgo.html")


###################################################
### chunk: htmlBrowse eval=FALSE
###################################################
## browseURL("ALL_hgo.html")


###################################################
### chunk: termGraphs
###################################################
sigSub = termGraphs(hgOver)


###################################################
### chunk: sigSub
###################################################
numG = length(sigSub)
sizes = sapply(sigSub, numNodes)
sizes


###################################################
### chunk: cellComFig
###################################################
plotGOTermGraph(sigSub[[1]], hgOver, max.nchar=100)


###################################################
### chunk: condHyperGeo
###################################################
paramsCond = params
conditional(paramsCond) = TRUE


###################################################
### chunk: condhgRun
###################################################
hgCond = hyperGTest(paramsCond)


###################################################
### chunk: condhgResult
###################################################
hgCond


###################################################
### chunk: 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: compareResults
###################################################
stdIds = sigCategories(hgOver)
condIds = sigCategories(hgCond)
setdiff(stdIds, condIds)


###################################################
### chunk: cellcomcond
###################################################
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, ]


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


###################################################
### chunk: basicTest
###################################################
hgans = hyperGTest(params)
hgansCond = hyperGTest(paramsCond)


###################################################
### chunk: ChromResult
###################################################
summary(hgans, categorySize=10)


###################################################
### chunk: kegg1
###################################################
kparams = new("KEGGHyperGParams",
    geneIds=selectedEntrezIds, 
    universeGeneIds=entrezUniverse,
    annotation="hgu95av2",
    pvalueCutoff=0.05,
    testDirection="over")

kans = hyperGTest(kparams)


###################################################
### chunk: kegg2
###################################################
summary(kans)
kparamsUnder = kparams
testDirection(kparamsUnder) = "under"


###################################################
### chunk: keggUnder
###################################################
kansUnder = hyperGTest(kparamsUnder)


###################################################
### chunk: kegg3
###################################################
summary(kansUnder)


###################################################
### chunk: pfam1
###################################################
pparams = new("PFAMHyperGParams",
    geneIds=selectedEntrezIds,
    universeGeneIds=entrezUniverse,
    annotation="hgu95av2",
    pvalueCutoff=hgCutoff,
    testDirection="over")


pans = hyperGTest(pparams)


###################################################
### chunk: pfam2
###################################################
summary(pans)