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

## ----install, eval = FALSE----------------------------------------------------
# if (!require("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# BiocManager::install("rnaSentry")

## ----data-intake--------------------------------------------------------------
library(rnaSentry)
library(SummarizedExperiment)

cnt <- matrix(rpois(600, lambda = 200), nrow = 60, ncol = 10,
              dimnames = list(paste0("gene", 1:60), paste0("S", 1:10)))
cnt <- cnt * 1.0
cnt[1, ] <- 0L                 # gene1 has no counts anywhere (empty row)
cnt[2, 1] <- cnt[2, 1] + 0.5   # a non-integer value
colnames(cnt)[10] <- "S1"      # a duplicate sample identifier
cd <- S4Vectors::DataFrame(condition = rep(c("tumor", "normal"), each = 5),
                           row.names = paste0("S", 1:10))

se_in <- load_counts(cnt, cd)
S4Vectors::metadata(se_in)$flags

## ----qc-explore---------------------------------------------------------------
set.seed(1)
counts <- matrix(rpois(600, lambda = 200), nrow = 60, ncol = 10,
                  dimnames = list(paste0("gene", 1:60), paste0("S", 1:10)))
counts[, 1] <- counts[, 1] * 50L
coldata <- S4Vectors::DataFrame(condition = rep(c("tumor", "normal"), each = 5),
                                row.names = paste0("S", 1:10))
se <- SummarizedExperiment(assays = list(counts = counts), colData = coldata)

qc <- qc_explore(se)
qc

## ----sex-check----------------------------------------------------------------
genes <- c("XIST", "RPS4Y1", "DDX3Y", "KDM5D")
expr <- matrix(c(950, 5, 5, 5,   # sample looks female
                  5, 400, 380, 410),  # sample looks male
                nrow = 4, dimnames = list(genes, c("S1", "S2")))
coldata2 <- S4Vectors::DataFrame(sex = c("F", "F"), row.names = c("S1", "S2")) # S2 mislabeled
se2 <- SummarizedExperiment(assays = list(counts = expr), colData = coldata2)

sex_check(se2)

## ----run-pipeline-------------------------------------------------------------
set.seed(1)
n_genes <- 120
n_samples <- 120
counts <- matrix(stats::rpois(n_genes * n_samples, lambda = 500),
                 nrow = n_genes, ncol = n_samples,
                 dimnames = list(paste0("gene", 1:n_genes),
                                 paste0("S", 1:n_samples)))
sig_expr <- colMeans(counts[1:5, , drop = FALSE])
risk <- scale(sig_expr)[, 1] * 0.4
event_time <- stats::rexp(n_samples, rate = 0.03 * exp(0.8 * risk))
censor_time <- stats::rexp(n_samples, rate = 0.02)
coldata <- S4Vectors::DataFrame(time = pmin(event_time, censor_time),
                                event = as.integer(event_time < censor_time),
                                condition = factor(rep(c("A", "B"), length.out = n_samples)),
                                batch = factor(rep(c("B1", "B2", "B3"), length.out = n_samples)),
                                row.names = colnames(counts))
se_all <- SummarizedExperiment(assays = list(counts = counts), colData = coldata)

discovery <- se_all[, 1:72]
external <- se_all[, 73:120]

run <- run_rnaSentry(discovery, time_col = "time", event_col = "event",
                     top_n = 5, repeats = 2, folds = 3, seed = 7,
                     design_vars = c("condition", "batch"),
                     design_terms = "condition",
                     render_report = FALSE)
run

## ----unlock-signature---------------------------------------------------------
# Release the lock so the pipeline can be run again in this session.
run$stages$build_signature <-
  lock_signature(run$stages$build_signature, lock = FALSE)

## ----stage-flags--------------------------------------------------------------
run$stages$build_signature$flags

## ----ext-validate-------------------------------------------------------------
val <- validate_external(run$stages$build_signature, external,
                         cutpoint = run$stages$km_curve$cutpoint)
val

## ----intro-ext-km, fig.width = 7, fig.height = 5, fig.align = "center", fig.cap = "External-cohort Kaplan-Meier curves for the low- and high-risk groups, scored against the locked discovery cutpoint without re-estimation. The log-rank p-value tests the split on held-out data."----
plot(val)

## ----wrong-scale, warning = TRUE----------------------------------------------
prelogged <- log2(as.matrix(assay(se_all, "counts")) + 1)
se_wrong <- se_all
assay(se_wrong, "counts") <- prelogged
build_signature(se_wrong, "time", "event", top_n = 5, repeats = 1,
                folds = 2, seed = 7, min_events_per_parameter = 1)

## ----underpowered, warning = TRUE---------------------------------------------
small <- se_all[, 1:40]
build_signature(small, "time", "event", top_n = 10, repeats = 1,
                folds = 2, seed = 7)

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