## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE, comment = "#>", message = FALSE, warning = FALSE,
  fig.width = 6, fig.height = 4, dpi = 90)

## ----load, message = FALSE----------------------------------------------------
library(PACE)
library(SpatialExperiment)
library(ggplot2)
library(dplyr)
library(tidyr)

spe <- readRDS(system.file("extdata", "bc_xenium_subset.rds", package = "PACE"))
spe
table(spe$cellType)

## ----fit----------------------------------------------------------------------
fit <- paceFit(spe,
               celltype_col  = "cellType",
               contamination = "percell_hc",   # per-cell contamination correction
               dispersion    = "nb1",
               verbose       = FALSE)
fit

## ----staged-------------------------------------------------------------------
staged <- paceModel(spe, celltype_col = "cellType", verbose = FALSE) |>
          paceShrink()            |>   # shrink the neighbour slopes
          paceDecompose(spe)      |>   # variance decomposition
          paceDrivers(spe)             # per-pair driver scores

# the same object, by either route
all.equal(varianceDecomposition(staged), varianceDecomposition(fit))

## ----primitives---------------------------------------------------------------
kern <- buildNeighbourhood(spe, "cellType")   # Gaussian K_bio + technical K_tech
W    <- ambientField(spe, "cellType")         # sparse E^tech contamination field

vapply(kern, function(x) paste(dim(x), collapse = " x "), character(1))
dim(W)
head(anchorGenes(fit, spe), 3)                # per-cell-type contamination anchors

## ----decomp, fig.width = 11, fig.height = 4.5---------------------------------
plotDecomposition(fit)

## ----contamination------------------------------------------------------------
cc <- cellContamination(fit)
head(cc, 3)
round(tapply(cc$contamFraction, cc$celltype, median), 3)

## ----contamination-zero-------------------------------------------------------
zero <- cc$contamFraction == 0
c(cells = sum(zero), distinct_rho = length(unique(round(cc$rho[zero], 12))))
round(tapply(zero, cc$celltype, mean), 2)

## ----pairs, fig.width = 6.5, fig.height = 5-----------------------------------
plotPairHeatmap(fit)

## ----pairvariance-------------------------------------------------------------
pairVariance(fit) |>
  arrange(desc(val)) |>
  head(5)

## ----drivers, fig.width = 10, fig.height = 4.5--------------------------------
plotDrivers(fit, "Macrophage", "Tumour")

## ----flagships----------------------------------------------------------------
neighbourSlopes(fit) |>
  filter(focal == "Macrophage", neighbour == "Tumour",
         gene %in% c("MRC1", "APOC1")) |>
  select(gene, estimate_shrunk, lfsr)

## ----gradient, fig.width = 11, fig.height = 4.6-------------------------------
plotProximity(fit, spe, c("MRC1", "APOC1"), "Macrophage", "Tumour")

## ----stash--------------------------------------------------------------------
spe <- stashPACE(spe, fit)
names(metadata(spe)$PACE)

## ----stash-size---------------------------------------------------------------
c(fit = object.size(fit), stashed = object.size(metadata(spe)$PACE)) |>
  format(units = "MB")

## ----fit-size-----------------------------------------------------------------
round(as.numeric(object.size(fit)) / 1e6, 1)   # MB

## ----fit-strip----------------------------------------------------------------
fit_small <- compactFit(fit)
round(as.numeric(object.size(fit_small)) / 1e6, 1)   # MB

## ----fit-save-----------------------------------------------------------------
path <- file.path(tempdir(), "fit.rds")
saveRDS(fit_small, path, compress = "xz")
round(file.size(path) / 1e6, 1)                      # MB on disk

## ----fit-strip-check----------------------------------------------------------
redone <- paceDecompose(fit_small, spe)
stored <- varianceDecomposition(fit)
num <- vapply(stored, is.numeric, logical(1))
max(abs(as.matrix(varianceDecomposition(redone)[num]) - as.matrix(stored[num])))

## ----sessioninfo--------------------------------------------------------------
sessionInfo()