## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)

## ----eval=FALSE---------------------------------------------------------------
# ?LimROTS

## ----echo=FALSE, message=FALSE, results="hide"--------------------------------
library(LimROTS, quietly = TRUE)
data("UPS1.Case4")

## ----echo=FALSE---------------------------------------------------------------
table(UPS1.Case4$fake.batch, UPS1.Case4$Conc.)

## ----load-packages, message=FALSE---------------------------------------------
library(LimROTS, quietly = TRUE)
library(BiocParallel, quietly = TRUE)
library(ggplot2, quietly = TRUE)
library(SummarizedExperiment, quietly = TRUE)

## ----load-dataset-------------------------------------------------------------
data("UPS1.Case4")
print(UPS1.Case4)

## ----set-seed-----------------------------------------------------------------
set.seed(1234, kind = "default", sample.kind = "default")

## ----run-limrots--------------------------------------------------------------
# Define analysis parameters
meta.info   <- c("Conc.", "tool", "fake.batch")
niter       <- 100   # bootstrap iterations (use >= 1000 for real analyses)
K           <- 100   # top list size for reproducibility optimization
formula.str <- "~ 0 + Conc. + tool + fake.batch"

# Run LimROTS
UPS1.Case4 <- LimROTS(
    x           = UPS1.Case4,
    niter       = niter,
    K           = K,
    meta.info   = meta.info,
    group       = "Conc.",
    formula.str = formula.str,
    trend       = TRUE,
    robust      = TRUE,
    permutating.group = FALSE
)

## ----bpparam-example----------------------------------------------------------
# Windows
# BPPARAM <- SnowParam(workers = 4)

# Linux / macOS
# BPPARAM <- MulticoreParam(workers = 4)

# Pass to LimROTS via: LimROTS(..., BPPARAM = BPPARAM)

## ----volcano-plot-------------------------------------------------------------
# Extract results
result_df <- data.frame(
    rowData(UPS1.Case4),
    row.names = rownames(UPS1.Case4)
)

# Annotate proteins
result_df$label <- ifelse(
    grepl("HUMAN", result_df$GeneID),
    "UPS1 (true positive)",
    "E. coli (true negative)"
)

ggplot(result_df, aes(
    x     = corrected.logfc,
    y     = -log10(qvalue),
    color = label
)) +
    geom_point(alpha = 0.7, size = 1.2) +
    geom_hline(
        yintercept = -log10(0.05),
        linetype   = "dashed",
        color      = "steelblue"
    ) +
    scale_color_manual(values = c(
        "E. coli (true negative)" = "grey60",
        "UPS1 (true positive)"    = "firebrick"
    )) +
    labs(
        title = "LimROTS \u2014 UPS1 Case Study",
        x     = "Corrected Log Fold Change",
        y     = expression(-log[10](q~value)),
        color = NULL
    ) +
    theme_bw(base_size = 12) +
    theme(legend.position = "top")

## ----quality-control, results="hide", message=FALSE, warning=FALSE------------
plot(metadata(UPS1.Case4)[["q_values"]])
hist(metadata(UPS1.Case4)[["q_values"]])
print(summary(metadata(UPS1.Case4)[["q_values"]]))

## ----session-info-------------------------------------------------------------
sessionInfo()