###################################################
### chunk number 1: setup
###################################################
options(width = 40)


###################################################
### chunk number 2: getting ALLfilt_bcrneg
###################################################
library(ALL)
library(hgu95av2.db)
data(ALL)
bcell <- grep("^B", as.character(ALL$BT))
types <- c("NEG", "BCR/ABL") 
moltyp <- which(as.character(ALL$mol.biol) %in% types)
# subsetting
ALL_bcrneg <- ALL[, intersect(bcell, moltyp)]
ALL_bcrneg$BT <- factor(ALL_bcrneg$BT)
ALL_bcrneg$mol.biol <- factor(ALL_bcrneg$mol.biol)
# nonspecific filter: remove genes that does not
## show much variation across samples
library(genefilter)  
filt_bcrneg <- nsFilter(ALL_bcrneg, 
                        var.cutoff=0.5)
ALLfilt_bcrneg <- filt_bcrneg$eset


###################################################
### chunk number 3: using KEGG
###################################################
library(KEGG.db)
library(GSEABase)
gsc <- GeneSetCollection(ALLfilt_bcrneg, 
                         setType=KEGGCollection())
Am <- incidence(gsc)


###################################################
### chunk number 4: nsF
###################################################
nsF <- ALLfilt_bcrneg[colnames(Am), ]


###################################################
### chunk number 5: exercise 1
###################################################
dim(nsF)
dim(Am)
nGene <- rowSums(Am)
rownames(Am)[nGene < 10]
sort(nGene, decreasing=TRUE)[1]
KEGGPATHID2NAME[["05200"]]


###################################################
### chunk number 6: compute per-gene test
###################################################
rtt <- rowttests(nsF, "mol.biol")
names(rtt)
rttStats <- rtt$statistic


###################################################
### chunk number 7: reduce am
###################################################
selectedRows <- (rowSums(Am) > 10)
Am2 <- Am[selectedRows, ]


###################################################
### chunk number 8: compute z score
###################################################
tA <- as.vector(Am2 %*% rttStats)
tAadj <- tA /sqrt(rowSums(Am2))
names(tAadj) <- rownames(Am2)


###################################################
### chunk number 9: exercise 2
###################################################
smPW <- tAadj[tAadj < -5]
mget(names(smPW),KEGGPATHID2NAME)
lgPW <- tAadj[tAadj > 5]
mget(names(lgPW), KEGGPATHID2NAME)


###################################################
### chunk number 10: KEGG2heatmap
###################################################
KEGG2heatmap("03010", nsF, "hgu95av2")


###################################################
### chunk number 11: permutation
###################################################
library(Category)
set.seed(123)
pvals <- gseattperm(nsF, nsF$mol.biol, Am2, 1000)
pvalCut <- 0.05
lowC <- rownames(pvals)[pvals[, 1] <= pvalCut]
unlist(getPathNames(lowC), use.names=FALSE)


###################################################
### chunk number 12: remove probes
###################################################
ALLfilt <- nsFilter(ALL_bcrneg, require.entez=TRUE,
                    remove.dupEntrez=TRUE,
                    require.CytoBand=TRUE,
                    var.func=IQR, 
                    var.cutoff=0.5)$eset


###################################################
### chunk number 13: t-stats for ALLfilt
###################################################
library(limma)
design <- model.matrix(~0 + ALLfilt$mol.biol)
colnames(design) <- c("BCR/ABL", "NEG")
contr <- c(1, -1)
fit1 <- lmFit(ALLfilt, design)
fit2 <- contrasts.fit(fit1, contr)
fit3 <- eBayes(fit2)
tlimma <- topTable(fit3, number=nrow(fit3), 
                   adjust.method="none")
## annotation
entrezUniverse <- unlist(mget(tlimma$ID, 
                              hgu95av2ENTREZID))
tstats <- tlimma$t
names(tstats) <- entrezUniverse


###################################################
### chunk number 14: Chrmaplinearmparams
###################################################
library(Category)
params <- new("ChrMapLinearMParams",
              conditional=FALSE,
              testDirection="up",
              universeGeneIds=entrezUniverse,
              geneStats=tstats,
              annotation="hgu95av2",
              pvalueCutoff=0.01,
              minSize=4L)    


###################################################
### chunk number 15: linearMTest
###################################################
lman <- linearMTest(params)
lman
summary(lman)


###################################################
### chunk number 16: conditional test
###################################################
slotNames(params)
paramsCond <- params
paramsCond@conditional <- TRUE
lmanCond <- linearMTest(paramsCond)
summary(lmanCond)


###################################################
### chunk number 17: gene selection by t-test
###################################################
ttests <- rowttests(ALLfilt_bcrneg, "mol.biol")
smPV <- ttests[ttests$p.value < 0.005, ]
selectedEntrezIds <- unlist(mget(rownames(smPV),  
                                 hgu95av2ENTREZID))
entrezUniverse=unlist(mget(featureNames(ALLfilt_bcrneg),  
  hgu95av2ENTREZID))

#library(limma)
#design <- model.matrix(~0 + ALLfilt_bcrneg$mol.biol)
#colnames(design) <- c("BCR/ABL", "NEG")
#contr <- c(1, -1)
#fit1 <- lmFit(ALLfilt_bcrneg, design)
#fit2 <- contrasts.fit(fit1, contr)
#fit3 <- eBayes(fit2)
#tlimma <- topTable(fit3, number=nrow(fit3), 
#                   adjust.method="none")
#smPV <- tlimma[tlimma$P.Value < 0.01, ]
#dim(smPV)
#selectedEntrezIds <- unlist(mget(smPV$ID, hgu95av2ENTREZID))
#entrezUniverse=unlist(mget(featureNames(ALLfilt_bcrneg),  hgu95av2ENTREZID))


###################################################
### chunk number 18: GOHyperGParams
###################################################
library(GOstats)
hgCutoff <- 0.001
GOparams <- new("GOHyperGParams", 
                geneIds=selectedEntrezIds, 
                universeGeneIds=entrezUniverse, 
                annotation="hgu95av2.db", 
                ontology="BP", 
                pvalueCutoff=0.001, 
                conditional=TRUE, 
                testDirection="over")


###################################################
### chunk number 19: hyperGTest eval=FALSE
###################################################
## hgOver <- hyperGTest(GOparams)
## class(hgOver)
## summary(hgOver)


###################################################
### chunk number 20: use htmlReport to generate an html report eval=FALSE
###################################################
## showMethods("htmlReport")
## htmlReport(hgOver, file="hgResult.html")