gage.R
In gage: Generally Applicable Gene-set Enrichment for Pathway Analysis

### R code from vignette source 'gage.Rnw'

###################################################
### code chunk number 1: synopsis1 (eval = FALSE)
###################################################
## library(gage)
## data(gse16873)
## hn=(1:6)*2-1
## dcis=(1:6)*2


###################################################
### code chunk number 2: synopsis2 (eval = FALSE)
###################################################
## data(kegg.gs)
## gse16873.kegg.p <- gage(gse16873, gsets = kegg.gs, ref = hn, samp = dcis)


###################################################
### code chunk number 3: synopsis3 (eval = FALSE)
###################################################
## library(gageData)
## data(go.sets.hs)
## data(go.subs.hs)
## gse16873.bp.p <- gage(gse16873, gsets = go.sets.hs[go.subs.hs$BP], ref = hn, samp = dcis)
## gse16873.mf.p <- gage(gse16873, gsets = go.sets.hs[go.subs.hs$MF], ref = hn, samp = dcis)
## gse16873.cc.p <- gage(gse16873, gsets = go.sets.hs[go.subs.hs$CC], ref = hn, samp = dcis)


###################################################
### code chunk number 4: install (eval = FALSE)
###################################################
## if (!requireNamespace("BiocManager", quietly=TRUE))
##     install.packages("BiocManager")
## BiocManager::install(c("gage","gageData"))


###################################################
### code chunk number 5: install (eval = FALSE)
###################################################
## install.packages(c("/your/local/directory/gage_2.9.1.tar.gz",
##     "/your/local/directory/gageData_1.0.0.tar.gz"), 
##     repos = NULL, type = "source")


###################################################
### code chunk number 6: start
###################################################
options(width=80)


###################################################
### code chunk number 7: start
###################################################
library(gage)


###################################################
### code chunk number 8: start (eval = FALSE)
###################################################
## library(help=gage)


###################################################
### code chunk number 9: start (eval = FALSE)
###################################################
## help(gage)
## ?gage


###################################################
### code chunk number 10: dataPrep
###################################################
data(gse16873)
cn=colnames(gse16873)
hn=grep('HN',cn, ignore.case =T)
adh=grep('ADH',cn, ignore.case =T)
dcis=grep('DCIS',cn, ignore.case =T)
print(hn)
print(dcis)


###################################################
### code chunk number 11: dataPrep
###################################################
data(kegg.gs)
data(go.gs)
lapply(kegg.gs[1:3],head)
head(rownames(gse16873))


###################################################
### code chunk number 12: kegg.gsets
###################################################
kg.hsa=kegg.gsets("hsa")
kegg.gs=kg.hsa$kg.sets[kg.hsa$sigmet.idx]
#save(kegg.gs, file="kegg.hsa.sigmet.gsets.RData")


###################################################
### code chunk number 13: korg
###################################################
#kegg.gsets works with 3000 KEGG species,for examples:
data(korg)
head(korg[,1:3])


###################################################
### code chunk number 14: go.gsets
###################################################
go.hs=go.gsets(species="human")
go.bp=go.hs$go.sets[go.hs$go.subs$BP]
go.mf=go.hs$go.sets[go.hs$go.subs$MF]
go.cc=go.hs$go.sets[go.hs$go.subs$CC]
#save(go.bp, go.mf, go.cc, file="go.hs.gsets.RData")


###################################################
### code chunk number 15: bods
###################################################
#for Bioconductor species supported by go.gsets function:
data(bods)
print(bods)


###################################################
### code chunk number 16: gage.1direction
###################################################
gse16873.kegg.p <- gage(gse16873, gsets = kegg.gs, 
    ref = hn, samp = dcis)
#go.gs here only the first 1000 entries as a fast example.
gse16873.go.p <- gage(gse16873, gsets = go.gs, 
    ref = hn, samp = dcis)
#or we can do analysis on BP, MF, or CC subcategories if we've
#generated the data as above.
#gse16873.bp.p <- gage(gse16873, gsets = go.bp,
#                        ref = hn, samp = dcis)
str(gse16873.kegg.p, strict.width='wrap')
head(gse16873.kegg.p$greater[, 1:5], 4)
head(gse16873.kegg.p$less[, 1:5], 4)
head(gse16873.kegg.p$stats[, 1:5], 4)


###################################################
### code chunk number 17: gage.2direction
###################################################
gse16873.kegg.2d.p <- gage(gse16873, gsets = kegg.gs, 
    ref = hn, samp = dcis, same.dir = F)
head(gse16873.kegg.2d.p$greater[,1:5], 4)
head(gse16873.kegg.2d.p$stats[,1:5], 4)


###################################################
### code chunk number 18: page
###################################################
gse16873.kegg.page.p <- gage(gse16873, gsets = kegg.gs, 
    ref = hn, samp = dcis, saaTest = gs.zTest)
head(gse16873.kegg.page.p$greater[, 1:5], 4)


###################################################
### code chunk number 19: gage.alternative (eval = FALSE)
###################################################
## gse16873.kegg.t.p <- gage(gse16873, gsets = kegg.gs, 
##     ref = hn, samp = dcis, use.fold = F)


###################################################
### code chunk number 20: gage.alternative (eval = FALSE)
###################################################
## gse16873.kegg.rk.p <- gage(gse16873, gsets = kegg.gs, 
##     ref = hn, samp = dcis, rank.test = T)


###################################################
### code chunk number 21: gage.alternative (eval = FALSE)
###################################################
## gse16873.kegg.ks.p <- gage(gse16873, gsets = kegg.gs, 
##     ref = hn, samp = dcis, saaTest = gs.KSTest)


###################################################
### code chunk number 22: gage.alternative (eval = FALSE)
###################################################
## gse16873.kegg.unpaired.p <- gage(gse16873, gsets = kegg.gs, 
##     ref = hn, samp = dcis, compare = "unpaired")


###################################################
### code chunk number 23: gage.alternative (eval = FALSE)
###################################################
## gse16873.kegg.gamma.p <- gage(gse16873, gsets = kegg.gs, 
##     ref = hn, samp = dcis, use.stouffer=FALSE)


###################################################
### code chunk number 24: symbol.ID
###################################################
head(rownames(gse16873))
gse16873.sym<-gse16873
data(egSymb)
rownames(gse16873.sym)<-eg2sym(rownames(gse16873.sym))
head(rownames(gse16873.sym))


###################################################
### code chunk number 25: symbol.ID
###################################################
kegg.gs.sym<-lapply(kegg.gs, eg2sym)
lapply(kegg.gs.sym[1:3],head)
gse16873.kegg.p2 <- gage(gse16873.sym, gsets = kegg.gs.sym, 
    ref = hn, samp = dcis)


###################################################
### code chunk number 26: full.table
###################################################
write.table(gse16873.kegg.2d.p$greater, file = "gse16873.kegg.2d.p.txt", 
    sep = "\t")
write.table(rbind(gse16873.kegg.p$greater, gse16873.kegg.p$less), 
    file = "gse16873.kegg.p.txt", sep = "\t")


###################################################
### code chunk number 27: significant.genesets
###################################################
gse16873.kegg.sig<-sigGeneSet(gse16873.kegg.p, outname="gse16873.kegg")
str(gse16873.kegg.sig, strict.width='wrap')
gse16873.kegg.2d.sig<-sigGeneSet(gse16873.kegg.2d.p, outname="gse16873.kegg")
str(gse16873.kegg.2d.sig, strict.width='wrap')
write.table(gse16873.kegg.2d.sig$greater, file = "gse16873.kegg.2d.sig.txt", 
    sep = "\t")
write.table(rbind(gse16873.kegg.sig$greater, gse16873.kegg.sig$less), 
    file = "gse16873.kegg.sig.txt", sep = "\t")


###################################################
### code chunk number 28: nonredundant.genesets
###################################################
gse16873.kegg.esg.up <- esset.grp(gse16873.kegg.p$greater, 
    gse16873, gsets = kegg.gs, ref = hn, samp = dcis, 
    test4up = T, output = T, outname = "gse16873.kegg.up", make.plot = F)
gse16873.kegg.esg.dn <- esset.grp(gse16873.kegg.p$less, 
    gse16873, gsets = kegg.gs, ref = hn, samp = dcis, 
    test4up = F, output = T, outname = "gse16873.kegg.dn", make.plot = F)
names(gse16873.kegg.esg.up)
head(gse16873.kegg.esg.up$essentialSets, 4)
head(gse16873.kegg.esg.up$setGroups, 4)
head(gse16873.kegg.esg.up$coreGeneSets, 4)


###################################################
### code chunk number 29: sel.gene.expression
###################################################
rownames(gse16873.kegg.p$greater)[1:3]
gs=unique(unlist(kegg.gs[rownames(gse16873.kegg.p$greater)[1:3]]))
essData=essGene(gs, gse16873, ref =hn, samp =dcis)
head(essData, 4)
ref1=1:6
samp1=7:12
for (gs in rownames(gse16873.kegg.p$greater)[1:3]) {
    outname = gsub(" |:|/", "_", substr(gs, 10, 100))
    geneData(genes = kegg.gs[[gs]], exprs = essData, ref = ref1, 
        samp = samp1, outname = outname, txt = T, heatmap = T, 
        Colv = F, Rowv = F, dendrogram = "none", limit = 3, scatterplot = T)
}


###################################################
### code chunk number 30: all.gene.expression
###################################################
for (gs in rownames(gse16873.kegg.p$greater)[1:3]) {
    outname = gsub(" |:|/", "_", substr(gs, 10, 100))
    outname = paste(outname, "all", sep=".")
    geneData(genes = kegg.gs[[gs]], exprs = gse16873, ref = hn, 
        samp = dcis, outname = outname, txt = T, heatmap = T, 
        Colv = F, Rowv = F, dendrogram = "none", limit = 3, scatterplot = T)
}


###################################################
### code chunk number 31: pathview
###################################################
library(pathview)
gse16873.d <- gse16873[ ,dcis] - gse16873[ ,hn]
path.ids=c("hsa04110 Cell cycle", "hsa00020 Citrate cycle (TCA cycle)")
path.ids2 <- substr(path.ids, 1, 8)
#native KEGG view
pv.out.list <- sapply(path.ids2, function(pid) pathview(gene.data = gse16873.d[, 
                        1:2], pathway.id = pid, species = "hsa"))
#Graphviz view
pv.out.list <- sapply(path.ids2, function(pid) pathview(gene.data = gse16873.d[, 
                        1:2], pathway.id = pid, species = "hsa", kegg.native=F,
                        sign.pos="bottomleft")) 


###################################################
### code chunk number 32: pathview.all (eval = FALSE)
###################################################
## sel <- gse16873.kegg.p$greater[, "q.val"] < 0.1 & !is.na(gse16873.kegg.p$greater[,
## "q.val"])
## path.ids <- rownames(gse16873.kegg.p$greater)[sel]
## path.ids2 <- substr(path.ids, 1, 8)
## pv.out.list <- sapply(path.ids2, function(pid) pathview(gene.data = gse16873.d[,
##                         1:2], pathway.id = pid, species = "hsa"))


###################################################
### code chunk number 33: gagePipe
###################################################
#introduce another half dataset
library(gageData)
data(gse16873.2)
cn2=colnames(gse16873.2)
hn2=grep('HN',cn2, ignore.case =T)
dcis2=grep('DCIS',cn2, ignore.case =T)
#batch GAGE analysis on the combined dataset
gse16873=cbind(gse16873, gse16873.2)
dataname='gse16873' #output data prefix
sampnames=c('dcis.1', 'dcis.2')
refList=list(hn, hn2+12)
sampList=list(dcis, dcis2+12)
gagePipe(gse16873, gsname = "kegg.gs", dataname = "gse16873", 
    sampnames = sampnames, ref.list = refList, samp.list = sampList, 
    comp.list = "paired")


###################################################
### code chunk number 34: gageComp
###################################################
load('gse16873.gage.RData')
gageComp(sampnames, dataname, gsname = "kegg.gs", 
    do.plot = TRUE)


###################################################
### code chunk number 35: heter.gage
###################################################
boxplot(data.frame(gse16873))
gse16873.kegg.heter.p <- heter.gage(gse16873, gsets = kegg.gs, 
    ref.list = refList, samp.list = sampList)
gse16873.kegg.heter.2d.p <- heter.gage(gse16873, gsets = kegg.gs, 
    ref.list = refList, samp.list = sampList, same.dir = F)