## ----echo=FALSE, results="hide"-----------------------------------------------
knitr::opts_chunk$set(error=FALSE, message=FALSE, warning=FALSE)
self <- BiocStyle::Biocpkg("augere.gsea")

## ----fig.show="hide"----------------------------------------------------------
library(augere.de)
se <- loadExampleDataset()

# Testing for DE between 'trt' and 'untrt' in the 'dex' grouping factor.
output.de <- tempfile()
res.de <- runVoom(se, group="dex", comparisons=c("trt", "untrt"), output.dir=output.de)
res.de$results[[1]]

## -----------------------------------------------------------------------------
library(org.Hs.eg.db)
sets <- select(org.Hs.eg.db, keys=keys(org.Hs.eg.db, "GO"), keytype="GO", columns="ENSEMBL")
sets <- split(sets$ENSEMBL, sets$GO)
sets <- head(sets, 100) # to cut down on vignette runtime.
length(sets)

## ----fig.show="hide"----------------------------------------------------------
library(augere.gsea)
output.gsea <- tempfile()

# Set the seed just in case any of the methods require randomization.
# runPrecomputed() also hard-codes some set.seed() calls in the generated
# Rmarkdown report to ensure that the report's results are reproducible.
set.seed(100000)

res.gsea <- runPrecomputed(
    res.de$results[[1]],
    sets,
    signif.field="FDR",
    signif.threshold=0.05,
    rank.field="t",
    output.dir=output.gsea
) 
res.gsea

## -----------------------------------------------------------------------------
fname <- file.path(output.gsea, "report.Rmd")
cat(head(readLines(fname), 50), sep="\n")

## -----------------------------------------------------------------------------
reloaded <- readResult(file.path(output.gsea, "results", "hypergeometric"))

# The result itself:
reloaded$x

# Along with some metadata:
str(reloaded$metadata)

## ----fig.show="hide"----------------------------------------------------------
output.up <- tempfile()
set.seed(100001)
res.up <- runPrecomputed(
    res.de$results[[1]],
    sets,
    alternative="up",
    signif.field="FDR",
    signif.threshold=0.05,
    rank.field="t",
    sign.field="LogFC", # supplying the sign for each gene. 
    output.dir=output.up,
    save.results=FALSE # skipping the save, to reduce run-time.
)
res.up$hypergeometric

## ----fig.show="hide"----------------------------------------------------------
output.either <- tempfile()
set.seed(100002)
res.either <- runPrecomputed(
    res.de$results[[1]],
    sets,
    alternative="either",
    signif.field="FDR",
    signif.threshold=0.05,
    rank.field="t",
    sign.field="LogFC",
    output.dir=output.either,
    save.results=FALSE
)
res.either$hypergeometric

## ----fig.show="hide"----------------------------------------------------------
output.de2 <- tempfile()
res.de2 <- runEdgeR(se, group="dex", comparisons=c("trt", "untrt"), output.dir=output.de2)
res.de2$results[[1]]

output.up2 <- tempfile()
set.seed(100003)
res.up2 <- runPrecomputed(
    res.de2$results[[1]],
    sets,
    signif.field="FDR",
    signif.threshold=0.05,
    rank.field="F",
    rank.sqrt=TRUE, # square-rooting the F statistic to mimic a t-test.
    sign.field="LogFC", # using this to restore the sign to the test statistic.
    alternative="up",
    output.dir=output.up2,
    save.results=FALSE,
)
res.up2$hypergeometric

## ----fig.show="hide"----------------------------------------------------------
output.dgst <- tempfile()
set.seed(200000)
res.dgst <- runContrast(
    se,
    sets,
    group="dex",
    comparison=c("trt", "untrt"),
    save.results=FALSE,
    output.dir=output.dgst
) 
res.dgst

## ----fig.show="hide"----------------------------------------------------------
output.dgst2 <- tempfile()
set.seed(200001)
res.dgst2 <- runContrast(
    se,
    sets,
    design=~dex,
    contrast="dexuntrt",
    output.dir=output.dgst2
) 

## -----------------------------------------------------------------------------
output.dry <- tempfile()
set.seed(300000)
runPrecomputed(
    res.de$results[[1]],
    sets,
    signif.field="FDR",
    signif.threshold=0.05,
    rank.field="t",
    output.dir=output.dry,
    dry.run=TRUE
)
list.files(output.dry)

## ----echo=FALSE---------------------------------------------------------------
fname <- file.path(output.dry, "report.Rmd")
lines <- readLines(fname)
has.stop <- grep("stop(", lines, fixed=TRUE)[1]
cat(lines[has.stop + (-5):5], sep="\n")

## ----fig.show="hide"----------------------------------------------------------
set.seed(300001)
runPrecomputed(
    wrapInput(res.de$results[[1]], commands=sprintf(
        # Pulling the result out of the directory saved by runVoom().
        "augere.core::readResult(%s)$x", deparse(file.path(output.de, "results", "de-1"))
    )),
    wrapInput(sets, commands="
library(org.Hs.eg.db)
sets <- select(org.Hs.eg.db, keys=keys(org.Hs.eg.db, 'GO'), keytype='GO', columns='ENSEMBL')
sets <- split(sets$ENSEMBL, sets$GO)
head(sets, 100)"), 
    signif.field="FDR",
    signif.threshold=0.05,
    rank.field="t",
    output.dir=output.dry,
    dry.run=TRUE
)

## ----echo=FALSE---------------------------------------------------------------
fname <- file.path(output.dry, "report.Rmd")
lines <- readLines(fname)
has.loader <- grep("readResult", lines, fixed=TRUE)[1]
cat(lines[has.loader + (-5):5], sep="\n")

## ----fig.show="hide"----------------------------------------------------------
library(GO.db)
annotated.sets <- List(sets)
mcols(annotated.sets)$name <- mapIds(
    GO.db,
    keys=names(annotated.sets),
    keytype="GOID",
    column="TERM"
)

output.anno.dir <- tempfile()
set.seed(300002)
res.anno <- runPrecomputed(
    res.de$results[[1]],
    annotated.sets,
    signif.field="FDR",
    signif.threshold=0.05,
    rank.field="t",
    save.results=FALSE,
    annotation="name",
    output.dir=output.anno.dir
)

# The desired columns of the 'mcols' are now included: 
res.anno$hypergeometric

## -----------------------------------------------------------------------------
sessionInfo()