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

## ----install, eval = FALSE----------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE)) {
#     install.packages("BiocManager")
# }
# BiocManager::install("thyroidBRS")

## ----setup--------------------------------------------------------------------
library(thyroidBRS)

## ----simulate-data------------------------------------------------------------
set.seed(1)
genes <- na.omit(brs_genes$current_symbol)
samples <- paste0("sample_", seq_len(60))

# Two groups with a real difference across the signature, so the example
# shows the classifier doing something rather than fitting noise. The 24
# samples without a label are a mix of both kinds, as a real cohort is.
group <- rep(c("BRAF_V600E", "RAS", "unknown"), c(24, 12, 24))
up_in_ras <- brs_genes$up_in[!is.na(brs_genes$current_symbol)] == "RAS"
shift <- ifelse(up_in_ras, 1, -1)
side <- ifelse(group == "BRAF_V600E", -0.5, 0.5)
side[group == "unknown"] <- rep(c(-0.5, 0.5), length.out = 24)

expr <- matrix(rnorm(length(genes) * length(samples)),
               nrow = length(genes), dimnames = list(genes, samples))
expr <- expr + outer(shift, side)

# Mutation labels for the 36 samples that have exome data...
labels <- setNames(group[seq_len(36)], samples[seq_len(36)])
# ...and 24 with expression only.
unlabeled <- setdiff(samples, names(labels))

## ----fit----------------------------------------------------------------------
fit <- brs_fit(expr, labels)
fit

## ----classify-new-------------------------------------------------------------
preds_new <- predict(fit, expr[, unlabeled])
head(preds_new)
table(preds_new$brs_class)

## ----validate-----------------------------------------------------------------
preds_all <- predict(fit, expr)
validate_brs(preds_all, labels, fit = fit)

## ----crossval-----------------------------------------------------------------
set.seed(2)
n_folds <- 4
fold <- sample(rep(seq_len(n_folds), length.out = length(labels)))

# Fit on all folds but one, score the held-out fold, and name each call by
# its sample.
cv_fold <- function(k) {
    f <- brs_fit(expr, labels[fold != k])
    held <- names(labels)[fold == k]
    setNames(predict(f, expr[, held, drop = FALSE])$brs_class, held)
}
cv <- unlist(lapply(seq_len(n_folds), cv_fold))[names(labels)]

truth <- ifelse(labels == "BRAF_V600E", "Braf-like", "Ras-like")
mean(cv == truth)

## ----one-shot-----------------------------------------------------------------
all_preds <- brs_score(expr, labels)
head(all_preds)

## ----containers, message = FALSE----------------------------------------------
library(SummarizedExperiment)

# colData() has to carry driver-mutation status. Whatever you name the
# column, it must not be the published BRAF-like/RAS-like class: that is
# the output of this same classifier, and fitting on it is circular.
# Samples without exome data get NA and are simply not used for the fit.
se <- SummarizedExperiment(
    assays = list(logexpr = expr),
    colData = DataFrame(driver = replace(group, group == "unknown", NA),
                        row.names = samples)
)

fit_se <- brs_fit(se, "driver", assay = "logexpr")
identical(fit_se, fit)
head(predict(fit_se, se, assay = "logexpr"))

## ----tcga, eval = FALSE-------------------------------------------------------
# library(TCGAbiolinks)
# 
# # 1. A log2(TPM + 1) matrix of primary tumors, gene symbols x patient barcode.
# q <- GDCquery(
#     project = "TCGA-THCA",
#     data.category = "Transcriptome Profiling",
#     data.type = "Gene Expression Quantification",
#     workflow.type = "STAR - Counts"
# )
# GDCdownload(q)
# se <- GDCprepare(q)
# # keep assay "tpm_unstranded", primary tumors (-01), collapse duplicated
# # gene symbols, then log2(x + 1)  ->  expr
# 
# # 2. Reference labels: driver mutation, NOT the published class.
# ref <- read.csv(system.file("extdata", "thca_reference.csv",
#                             package = "thyroidBRS"))
# mut <- setNames(ref$driver_group, ref$patient)
# mut <- mut[mut %in% c("BRAF_V600E", "RAS")]
# 
# # 3. Fit, then score everyone.
# fit <- brs_fit(expr, mut)
# preds <- predict(fit, expr)
# 
# # 4. Validate against something that did not enter the fit.
# published <- setNames(ref$published_class, ref$patient)
# validate_brs(preds, published[!is.na(published)], fit = fit)

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