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

## ----installation, eval=FALSE-------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# BiocManager::install("IntrinsicScan")

## ----minimal-example----------------------------------------------------------
library(IntrinsicScan)

x <- rbind(
    feature_A = c(1, 1.1, 1.2, 1.3, 1.4, 5, 5.1, 5.2, 5.3, 5.4),
    feature_B = c(1, 2, 3, 2, 3, 2, 3, 4, 3, 4),
    feature_C = rep(7, 10)
)
colnames(x) <- paste0("S", seq_len(ncol(x)))
y <- rep(c("control", "case"), each = 5)

fit <- intrinsicScan(
    x,
    phenotype = y,
    positive = "case",
    verbose = FALSE
)

res <- intrinsicResults(fit)

res[, c(
    "feature_id",
    "valid",
    "reason",
    "threshold",
    "low_medoid",
    "high_medoid",
    "J",
    "positive_state",
    "exact_J_p",
    "exact_J_fdr",
    "relative_gap"
)]

## ----output-structure---------------------------------------------------------
SummarizedExperiment::assayNames(fit)

rd <- SummarizedExperiment::rowData(fit)
rd[, c(
    "scan_feature_id",
    "scan_threshold",
    "scan_J",
    "scan_exact_J_fdr",
    "scan_rank_J"
)]

## ----metadata-structure-------------------------------------------------------
names(S4Vectors::metadata(fit)$IntrinsicScan)

## ----biological-example, message=FALSE----------------------------------------
library(leukemiasEset)
data("leukemiasEset")

keep <- leukemiasEset$LeukemiaType %in% c("AML", "NoL")
aml_nol <- leukemiasEset[, keep]
aml_nol <- aml_nol[seq_len(min(1000L, nrow(aml_nol))), ]

table(as.character(aml_nol$LeukemiaType))

bio_fit <- intrinsicScan(
    aml_nol,
    phenotype = "LeukemiaType",
    positive = "AML",
    nResamples = 0L,
    keepStates = FALSE,
    verbose = FALSE
)

bio_res <- intrinsicResults(bio_fit)

utils::head(
    bio_res[, c(
        "feature_id",
        "threshold",
        "J",
        "positive_state",
        "exact_J_p",
        "exact_J_fdr",
        "relative_gap"
    )],
    10
)

## ----biological-example-plot--------------------------------------------------
top_id <- as.character(bio_res$feature_id[which(bio_res$valid)[1]])
top_threshold <- bio_res$threshold[match(top_id, bio_res$feature_id)]
top_values <- Biobase::exprs(aml_nol)[top_id, ]

group <- factor(
    as.character(aml_nol$LeukemiaType),
    levels = c("NoL", "AML")
)

boxplot(
    top_values ~ group,
    xlab = "Phenotype",
    ylab = "Expression",
    main = paste("Top IntrinsicScan feature:", top_id)
)
abline(h = top_threshold, lty = 2)

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