Some gsea

###################################################
### chunk: ALL
###################################################
if(!exists("biocLite")) {
source("http://www.bioconductor.org/biocLite.R")
}

is.ALL.found=require("ALL")
if(!is.ALL.found) {
biocLite("ALL")
}

is.genefilter.found = require("genefilter")
if(!is.genefilter.found) {
biocLite("genefilter")
}

is.KEGGdb.found = require("KEGG.db")
if(!is.KEGGdb.found) {
biocLite("KEGG.db")
}

library("ALL")
data("ALL")


###################################################
### chunk: bcrabl
###################################################
bcell = grep("^B", as.character(ALL$BT))
moltyp = which(as.character(ALL$mol.biol)
%in% c("NEG", "BCR/ABL"))
ALL_bcrneg = ALL[, intersect(bcell, moltyp)]
ALL_bcrneg$mol.biol = factor(ALL_bcrneg$mol.biol)


###################################################
### chunk: nsFilter
###################################################
library("genefilter")
ALLfilt_bcrneg = nsFilter(ALL_bcrneg, var.cutoff=0.5)$eset


###################################################
### chunk: KEGGimat
###################################################
library("GSEABase")
gsc = GeneSetCollection(ALLfilt_bcrneg,
setType=KEGGCollection())
Am = incidence(gsc)
dim(Am)


###################################################
### chunk: EsetfromKEGG
###################################################
nsF = ALLfilt_bcrneg[colnames(Am),]


###################################################
### chunk: compttests
###################################################
rtt = rowttests(nsF, "mol.biol")
rttStat = rtt$statistic


###################################################
### chunk: reducetoInt
###################################################
selectedRows = (rowSums(Am)>10)
Am2 = Am[selectedRows, ]


###################################################
### chunk: pctests
###################################################
tA = as.vector(Am2 %*% rttStat)
tAadj = tA/sqrt(rowSums(Am2))
names(tA) = names(tAadj) = rownames(Am2)


###################################################
### chunk: qqplot
###################################################
x11()
qqpos=qqnorm(tAadj)
o = order(tAadj, decreasing=FALSE)
scan()
xpos = qqpos$x[o[1]]
ypos = qqpos$y[o[1]]
points(x = xpos, y = ypos, col = "red",
pch = 19, cex = 1.5)
legend(legend=KEGGPATHID2NAME[[names(tAadj)[o[1]]]],
x = xpos,
y = ypos + 1,
text.col = "red",
bty="n", cex = 1.5)
scan()

###################################################
### chunk: findSmPW
###################################################
library("KEGG.db")
smPW = tAadj[tAadj == min(tAadj, na.rm=T)]
pwName = KEGGPATHID2NAME[[names(smPW)]]
pwName


###################################################
### chunk: mnplot
###################################################
x11()
KEGGmnplot(names(smPW), nsF, "hgu95av2", nsF$"mol.biol",
pch=16, col="darkblue")
scan()

###################################################
### chunk: heatmap
###################################################
x11()
sel = as.integer(nsF$mol.biol)
KEGG2heatmap(names(smPW), nsF, "hgu95av2",
col=colorRampPalette(c("white", "darkblue"))(256),
ColSideColors=c("black", "white")[sel])
scan()