---
title: "PACE package manual"
author: "Elijah Willie"
date: "`r Sys.Date()`"
output:
    BiocStyle::html_document:
        toc: true
        toc_float: true
vignette: >
  %\VignetteIndexEntry{PACE package manual}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

```{css scrollable-output, echo = FALSE}
/* Long console output (the shrinkage progress log, sessionInfo) scrolls in
   place rather than pushing the prose off the page. Short blocks are
   unaffected: max-height only engages once the block is taller. */
pre {
  max-height: 340px;
  overflow-y: auto;
}
```

# Installation

```r
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("PACE")
```

# Introduction

`r Biocpkg("PACE")` (Proximity-Associated Changes in Expression) is a framework
for imaging-based spatial transcriptomics (such as Xenium and CosMx) that
quantifies how a cell's gene expression changes with proximity to *specific*
neighbouring cell types. Rather than asking only which genes are spatially
variable, PACE asks a directed question: when a macrophage sits next to tumour
cells, which of its genes go up or down, and by how much?

PACE fits a hierarchical negative binomial mixed model with partial pooling
across cell types, so that each focal-neighbour relationship borrows strength
from the others and noisy pairs are regularised towards the consensus. A per-cell
contamination term absorbs the short-range ambient signal that leaks between
adjacent cells through segmentation error, so the estimated neighbour effect
reflects biology rather than transcript misassignment.

# Overview of the PACE workflow

A PACE analysis has three stages, all driven by the single call `paceFit()`:

1. **Neighbourhood construction.** For every cell, a Gaussian biological kernel
   summarises the abundance of each neighbouring cell type, and a short-range
   exponential technical kernel captures the ambient contamination field.
2. **Model fitting.** A streaming penalised quasi-likelihood fit estimates, for
   every gene and every focal cell type, a partially pooled slope on each
   neighbour type's abundance, jointly with a per-cell contamination loading and
   gene-wise overdispersions. The neighbour slopes are then stabilised by
   multivariate adaptive shrinkage.
3. **Interpretation.** The fit is read out three ways: shrunken neighbour slopes
   with a local false sign rate (`neighbourSlopes()`), a per-gene variance
   decomposition (`varianceDecomposition()`), and per-pair driver-score tables
   ranking the genes that mediate each relationship (`topDrivers()`).

# The breast cancer dataset

The package ships a small worked example: a spatial crop of a human breast cancer
Xenium section (10x Genomics; nucleus segmentation, cell types annotated with
scClassify), centred on a tumour-macrophage interface. It contains `r "7,898"`
cells across eight cell types (B cell, dendritic cell, endothelial, macrophage,
myoepithelial, stromal, T cell, and tumour) and a 278-gene panel, packaged as a
`r Biocpkg("SpatialExperiment")`. This region reproduces the macrophage
reprogramming that PACE identifies on the full section.

```{r 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)
```

# Fitting the model

`paceFit()` reads the counts, spatial coordinates, and cell-type labels from the
`SpatialExperiment` and runs the whole pipeline. On this crop it takes about half
a minute.

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

## The pipeline step by step

`paceFit()` is a convenience wrapper. The same result is produced by calling each
step explicitly, which lets you inspect the fitted model and re-run the
downstream layers (shrinkage, decomposition, drivers) without refitting:

```{r 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))
```

The fit-construction primitives expose the neighbourhood and contamination field
directly, for inspection independently of a fit:

```{r 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
```

# Variance decomposition

`plotDecomposition()` shows, for each focal cell type, the share of expression
variance attributable to cell-type identity, spatial cell state (proximity
effects), and contamination (spillover), pooled over genes. The left panel is the
full stacked bar; the right zooms in on the two small spatial components.

```{r decomp, fig.width = 11, fig.height = 4.5}
plotDecomposition(fit)
```

Beyond the dominant cell-type identity component, B cells, stromal cells and
macrophages carry the largest spatial cell-state contributions, while dendritic
and myoepithelial cells carry none on this crop. The underlying per-gene table
is available with `varianceDecomposition(fit)`.

# Per-cell contamination

The spillover component above is a variance share per gene. `cellContamination()`
gives the per-cell view of the same term: the loading `rho_i` and, more usefully,
`contamFraction`, the share of a cell's expected counts attributed to its local
ambient field. A high fraction marks a cell whose profile is largely explained by
its neighbours rather than by its own type, which is what a segmentation or
transcript-assignment error looks like.

```{r contamination}
cc <- cellContamination(fit)
head(cc, 3)
round(tapply(cc$contamFraction, cc$celltype, median), 3)
```

Read the fraction rather than the loading. The ambient field is cross-cell-type,
so a cell with no differently-typed neighbour inside the technical kernel has no
ambient signal at all; its `rho_i` is unidentified and shrinks to the empirical
Bayes prior mean, giving every such cell the same apparently middling loading.
On this crop that is most of the tumour interior:

```{r 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)
```

`contamFraction` reports those cells as zero, which is correct; `rho` alone does
not.

# Pairwise spatial interactions

`plotPairHeatmap()` gives the percentage of each focal type's total variance
contributed by spatial interaction with each neighbour, computed as the focal's
spatial share split across neighbours by a normalised Pratt attribution.
**Tumour as a neighbour** acts the most widely across the microenvironment,
reaching four focal types where no other neighbour reaches more than two. The
single strongest pair is B cells next to T cells.

```{r pairs, fig.width = 6.5, fig.height = 5}
plotPairHeatmap(fit)
```

`pairVariance()` returns those same numbers as a long table. The heatmap is
drawn from it, so the two cannot disagree, and the table is the easier thing to
rank or pass on:

```{r pairvariance}
pairVariance(fit) |>
  arrange(desc(val)) |>
  head(5)
```

For a cohort with a condition, `block = "responder"` gives the same attribution
for the responder-by-proximity block rather than the baseline spatial one.

# Macrophage–tumour drivers

`plotDrivers()` ranks the genes mediating a relationship by a driver score
(combining the shrunken slope, its cell-type specificity, and its expression
level), alongside their per-gene single-frame decomposition. For the
macrophage-tumour pair, the top drivers are the tissue-resident marker **MRC1**
(CD206), reduced near tumour, and the lipid-associated marker **APOC1**, elevated
near tumour.

```{r drivers, fig.width = 10, fig.height = 4.5}
plotDrivers(fit, "Macrophage", "Tumour")
```

The shrunken slopes confirm the opposing directions:

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

# Visualising the proximity effect

`plotProximity()` bins each macrophage by its number of tumour neighbours
(within 30 um) and shows each gene's raw counts per bin. Because these genes are
zero-inflated the box median sits at zero in most bins, so the per-bin mean is
overlaid as a point to make the trend explicit. Macrophage MRC1 falls and APOC1
rises with tumour proximity.

```{r gradient, fig.width = 11, fig.height = 4.6}
plotProximity(fit, spe, c("MRC1", "APOC1"), "Macrophage", "Tumour")
```

# Biological interpretation

On this section PACE recovers a coordinated shift in macrophage state at the
tumour interface: the tissue-resident marker MRC1 (CD206) is reduced and the
lipid-associated, immunosuppressive marker APOC1 is elevated in macrophages near
tumour cells. Because PACE separates this proximity effect from the transcript
contamination that would otherwise inflate tumour-marker signal in the
macrophages, the drivers it promotes are genuine macrophage genes rather than
tumour genes bleeding across cell boundaries.

# Keeping results with the data

A `PACEFit` is its own object, because a fit on a real cohort is large and
usually outlives the session that made it. When you would rather the results
travelled with the data, `stashPACE()` copies the readouts into
`metadata(spe)`, leaving the solver state behind:

```{r stash}
spe <- stashPACE(spe, fit)
names(metadata(spe)$PACE)
```

The three reporting layers are there, along with enough provenance to say which
fit produced them. The solver state is not, which is why the stash is a
fraction of the size of the fit:

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

# Saving a fit

A fit holds the fitted means and the contamination log-offset as full
`cells x genes` matrices. They make the object large: even this small crop is

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

and a whole section runs to several hundred megabytes. They do not have to be
kept. `paceDecompose()` and `paceDrivers()` read only per-cell-type summaries of
the fitted means, which the fit also stores, so a fit saved without the
matrices can still be re-decomposed (from the counts in `spe`) and re-scored. (For an older fit saved
without those summaries, pass `spe` to both, and they are rebuilt exactly from
the fit and the counts.)

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

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

Every readout above still works on the stripped fit, and re-decomposing it
reproduces the table it already carries:

```{r 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])))
```

# Session info

```{r sessioninfo}
sessionInfo()
```