## ----opts, include = FALSE----------------------------------------------------
box::use(knitr[opts_chunk])
opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

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

## ----overview, echo = FALSE, fig.align = "center", out.width = "80%"----------
box::use(knitr[include_graphics])
include_graphics("Graphical_Abstract.png")

## ----fetch, message = FALSE---------------------------------------------------
# Cleanly add functions from various packages into the current namespace
box::use(
  TENxPBMCData[TENxPBMCData],
  SingleCellExperiment[counts, rowData],
  data.table[fread],
  tibble[deframe],
  stringr[str_split_1])

# Gene expression data (with gene symbols)
dat <- TENxPBMCData(dataset = "pbmc3k")
gex <- counts(dat) |> as("dgCMatrix")
dimnames(gex) <- list(rowData(dat)[["Symbol_TENx"]], dat[["Barcode"]])

# Signatures (consisting of gene symbols)
sig <- fread("signatures.tsv") |>
  deframe() |>
  lapply(\(x) stringr::str_split_1(x, ","))

## ----pre----------------------------------------------------------------------
box::use(
  scuttle[perCellQCMetrics, perCellQCFilters, normalizeCounts],
  Seqtometry[impute])

qc <- function(mat) {
  # Indices of mitochondrial genes
  mito_idx <- grep("^MT-", rownames(mat))

  # Calculate QC metrics
  res <- mat |>
    perCellQCMetrics(subsets = list("Mito" = mito_idx)) |>
    perCellQCFilters(sub.fields = "subsets_Mito_percent")

  # Remove low quality cells
  mat[, !res$discard]
}

imputed_gex <- gex |>
  qc() |>
  normalizeCounts() |>
  impute()

## ----score--------------------------------------------------------------------
box::use(
  future[plan],
  Seqtometry[score])

# May set up parallel backend if it is supported
# Will default to non-parallel (sequential) execution
plan("sequential")

# Allot 1 GiB of memory for global variables
options(future.globals.maxSize = 1024 ^ 3)

# Scoring with lymphocyte signatures
scores <- sig[c("T cell", "B cell")] |>
  score(imputed_gex, signatures = _)

## ----plot, fig.align = "center", out.width = "50%"----------------------------
box::use(ggplot2[...])

# A flow cytometry style theme for plots
theme_custom <- theme_bw() + theme(
  aspect.ratio = 1,
  panel.grid = element_blank(),
  panel.background = element_blank(),
  panel.border = element_rect(color = "black", fill = NA, linewidth = 2),
  axis.title = element_text(size = 32, face = "italic"),
  axis.text = element_text(size = 16))

# Utility function for generating flow cytometry style plots
biaxial_plot <- function(d, x, y) {
  lim <- c(-0.15, 1.15)
  brk <- seq(0, 1, 0.25)
  lbl <- c("0", "0.25", "0.5", "0.75", "1")
  ggplot(d, aes(x = .data[[x]], y = .data[[y]])) +
    geom_density_2d(color = "black", h = c(0.2, 0.2),
      contour_var = "ndensity", binwidth = 0.125) +
    scale_x_continuous(limits = lim, breaks = brk, labels = lbl) +
    scale_y_continuous(limits = lim, breaks = brk, labels = lbl) +
    theme_custom
}

# Utility function to draw a rectangular region
rect_region <- function(lims) {
  do.call(annotate, c(lims, geom = "rect",
    color = "blue", fill = NA, linewidth = 2))
}

# Limits defining a rectangular region for double negative cells
dn_lims <- list(xmin = -0.075, xmax = 0.5, ymin = 0, ymax = 0.5)

# Draw the biaxial plot with the rectangular region
biaxial_plot(scores, "T cell", "B cell") + rect_region(dn_lims)

## ----gate, fig.align = "center", out.width = "50%"----------------------------
# Subset the imputed data, keeping only the cells in the blue region
gated_imputed_gex <- with(scores,
  imputed_gex[, `T cell` < dn_lims$xmax & `B cell` < dn_lims$ymax])

# Score the subsetted data with the remaining signatures
gated_scores <- sig[c("NK cell", "Monocyte")] |>
  Seqtometry::score(gated_imputed_gex, signatures = _)

# Plot the new scores
biaxial_plot(gated_scores, "Monocyte", "NK cell")

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