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

## -----------------------------------------------------------------------------
make_mock_dataset <- function() {
    set.seed(100)
    mu <- 2^runif(1000, 3, 6)
    grouping <- rep(c("A", "B", "ctrl"), each=3)
    counts <- matrix(rnbinom(length(mu)* length(grouping), mu=mu, size=20), ncol=length(grouping))
    rownames(counts) <- sprintf("BARCODE-%s", seq_along(mu))

    library(SummarizedExperiment)
    se <- SummarizedExperiment(list(counts=counts), colData=DataFrame(group=grouping))
    rowData(se)$type <- "other"
    rowData(se)$gene <- paste0("GENE-", sample(LETTERS, length(mu), replace=TRUE))
    rowData(se)$type[1:50] <- "NEG"
    rowData(se)$gene[1:50] <- "NEG"

    se
}

se <- make_mock_dataset()
se

## ----fig.show="hide"----------------------------------------------------------
library(augere.screen)
output.dir <- tempfile()
res <- runScreen(
    se,
    group="group",
    comparison=c("A", "B"),
    gene.field="gene",
    filter.reference.factor="group",
    filter.reference.levels="ctrl",
    norm.control.type.field="type",
    norm.control.types="NEG",
    output.dir=output.dir
)

# List of tables of barcode-level differential results.
names(res$barcode)
res$barcode[[1]]

# List of tables of gene-level statistics.
names(res$gene)

# Copy of 'se' with normalized abundance values.
res$normalized

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

## -----------------------------------------------------------------------------
reloaded <- readResult(file.path(output.dir, "results", "barcode-1"))
reloaded$x # the result itself
str(reloaded$metadata) # along with some metadata

## -----------------------------------------------------------------------------
names(res$gene$simes)
res$gene$simes[[1]]

## -----------------------------------------------------------------------------
res$gene$`holm-min`[[1]]

## -----------------------------------------------------------------------------
res$gene$fry[[1]]

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

## -----------------------------------------------------------------------------
output.dry.dir <- tempfile()
runScreen(
    wrapInput(se, commands=as.character(body(make_mock_dataset))[-1]),
    group="group",
    comparison=c("A", "B"),
    gene.field="gene",
    filter.reference.factor="group",
    filter.reference.levels="ctrl",
    norm.control.type.field="type",
    norm.control.types="NEG",
    output.dir=output.dry.dir,
    dry.run=TRUE
)

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

## ----fig.show="hide"----------------------------------------------------------
output.custom.dir <- tempfile()
res <- runScreen(
    se,
    group="group",
    comparison=c("A", "B"),
    gene.field="gene",
    filter.reference.factor="group",
    filter.reference.levels="ctrl",
    norm.control.type.field="type",
    norm.control.types="NEG",
    row.data=c("gene", "type"),
    gene.data="type",
    metadata=list(
        project="my_phd_project",
        collaborators=list("anne", "bob", "charlie"),
        data_source="https://foo.bar.com"
    ),
    output.dir=output.custom.dir
)

# rowData columns are now included in the barcode- and gene-level DFss.
res$barcode[[1]]
res$gene$simes[[1]]

# The custom metadata is also saved to disk.
str(readResult(file.path(output.custom.dir, "results", "barcode-1"))$metadata)

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