## ----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)

spe <- readRDS(system.file("extdata", "mel_cosmx_subset.rds", package = "PACE"))
spe
table(spe$cellType)
table(unique(data.frame(image = spe$image, arm = spe$Responder))$arm)

## ----fit----------------------------------------------------------------------
fit <- paceFit(spe,
               celltype_col     = "cellType",
               condition_col    = "Responder",
               image_col        = "image",
               kernel_per_image = TRUE,
               image_re         = "intercept",
               dispersion       = "nb1",
               verbose          = FALSE)
fit

## ----respterm-----------------------------------------------------------------
paceParams(fit)$resp_term

## ----decomp-------------------------------------------------------------------
varianceDecomposition(fit) |>
  group_by(focal) |>
  summarise(`Spatial %`           = round(median(`Spatial %`), 4),
            `Responder spatial %` = round(median(`Responder spatial %`), 4),
            `Spillover %`         = round(median(`Spillover %`), 2),
            .groups = "drop")

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

## ----pairs-responder, fig.width = 6, fig.height = 4.6-------------------------
pairVariance(fit, block = "responder") |>
  arrange(desc(val)) |>
  head(5)

plotPairHeatmap(fit, block = "responder")

## ----pairs-spatial, fig.width = 6, fig.height = 4.6---------------------------
plotPairHeatmap(fit, block = "spatial")

## ----calls--------------------------------------------------------------------
ns <- neighbourSlopes(fit)
resp <- subset(ns, grepl("^ResponderPD:", term) & lfsr < 0.05)
nrow(resp)

resp[order(resp$lfsr), c("gene", "focal", "neighbour", "estimate_shrunk", "lfsr")] |>
  head(8)

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

## ----drivers-table------------------------------------------------------------
topDrivers(fit)[["Tumour_T_Cell"]]$scores |>
  select(rank, gene, driver_score, b_clean, spec, focal_mean, lfsr) |>
  head(8)

## ----spp1---------------------------------------------------------------------
subset(ns, gene == "SPP1" & focal == "Macrophage" & neighbour == "Tumour")

## ----proximity, fig.width = 7, fig.height = 4---------------------------------
plotProximity(fit, spe, "SPP1", "Macrophage", "Tumour",
              breaks = c(0, 10, 20, 30, 40, Inf), condition = TRUE)

## ----curve, fig.width = 7, fig.height = 5-------------------------------------
plotResponseCurve(fit, spe, "SPP1", "Macrophage", "Tumour")

## ----map, fig.width = 9, fig.height = 8---------------------------------------
plotResponseMap(fit, spe, "SPP1", "Macrophage", "Tumour",
                images = c(nonPD = "32157_18", PD = "32156_17"))

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