## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
    collapse  = TRUE,
    comment   = "#>",
    fig.align = "center",
    fig.width = 7,
    fig.height = 4.5,
    message   = FALSE
)
library(levi)
library(ggplot2)
extdata <- function(...) system.file("extdata", ..., package = "levi")

## ----calibration-files--------------------------------------------------------
val <- extdata("validation")
list.files(val)

## ----type1--------------------------------------------------------------------
type1 <- read.csv(file.path(val, "type1_error.csv"))
type1$fwer_ci <- sprintf("%.3f (%.3f - %.3f)", type1$fwer,
    pmax(0, type1$fwer - 1.96 * type1$fwer_se),
    pmin(1, type1$fwer + 1.96 * type1$fwer_se))
knitr::kable(type1[, c("test", "network", "design", "replicates",
                       "fwer_ci", "no_region")],
    col.names = c("Test", "Network", "Design", "Replicates",
                  "FWER at 0.05 (95% CI)", "No region detected"),
    caption = "Family-wise error rate under the global null. The last column
    is the fraction of null replicates in which no region beyond the
    threshold existed, so the landscape test had nothing to reject.")

## ----null-ecdf, fig.height=4--------------------------------------------------
files <- list.files(val, pattern = "^null_pvalues_", full.names = TRUE)
nulls <- do.call(rbind, lapply(files, read.csv))
nulls <- nulls[is.finite(nulls$min_p), ]
ggplot(nulls, aes(x = min_p, colour = test)) +
    stat_ecdf(geom = "step", linewidth = 0.8) +
    geom_abline(slope = 1, intercept = 0, linetype = "dashed",
                colour = "grey40") +
    geom_vline(xintercept = 0.05, linetype = "dotted", colour = "grey40") +
    facet_wrap(~ network) +
    coord_cartesian(xlim = c(0, 1), ylim = c(0, 1)) +
    labs(x = "smallest p-value in a null replicate",
         y = "cumulative fraction of replicates", colour = NULL) +
    theme_minimal(base_size = 12) +
    theme(legend.position = "bottom")

## ----power--------------------------------------------------------------------
power <- read.csv(file.path(val, "power.csv"))
ggplot(power, aes(x = effect, y = power, colour = test, group = test)) +
    geom_line(linewidth = 0.8) +
    geom_point(size = 2) +
    geom_hline(yintercept = 0.05, linetype = "dotted", colour = "grey40") +
    facet_wrap(~ network) +
    scale_y_continuous(limits = c(0, 1)) +
    labs(x = "shift of the hub module (standard deviations)",
         y = "power at alpha = 0.05", colour = NULL) +
    theme_minimal(base_size = 12) +
    theme(legend.position = "bottom")

## ----layout-------------------------------------------------------------------
lay <- read.csv(file.path(val, "layout_sensitivity.csv"))
lay_summary <- read.csv(file.path(val, "layout_sensitivity_summary.csv"))
knitr::kable(lay_summary, digits = 3,
    col.names = c("Mean pairwise Jaccard of significant genes", "Layouts",
                  "TFCE smallest p (layout-free)"))
summary(lay[, c("min_p", "n_regions", "n_sig_genes", "module_recovered")])

## ----airway-load--------------------------------------------------------------
aw <- function(f) extdata("airway", f)
genes   <- read.delim(aw("airway_dex_genes.tsv"))
logcpm  <- as.matrix(read.delim(aw("airway_dex_logcpm.tsv"), row.names = 1))
samples <- read.delim(aw("airway_dex_samples.tsv"))
nodes   <- read.delim(aw("airway_string_nodes.tsv"))
edges   <- read.delim(aw("airway_string_edges.tsv"))
c(genes = nrow(genes), nodes = nrow(nodes), edges = nrow(edges))
table(samples$cell, samples$dex)

## ----airway-landscape, fig.height=6-------------------------------------------
set.seed(1)
land <- levi(
    expressionInput          = genes,
    networkCoordinatesInput  = nodes,
    networkInteractionsInput = edges,
    fileTypeInput            = "stg",
    geneSymbolInput          = "Symbol",
    readExpColumn            = readExpColumn("log2FoldChange-log2FoldChange"),
    signal_mode              = "logfc",
    logfc_k                  = 0.7,
    resolutionValueInput     = 30,
    smoothValueInput         = 40,
    n_perm                   = 199)
land$regions$summary[, c("Region", "Direction", "Cells", "Mass", "PSpatial",
                         "Significant")]
head(leviRegionGenes(land, top_n = 5), 10)

## ----airway-blocked-----------------------------------------------------------
groups <- samples$dex
tfce_blocked <- tryCatch(
    leviGraphTFCEInference(logcpm, groups, nodes, edges, fileTypeInput = "stg",
        test = "trt", control = "untrt", blocks = samples$cell),
    warning = function(w) conditionMessage(w))
tfce_blocked

## ----airway-unblocked---------------------------------------------------------
tfce <- leviGraphTFCEInference(logcpm, groups, nodes, edges,
    fileTypeInput = "stg", test = "trt", control = "untrt")
head(tfce$statistic[order(tfce$statistic$PGlobal), ], 8)

set.seed(1)
rep_land <- leviReplicateInference(logcpm, groups, test = "trt",
    control = "untrt", networkCoordinatesInput = nodes,
    networkInteractionsInput = edges, fileTypeInput = "stg",
    logfc_k = 0.7, resolutionValueInput = 30, smoothValueInput = 40)
rep_land$regions$summary[, c("Region", "Direction", "Mass", "PSpatial",
                             "Significant")]

## ----airway-rewiring----------------------------------------------------------
scores <- setNames(genes$stat, genes$Symbol)
set.seed(1)
rew <- leviGraphRewiringInference(scores, nodes, edges, fileTypeInput = "stg",
                                  threshold = 4, n_perm = 199)
# The STRING network of 78 responders is sparse (80 edges), so most
# clusters are single genes; show the connected ones.
clusters <- rew$regions$summary
clusters[clusters$Nodes >= 2,
         c("Region", "Direction", "Nodes", "Mass", "PeakGene", "PSpatial")]

## ----scripts, eval=FALSE------------------------------------------------------
# scripts <- system.file("scripts", package = "levi")
# # About an hour on four cores; writes inst/extdata/validation/*.csv
# system2("Rscript", c(file.path(scripts, "10-calibration.R"), "4", "300"))
# # A few minutes plus the STRING download; writes inst/extdata/airway/*.tsv
# system2("Rscript", file.path(scripts, "11-airway-preprocess.R"))

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