###################################################
### chunk number 1: 
###################################################
options(width = 70)


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


###################################################
### chunk number 3: 
###################################################
BPterms <- ls(GOBPTerm)
str(BPterms)


###################################################
### chunk number 4: 
###################################################
affyLib <- annotation(ALL)
library(package = affyLib, character.only = TRUE)


###################################################
### chunk number 5: 
###################################################
library(genefilter)
f1 <- pOverA(0.25, log2(100))
f2 <- function(x) (IQR(x) > 0.5)
ff <- filterfun(f1, f2)
eset <- ALL[genefilter(ALL, ff), ]


###################################################
### chunk number 6: 
###################################################
geneNames <- featureNames(eset)
length(geneNames)


###################################################
### chunk number 7: 
###################################################
myInterestedGenes <- sample(geneNames, 100)
geneList <- factor(as.integer(geneNames %in% myInterestedGenes))
names(geneList) <- geneNames
str(geneList)


###################################################
### chunk number 8: 
###################################################
GOdata <- new("topGOdata", ontology = "MF", allGenes = geneList,
              annot = annFUN.hgu, affyLib = affyLib)


###################################################
### chunk number 9: 
###################################################
GOdata


###################################################
### chunk number 10: 
###################################################
y <- as.integer(sapply(eset$BT, function(x) return(substr(x, 1, 1) == 'T')))
str(y)


###################################################
### chunk number 11: 
###################################################
geneList <- getPvalues(exprs(eset), classlabel = y, alternative = "greater")
hist(geneList, br = 50)


###################################################
### chunk number 12: 
###################################################
par(mfrow = c(1, 2))
hist(geneList, br = 50, xlab = "p-values", main = "")
hist(geneList[geneList < 0.1], br = 50, xlab = "p-values < 0.1", main ="")


###################################################
### chunk number 13: 
###################################################
topDiffGenes <- function(allScore) {
  return(allScore < 0.01)
}
x <- topDiffGenes(geneList)
sum(x) ## the number of selected genes


###################################################
### chunk number 14: 
###################################################
GOdata <- new("topGOdata", ontology = "BP", allGenes = geneList, geneSel = topDiffGenes,
              description = "GO analysis of ALL data based on diff. expression.", 
              annot = annFUN.hgu, affyLib = affyLib)


###################################################
### chunk number 15: 
###################################################
description(GOdata)
description(GOdata) <- paste(description(GOdata), "Object modified on:", format(Sys.time(), "%d %b %Y"), sep = " ")
description(GOdata)


###################################################
### chunk number 16: 
###################################################
a <- genes(GOdata) ## obtain the list of genes
str(a)
numGenes(GOdata)


###################################################
### chunk number 17: 
###################################################
selGenes <- sample(a, 10)
gs <- geneScore(GOdata, whichGenes = selGenes) 
print(gs)


###################################################
### chunk number 18: 
###################################################
gs <- geneScore(GOdata, whichGenes = selGenes, use.names = FALSE)
print(gs)

gs <- geneScore(GOdata, use.names = FALSE)
str(gs)


###################################################
### chunk number 19: 
###################################################
sg <- sigGenes(GOdata)
str(sg)
numSigGenes(GOdata)


###################################################
### chunk number 20: 
###################################################
.geneList <- geneScore(GOdata, use.names = TRUE)
GOdata ## more available genes
GOdata <- updateGenes(GOdata, .geneList, topDiffGenes)
GOdata ## the available genes are now the feasible genes


###################################################
### chunk number 21: 
###################################################
graph(GOdata) ## returns the GO graph

ug <- usedGO(GOdata)
str(ug) 


###################################################
### chunk number 22: 
###################################################
sel.terms <- sample(usedGO(GOdata), 10)

num.ann.genes <- countGenesInTerm(GOdata, sel.terms) ## the number of annotated genes
num.ann.genes

ann.genes <- genesInTerm(GOdata, sel.terms) ## get the annotations
str(ann.genes)


###################################################
### chunk number 23: 
###################################################
ann.score <- scoresInTerm(GOdata, sel.terms)
str(ann.score)
ann.score <- scoresInTerm(GOdata, sel.terms, use.names = TRUE)
str(ann.score)


###################################################
### chunk number 24: 
###################################################
termStat(GOdata, sel.terms)


###################################################
### chunk number 25: 
###################################################
test.stat <- new("classicCount", testStatistic = GOFisherTest, name = "Fisher test")
resultFis <- getSigGroups(GOdata, test.stat)


###################################################
### chunk number 26: 
###################################################
test.stat <- new("classicScore", testStatistic = GOKSTest, name = "KS tests")
resultKS <- getSigGroups(GOdata, test.stat)


###################################################
### chunk number 27: 
###################################################
test.stat <- new("elimCount", testStatistic = GOFisherTest, name = "Fisher test", cutOff = 0.01)
resultElim <- getSigGroups(GOdata, test.stat)


###################################################
### chunk number 28: 
###################################################
test.stat <- new("weightCount", testStatistic = GOFisherTest, name = "Fisher test", sigRatio = "ratio")
resultWeight <- getSigGroups(GOdata, test.stat)


###################################################
### chunk number 29: 
###################################################
l <- list(classic = score(resultFis),
          KS = score(resultKS),
          elim = score(resultElim),
          weight = score(resultWeight))

allRes <- genTable(GOdata, l, orderBy = "weight", ranksOf = "classic", top = 20)


###################################################
### chunk number 30: 
###################################################



###################################################
### chunk number 31: 
###################################################
if(require(xtable))
  print(xtable(apply(allRes, 2, as.character)), floating = FALSE)


###################################################
### chunk number 32: 
###################################################
par(mfrow = c(2, 2))
for(nn in names(l)) {
  p.val <- l[[nn]]
  hist(p.val[p.val < 1], br = 50, xlab = "p values",
       main = paste("Histogram for method:", nn, sep = " "))
}


###################################################
### chunk number 33: 
###################################################
par(mfrow = c(2, 2))
for(nn in names(l)) {
  p.val <- l[[nn]]
  hist(p.val[p.val < 1], br = 50, xlab = "p values",
       main = paste("Histogram for method:", nn, sep = " "))
}


###################################################
### chunk number 34: 
###################################################
showSigOfNodes(GOdata, score(resultFis), firstTerms = 5, useInfo = 'all')
showSigOfNodes(GOdata, score(resultWeight), firstTerms = 5, useInfo = 'def')


###################################################
### chunk number 35: 
###################################################
printGraph(GOdata, resultFis, firstSigNodes = 5, fn.prefix = "tGO", useInfo = "all", pdfSW = TRUE)
printGraph(GOdata, resultWeight, firstSigNodes = 5, fn.prefix = "tGO", useInfo = "def", pdfSW = TRUE)


###################################################
### chunk number 36: 
###################################################
printGraph(GOdata, resultWeight, firstSigNodes = 5, fn.prefix = "tGO", pdfSW = TRUE)


###################################################
### chunk number 37: 
###################################################
printGraph(GOdata, resultWeight, firstSigNodes = 10, resultFis, fn.prefix = "tGO", useInfo = "def")
printGraph(GOdata, resultElim, firstSigNodes = 15, resultFis, fn.prefix = "tGO", useInfo = "all")