### R code from vignette source 'EnrichmentBrowser.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: style
###################################################
BiocStyle::latex()


###################################################
### code chunk number 2: read-eset
###################################################
library(EnrichmentBrowser)
data.dir <- system.file("extdata", package="EnrichmentBrowser")
exprs.file <- file.path(data.dir, "exprs.tab")
pdat.file <- file.path(data.dir, "pData.tab")
fdat.file <- file.path(data.dir, "fData.tab")
eset <- read.eset(exprs.file, pdat.file, fdat.file)


###################################################
### code chunk number 3: help (eval = FALSE)
###################################################
## ?read.eset
## ?ExpressionSet


###################################################
### code chunk number 4: load-ALL
###################################################
library(ALL)
data(ALL)


###################################################
### code chunk number 5: subset-ALL
###################################################
ind.bs <- grep("^B", ALL$BT)
ind.mut <- which(ALL$mol.biol %in% c("BCR/ABL", "NEG"))
sset <- intersect(ind.bs, ind.mut)
all.eset <- ALL[, sset]


###################################################
### code chunk number 6: show-ALL
###################################################
dim(all.eset)
exprs(all.eset)[1:4,1:4]


###################################################
### code chunk number 7: probe2gene
###################################################
all.eset <- probe.2.gene.eset(all.eset)
head(featureNames(all.eset))


###################################################
### code chunk number 8: show-probe2gene
###################################################
head(fData(eset))


###################################################
### code chunk number 9: load-airway
###################################################
library(airway)
data(airway)


###################################################
### code chunk number 10: preproc-airway
###################################################
expr <- assays(airway)[[1]]
expr <- expr[grep("^ENSG", rownames(expr)),]
expr <- expr[rowMeans(expr) > 10,]
air.eset <- new("ExpressionSet", exprs=expr, annotation="hsa")
dim(air.eset)
exprs(air.eset)[1:4,1:4]


###################################################
### code chunk number 11: norm-ma
###################################################
before.norm <- exprs(all.eset)
all.eset <- normalize(all.eset, norm.method="quantile")
after.norm <- exprs(all.eset)


###################################################
### code chunk number 12: plot-norm
###################################################
par(mfrow=c(1,2))
boxplot(before.norm)
boxplot(after.norm)


###################################################
### code chunk number 13: norm-rseq
###################################################
norm.air <- normalize(air.eset, norm.method="quantile")


###################################################
### code chunk number 14: lgc (eval = FALSE)
###################################################
## ids <- head(featureNames(air.eset))
## lgc <- EDASeq::getGeneLengthAndGCContent(ids, org="hsa", mode="biomart") 


###################################################
### code chunk number 15: norm2-rseq
###################################################
lgc.file <- file.path(data.dir, "air_lgc.tab")
fData(air.eset) <- read.delim(lgc.file)
norm.air <- normalize(air.eset, within=TRUE)


###################################################
### code chunk number 16: sample-groups-ALL
###################################################
pData(all.eset)$GROUP <- ifelse(all.eset$mol.biol == "BCR/ABL", 1, 0)
table(pData(all.eset)$GROUP)


###################################################
### code chunk number 17: sample-groups-airway
###################################################
pData(air.eset)$GROUP <- ifelse(colData(airway)$dex == "trt", 1, 0)
table(pData(air.eset)$GROUP)


###################################################
### code chunk number 18: sample-blocks
###################################################
pData(air.eset)$BLOCK <- colData(airway)$cell
table(pData(air.eset)$BLOCK)


###################################################
### code chunk number 19: DE-ana-ALL
###################################################
all.eset <- de.ana(all.eset)
head(fData(all.eset), n=4)


###################################################
### code chunk number 20: plot-DE
###################################################
par(mfrow=c(1,2))
pdistr(fData(all.eset)$ADJ.PVAL)
volcano(fData(all.eset)$FC, fData(all.eset)$ADJ.PVAL)


###################################################
### code chunk number 21: DE-exmpl
###################################################
fData(all.eset)[ which.min(fData(all.eset)$ADJ.PVAL), ]


###################################################
### code chunk number 22: DE-ana-airway
###################################################
air.eset <- de.ana(air.eset, de.method="edgeR")
head(fData(air.eset), n=4)


###################################################
### code chunk number 23: idmap-idtypes
###################################################
id.types("hsa")


###################################################
### code chunk number 24: idmap-airway
###################################################
head(featureNames(air.eset))
air.eset <- map.ids(air.eset, from="ENSEMBL", to="ENTREZID")
head(featureNames(air.eset))


###################################################
### code chunk number 25: get-kegg-gs (eval = FALSE)
###################################################
## kegg.gs <- get.kegg.genesets("hsa")


###################################################
### code chunk number 26: get-go-gs (eval = FALSE)
###################################################
## go.gs <- get.go.genesets(org="hsa", onto="BP", mode="GO.db")


###################################################
### code chunk number 27: parseGMT
###################################################
gmt.file <- file.path(data.dir, "hsa_kegg_gs.gmt")
hsa.gs <- parse.genesets.from.GMT(gmt.file)
length(hsa.gs)
hsa.gs[1:2]


###################################################
### code chunk number 28: sbea-methods
###################################################
sbea.methods()


###################################################
### code chunk number 29: sbea
###################################################
sbea.res <- sbea(method="ora", eset=all.eset, gs=hsa.gs, perm=0, alpha=0.05)
gs.ranking(sbea.res)


###################################################
### code chunk number 30: ea-browse (eval = FALSE)
###################################################
## ea.browse(sbea.res)


###################################################
### code chunk number 31: dummy-sbea
###################################################
dummy.sbea <- function(eset, gs, alpha, perm)
{
        sig.ps <- sample(seq(0,0.05, length=1000),5)
        insig.ps <- sample(seq(0.1,1, length=1000), length(gs)-5)
        ps <- sample(c(sig.ps, insig.ps), length(gs))
        names(ps) <- names(gs)
        return(ps)
}


###################################################
### code chunk number 32: sbea2
###################################################
sbea.res2 <- sbea(method=dummy.sbea, eset=all.eset, gs=hsa.gs)
gs.ranking(sbea.res2)


###################################################
### code chunk number 33: dwnld-pwys (eval = FALSE)
###################################################
## pwys <- download.kegg.pathways("hsa")


###################################################
### code chunk number 34: compile-grn
###################################################
pwys <- file.path(data.dir, "hsa_kegg_pwys.zip")
hsa.grn <- compile.grn.from.kegg(pwys)
head(hsa.grn)


###################################################
### code chunk number 35: nbea-methods
###################################################
nbea.methods()


###################################################
### code chunk number 36: nbea
###################################################
nbea.res <- nbea(method="ggea", eset=all.eset, gs=hsa.gs, grn=hsa.grn)
gs.ranking(nbea.res)


###################################################
### code chunk number 37: ggea-graph
###################################################
par(mfrow=c(1,2))
ggea.graph(
    gs=hsa.gs[["hsa05217_Basal_cell_carcinoma"]], 
    grn=hsa.grn, eset=all.eset)
ggea.graph.legend()


###################################################
### code chunk number 38: combine
###################################################
res.list <- list(sbea.res, nbea.res)
comb.res <- comb.ea.results(res.list)


###################################################
### code chunk number 39: browse-comb (eval = FALSE)
###################################################
## ea.browse(comb.res, graph.view=hsa.grn, nr.show=5)


###################################################
### code chunk number 40: all-in-one (eval = FALSE)
###################################################
## ebrowser(   meth=c("ora", "ggea"), 
##         exprs=exprs.file, pdat=pdat.file, fdat=fdat.file, 
##         org="hsa", gs=hsa.gs, grn=hsa.grn, comb=TRUE, nr.show=5)