PACE package manual
2026-10-15
1 Installation
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("PACE")2 Introduction
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.
3 Overview of the PACE workflow
A PACE analysis has three stages, all driven by the single call paceFit():
- 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.
- 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.
- 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()).
4 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 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 SpatialExperiment. This region reproduces the macrophage reprogramming that PACE identifies on the full section.
library(PACE)
library(SpatialExperiment)
library(ggplot2)
library(dplyr)
library(tidyr)
spe <- readRDS(system.file("extdata", "bc_xenium_subset.rds", package = "PACE"))
spe
#> class: SpatialExperiment
#> dim: 278 7898
#> metadata(0):
#> assays(1): counts
#> rownames(278): SEC11C DAPK3 ... NOSTRIN CD1C
#> rowData names(0):
#> colnames(7898): 442 444 ... 99842 99848
#> colData names(2): cellType sample_id
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):
#> spatialCoords names(2) : x y
#> imgData names(0):
table(spe$cellType)
#>
#> B_Cell Dendritic_Cell Endothelial Macrophage Myoepithelial
#> 330 32 499 591 47
#> Stromal T_Cell Tumour
#> 681 1076 46425 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.
fit <- paceFit(spe,
celltype_col = "cellType",
contamination = "percell_hc", # per-cell contamination correction
dispersion = "nb1",
verbose = FALSE)
fit
#> class: PACEFit
#> cell types (8): B_Cell, Dendritic_Cell, Endothelial, Macrophage, Myoepithelial, Stromal, T_Cell, Tumour
#> kernels: h_bio = 30 um, h_tech = 5 um | contamination: percell_hc; dispersion: nb1
#> pipeline: model -> shrink -> decompose -> drivers
#> neighbour slopes: 17792 rows (61 at lfsr < 0.05, 0.3 expected to have the wrong sign)5.1 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:
staged <- paceModel(spe, celltype_col = "cellType", verbose = FALSE) |>
paceShrink() |> # shrink the neighbour slopes
paceDecompose(spe) |> # variance decomposition
paceDrivers(spe) # per-pair driver scores
#> - Computing 278 x 313 likelihood matrix.
#> - Likelihood calculations took 0.19 seconds.
#> - Fitting model with 313 mixture components.
#> - Model fitting took 0.35 seconds.
#> - Computing posterior matrices.
#> - Computation allocated took 0.00 seconds.
#> - Computing 278 x 92 likelihood matrix.
#> - Likelihood calculations took 0.05 seconds.
#> - Fitting model with 92 mixture components.
#> - Model fitting took 0.07 seconds.
#> - Computing posterior matrices.
#> - Computation allocated took 0.00 seconds.
#> - Computing 278 x 404 likelihood matrix.
#> - Likelihood calculations took 0.38 seconds.
#> - Fitting model with 404 mixture components.
#> - Model fitting took 0.58 seconds.
#> - Computing posterior matrices.
#> - Computation allocated took 0.00 seconds.
#> - Computing 278 x 417 likelihood matrix.
#> - Likelihood calculations took 0.25 seconds.
#> - Fitting model with 417 mixture components.
#> - Model fitting took 1.02 seconds.
#> - Computing posterior matrices.
#> - Computation allocated took 0.01 seconds.
#> - Computing 278 x 92 likelihood matrix.
#> - Likelihood calculations took 0.01 seconds.
#> - Fitting model with 92 mixture components.
#> - Model fitting took 0.15 seconds.
#> - Computing posterior matrices.
#> - Computation allocated took 0.00 seconds.
#> - Computing 278 x 391 likelihood matrix.
#> - Likelihood calculations took 0.47 seconds.
#> - Fitting model with 391 mixture components.
#> - Model fitting took 0.46 seconds.
#> - Computing posterior matrices.
#> - Computation allocated took 0.00 seconds.
#> - Computing 278 x 430 likelihood matrix.
#> - Likelihood calculations took 0.27 seconds.
#> - Fitting model with 430 mixture components.
#> - Model fitting took 0.93 seconds.
#> - Computing posterior matrices.
#> - Computation allocated took 0.01 seconds.
#> - Computing 278 x 628 likelihood matrix.
#> - Likelihood calculations took 0.41 seconds.
#> - Fitting model with 628 mixture components.
#> - Model fitting took 1.26 seconds.
#> - Computing posterior matrices.
#> - Computation allocated took 0.01 seconds.
#> - Computing 278 x 404 likelihood matrix.
#> - Likelihood calculations took 0.25 seconds.
#> - Fitting model with 404 mixture components.
#> - Model fitting took 2.48 seconds.
#> - Computing posterior matrices.
#> - Computation allocated took 0.02 seconds.
#> - Computing 278 x 590 likelihood matrix.
#> - Likelihood calculations took 0.39 seconds.
#> - Fitting model with 590 mixture components.
#> - Model fitting took 3.04 seconds.
#> - Computing posterior matrices.
#> - Computation allocated took 0.01 seconds.
# the same object, by either route
all.equal(varianceDecomposition(staged), varianceDecomposition(fit))
#> [1] TRUEThe fit-construction primitives expose the neighbourhood and contamination field directly, for inspection independently of a fit:
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))
#> K_bio K_tech
#> "7898 x 8" "7898 x 8"
dim(W)
#> NULL
head(anchorGenes(fit, spe), 3) # per-cell-type contamination anchors
#> cellType n_anchors
#> 1 B_Cell 129
#> 2 Dendritic_Cell 141
#> 3 Endothelial 129
#> top_anchors
#> 1 ERBB2, KRT7, POSTN, ANKRD30A, FOXA1, EPCAM, SCD, GATA3, CCND1, MLPH
#> 2 ERBB2, KRT7, ANKRD30A, FOXA1, EPCAM, SCD, GATA3, CCND1, MLPH, FASN
#> 3 ERBB2, KRT7, ANKRD30A, FOXA1, EPCAM, SCD, GATA3, MLPH, FASN, TACSTD26 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.
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).
7 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.
cc <- cellContamination(fit)
head(cc, 3)
#> cell celltype rho contamFraction
#> 1 442 Tumour 0.11504406 0.000000000
#> 2 444 B_Cell 0.00282706 0.001557513
#> 3 446 Macrophage 0.58779262 0.123800291
round(tapply(cc$contamFraction, cc$celltype, median), 3)
#> B_Cell Dendritic_Cell Endothelial Macrophage Myoepithelial
#> 0.091 0.083 0.018 0.073 0.038
#> Stromal T_Cell Tumour
#> 0.038 0.072 0.000Read 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:
zero <- cc$contamFraction == 0
c(cells = sum(zero), distinct_rho = length(unique(round(cc$rho[zero], 12))))
#> cells distinct_rho
#> 4589 1
round(tapply(zero, cc$celltype, mean), 2)
#> B_Cell Dendritic_Cell Endothelial Macrophage Myoepithelial
#> 0.10 0.16 0.26 0.19 0.23
#> Stromal T_Cell Tumour
#> 0.30 0.18 0.84contamFraction reports those cells as zero, which is correct; rho alone does
not.
8 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.
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:
pairVariance(fit) |>
arrange(desc(val)) |>
head(5)
#> # A tibble: 5 × 3
#> focal neighbour val
#> <chr> <chr> <dbl>
#> 1 B_Cell T_Cell 2.73
#> 2 Stromal Tumour 1.59
#> 3 Macrophage Tumour 0.756
#> 4 Endothelial Tumour 0.372
#> 5 Tumour T_Cell 0.342For a cohort with a condition, block = "responder" gives the same attribution
for the responder-by-proximity block rather than the baseline spatial one.
9 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.
plotDrivers(fit, "Macrophage", "Tumour")
The shrunken slopes confirm the opposing directions:
neighbourSlopes(fit) |>
filter(focal == "Macrophage", neighbour == "Tumour",
gene %in% c("MRC1", "APOC1")) |>
select(gene, estimate_shrunk, lfsr)
#> gene estimate_shrunk lfsr
#> 1 APOC1 0.1301689 1.077054e-17
#> 2 MRC1 -0.1494760 1.554312e-1510 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.
plotProximity(fit, spe, c("MRC1", "APOC1"), "Macrophage", "Tumour")
11 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.
12 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:
spe <- stashPACE(spe, fit)
names(metadata(spe)$PACE)
#> [1] "neighbourSlopes" "varianceDecomposition" "topDrivers"
#> [4] "cellTypes" "params" "packageVersion"
#> [7] "stashedAt"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:
c(fit = object.size(fit), stashed = object.size(metadata(spe)$PACE)) |>
format(units = "MB")
#> fit stashed
#> "7162712" "2739200"13 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
round(as.numeric(object.size(fit)) / 1e6, 1) # MB
#> [1] 7.2and 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.)
fit_small <- compactFit(fit)
round(as.numeric(object.size(fit_small)) / 1e6, 1) # MB
#> [1] 7.2path <- file.path(tempdir(), "fit.rds")
saveRDS(fit_small, path, compress = "xz")
round(file.size(path) / 1e6, 1) # MB on disk
#> [1] 1Every readout above still works on the stripped fit, and re-decomposing it reproduces the table it already carries:
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])))
#> [1] 014 Session info
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.5 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] tidyr_1.3.2 dplyr_1.2.1
#> [3] ggplot2_4.0.3 SpatialExperiment_1.23.0
#> [5] SingleCellExperiment_1.35.2 SummarizedExperiment_1.43.0
#> [7] Biobase_2.73.2 GenomicRanges_1.65.4
#> [9] Seqinfo_1.3.2 IRanges_2.47.5
#> [11] S4Vectors_0.51.10 BiocGenerics_0.59.12
#> [13] generics_0.1.4 MatrixGenerics_1.25.0
#> [15] matrixStats_1.5.0 PACE_0.99.20
#> [17] BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] gtable_0.3.6 rjson_0.2.23 xfun_0.61
#> [4] bslib_0.12.0 lattice_0.23-1 vctrs_0.7.3
#> [7] tools_4.6.1 tibble_3.3.1 pkgconfig_2.0.3
#> [10] Matrix_1.7-6 SQUAREM_2026.1 RColorBrewer_1.1-3
#> [13] S7_0.2.2 assertthat_0.2.1 lifecycle_1.0.5
#> [16] truncnorm_1.0-9 stringr_1.6.0 compiler_4.6.1
#> [19] farver_2.1.2 tinytex_0.61 htmltools_0.5.9
#> [22] sass_0.4.10 yaml_2.3.12 pillar_1.11.1
#> [25] jquerylib_0.1.4 rmeta_3.0 cachem_1.1.0
#> [28] DelayedArray_0.39.8 dbscan_1.2.7 magick_2.9.1
#> [31] abind_1.4-8 tidyselect_1.2.1 digest_0.6.39
#> [34] stringi_1.8.9 mvtnorm_1.4-2 purrr_1.2.2
#> [37] bookdown_0.48 ashr_2.2-63 labeling_0.4.3
#> [40] fastmap_1.2.0 grid_4.6.1 cli_3.6.6
#> [43] invgamma_1.2 SparseArray_1.13.4 magrittr_2.0.5
#> [46] patchwork_1.3.2 S4Arrays_1.13.2 utf8_1.2.6
#> [49] dichromat_2.0-1 withr_3.0.3 scales_1.4.0
#> [52] rmarkdown_2.32 XVector_0.53.0 otel_0.2.0
#> [55] evaluate_1.0.5 knitr_1.52 viridisLite_0.4.3
#> [58] irlba_2.4.1 rlang_1.3.0 Rcpp_1.1.2
#> [61] mixsqp_0.3-54 glue_1.8.1 BiocManager_1.30.27
#> [64] jsonlite_2.0.0 plyr_1.8.9 mashr_0.2.79
#> [67] R6_2.6.1