## ----style-knitr, eval=TRUE, echo=FALSE, results="asis"--------------------
BiocStyle::latex()

## ----setup, echo=FALSE-----------------------------------------------------
suppressMessages(suppressWarnings(suppressPackageStartupMessages({
library(EnrichmentBrowser)
library(ALL)
library(airway)
library(edgeR)
library(limma)
})))

## ----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, "colData.tab")
fdat.file <- file.path(data.dir, "rowData.tab")
eset <- read.eset(exprs.file, pdat.file, fdat.file)

## ----help, eval=FALSE------------------------------------------------------
#  ?read.eset
#  ?SummarizedExperiment

## ----sexp2eset-------------------------------------------------------------
eset <- as(eset, "ExpressionSet")

## ----eset2sexp-------------------------------------------------------------
eset <- as(eset, "RangedSummarizedExperiment")

## ----load-ALL--------------------------------------------------------------
library(ALL)
data(ALL)

## ----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]

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

## ----probe2gene------------------------------------------------------------
all.eset <- probe.2.gene.eset(all.eset)
head(rownames(all.eset))

## ----show-probe2gene-------------------------------------------------------
rowData(eset, use.names=TRUE)

## ----load-airway-----------------------------------------------------------
library(airway)
data(airway)

## ----preproc-airway--------------------------------------------------------
air.eset <- airway[grep("^ENSG", rownames(airway)),]
air.eset <- air.eset[rowMeans(assay(air.eset)) > 10,]
dim(air.eset)
assay(air.eset)[1:4,1:4]

## ----norm-ma---------------------------------------------------------------
before.norm <- assay(all.eset)
all.eset <- normalize(all.eset, norm.method="quantile")
after.norm <- assay(all.eset)

## ----plot-norm, fig.width=12, fig.height=6---------------------------------
par(mfrow=c(1,2))
boxplot(before.norm)
boxplot(after.norm)

## ----norm-rseq-------------------------------------------------------------
norm.air <- normalize(air.eset, norm.method="quantile")

## ----lgc, eval=FALSE-------------------------------------------------------
#  ids <- rownames(air.eset)
#  lgc <- EDASeq::getGeneLengthAndGCContent(ids, org="hsa", mode="biomart")

## ----norm2-rseq------------------------------------------------------------
lgc.file <- file.path(data.dir, "air_lgc.tab")
rowData(air.eset) <- read.delim(lgc.file)
norm.air <- normalize(air.eset, within=TRUE)

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

## ----sample-groups-airway--------------------------------------------------
air.eset$GROUP <- ifelse(airway$dex == "trt", 1, 0)
table(air.eset$GROUP)

## ----sample-blocks---------------------------------------------------------
air.eset$BLOCK <- airway$cell
table(air.eset$BLOCK)

## ----DE-ana-ALL------------------------------------------------------------
all.eset <- de.ana(all.eset)
rowData(all.eset, use.names=TRUE)

## ----plot-DE, fig.width=12, fig.height=6-----------------------------------
par(mfrow=c(1,2))
pdistr(rowData(all.eset)$ADJ.PVAL)
volcano(rowData(all.eset)$FC, rowData(all.eset)$ADJ.PVAL)

## ----DE-exmpl--------------------------------------------------------------
ind.min <- which.min( rowData(all.eset)$ADJ.PVAL )
rowData(all.eset, use.names=TRUE)[ ind.min, ]

## ----DE-ana-airway---------------------------------------------------------
air.eset <- de.ana(air.eset, de.method="edgeR")
rowData(air.eset, use.names=TRUE)

## ----idmap-idtypes---------------------------------------------------------
id.types("hsa")

## ----idmap-airway----------------------------------------------------------
head(rownames(air.eset))
air.eset <- map.ids(air.eset, org="hsa", from="ENSEMBL", to="ENTREZID")
head(rownames(air.eset))

## ----get-kegg-gs, eval=FALSE-----------------------------------------------
#  kegg.gs <- get.kegg.genesets("hsa")

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

## ----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]

## ----sbea-methods----------------------------------------------------------
sbea.methods()

## ----sbea------------------------------------------------------------------
sbea.res <- sbea(method="ora", eset=all.eset, gs=hsa.gs, perm=0, alpha=0.1)
gs.ranking(sbea.res)

## ----ea-browse, eval=FALSE-------------------------------------------------
#  ea.browse(sbea.res)

## ----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)
}

## ----sbea2-----------------------------------------------------------------
sbea.res2 <- sbea(method=dummy.sbea, eset=all.eset, gs=hsa.gs)
gs.ranking(sbea.res2)

## ----dwnld-pwys, eval=FALSE------------------------------------------------
#  pwys <- download.kegg.pathways("hsa")

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

## ----nbea-methods----------------------------------------------------------
nbea.methods()

## ----nbea------------------------------------------------------------------
nbea.res <- nbea(method="ggea", eset=all.eset, gs=hsa.gs, grn=hsa.grn)
gs.ranking(nbea.res)

## ----ggea-graph, fig.width=12, fig.height=6--------------------------------
par(mfrow=c(1,2))
ggea.graph(
    gs=hsa.gs[["hsa05217_Basal_cell_carcinoma"]], 
    grn=hsa.grn, eset=all.eset)
ggea.graph.legend()

## ----combine---------------------------------------------------------------
res.list <- list(sbea.res, nbea.res)
comb.res <- comb.ea.results(res.list)

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

## ----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)

## ----config-set------------------------------------------------------------
config.ebrowser(key="OUTDIR.DEFAULT", value="/my/out/dir")

## ----config-get------------------------------------------------------------
config.ebrowser("OUTDIR.DEFAULT")

## ----config-man, eval=FALSE------------------------------------------------
#  ?config.ebrowser

## ----deTbl-----------------------------------------------------------------
deTable <-
     matrix(c(28, 142, 501, 12000),
            nrow = 2,
            dimnames = list(c("DE", "Not.DE"),
                            c("In.gene.set", "Not.in.gene.set")))
deTable

## ----fisher----------------------------------------------------------------
fisher.test(deTable, alternative = "greater")