---
title: "Advanced Spatial Analysis"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Advanced Spatial Analysis}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
resource_files:
- figures/logo.svg
- figures/logo-dark.svg
---
[](https://github.com/CTTIR/phenoscapR/actions/workflows/R-CMD-check.yaml)
[](https://cttir.github.io/phenoscapR/)
[](https://app.codecov.io/gh/CTTIR/phenoscapR?branch=main)
[](https://opensource.org/licenses/MIT)
[](https://lifecycle.r-lib.org/articles/stages.html#experimental)
```{r setup, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.width = 7,
fig.height = 4.5,
dpi = 96
)
has_deldir <- requireNamespace("deldir", quietly = TRUE)
```
## Overview
Beyond nearest neighbours and density, phenoscapR provides a suite of point-
pattern statistics for asking *how* cells are spatially organised:
- **Neighbourhood enrichment** — do phenotype pairs co-localise more than chance?
- **Ripley's K / L** — clustering or dispersion across many scales at once
- **Pair correlation function** — clustering at a *specific* distance
- **Moran's I** — spatial autocorrelation of a continuous marker
- **Quadrat analysis** — a quick chi-squared test of spatial randomness
- **Cross nearest-neighbour distance** — directed proximity between two phenotypes
- **Delaunay network** — the cell-contact graph
- **Expression clustering** — marker-defined communities, ignoring position
### One sample at a time
These statistics describe a single point pattern in one observation window.
Pooling several tissues would compute distances across samples and return
nonsense, so phenoscapR **refuses multi-sample objects** for these functions —
subset to one sample first. (Per-cell summaries like `FindNeighbours()`,
`CellDensity()`, and `InteractionMatrix()` are sample-aware and need no
subsetting.)
```{r setup-data}
library(phenoscapR)
data(phenoscapR_example)
spe <- phenoscapR_example
spe@meta_data$phenotype <- spe@meta_data$phenotype_true
# A single tissue section.
a <- spe[spe$sample_id == "tonsil_A", ]
a
```
## Neighbourhood enrichment
A permutation test: it shuffles phenotype labels many times and asks whether
each phenotype pair is observed as neighbours more (positive z) or less
(negative z) often than the shuffles. The planted B-cell follicle shows up as
strong B–B self-enrichment.
```{r enrichment}
ne <- NeighbourhoodEnrichment(a, radius = 40, n_perm = 199, seed = 1)
ne[order(-ne$z_score), ][1:6, ]
```
## Ripley's K and L
`RipleysK()` accumulates, for each radius `r`, the average number of other
cells within `r`, normalised by intensity. Under complete spatial randomness
`K(r) = πr²`; the `L` transform makes departures easy to read (`L(r) > 0`
means clustering). The border correction removes edge bias.
```{r ripleysk, fig.alt = "Border-corrected Ripley L across radius for synthetic tonsil A, with the zero reference line."}
rk <- RipleysK(a, r_seq = seq(0, 120, length.out = 40), correction = "border")
head(rk)
with(rk, {
plot(r, L, type = "l", lwd = 2, col = "#3366cc",
xlab = "r", ylab = "L(r) (>0 = clustered)",
main = "Ripley's L, border-corrected")
abline(h = 0, lty = 2)
})
```
## Pair correlation function
`PairCorrelation()` is the scale-resolved companion to Ripley's K: `g(r) > 1`
means more pairs at distance `r` than expected, `g(r) < 1` means inhibition.
```{r pcf, fig.alt = "Pair correlation g across radius for synthetic tonsil A."}
pcf <- PairCorrelation(a, r_seq = seq(5, 120, length.out = 40))
with(pcf, {
plot(r, g, type = "l", lwd = 2, col = "#cc3366",
xlab = "r", ylab = "g(r)", main = "Pair correlation function")
abline(h = 1, lty = 2)
})
```
## Moran's I
Spatial autocorrelation of a continuous variable. A lineage marker like CD20,
concentrated in the follicle, autocorrelates strongly (`I` well above its
expectation, tiny p-value).
```{r moransi}
MoransI(a, feature = "CD20", radius = 50)
```
## Quadrat analysis
A fast omnibus test: tile the window into an `nx × ny` grid, count cells per
tile, and compare to a uniform expectation with a chi-squared test. The
variance-to-mean ratio (`VMR`) exceeds 1 for clustered patterns.
```{r quadrat}
qa <- QuadratAnalysis(a, nx = 5, ny = 5)
qa[c("chi_sq", "p_value", "VMR")]
```
## Cross nearest-neighbour distance
For each cell of phenotype `from`, the distance to the nearest cell of
phenotype `to` — a directed measure of how close one population sits to
another. It returns one distance per source cell.
```{r crossnn}
d_help_to_b <- CrossNNDistance(a, from = "T helper", to = "B cell")
summary(d_help_to_b)
```
## Delaunay contact network
The Delaunay triangulation approximates which cells physically touch. Capping
edge length with `max_edge` drops the long edges that span empty space.
```{r delaunay, eval = has_deldir, fig.alt = "Delaunay connections between synthetic tonsil A cells after excluding edges longer than 60."}
a <- DelaunayNetwork(a, max_edge = 60)
nrow(a@spatial$delaunay_edges)
SpatialNetworkPlot(a)
```
```{r delaunay-note, echo = FALSE, results = "asis", eval = !has_deldir}
cat("> The `deldir` package is not installed, so the Delaunay example is",
"skipped. Install it with `install.packages(\"deldir\")`.\n")
```
## Expression clustering
Cluster cells by their normalised marker profiles, independent of position,
then look at where those clusters fall in the tissue.
```{r exprclusters, fig.alt = "Synthetic tonsil A cell positions coloured by expression-derived cluster."}
a <- ExpressionClusters(a, k = 6)
table(a$expr_cluster)
CellMap(a, colour_by = "expr_cluster")
```
## Choosing a statistic
| Question | Use | Output |
|---|---|---|
| Do phenotypes A and B sit together? | `NeighbourhoodEnrichment()` | z-score / p per pair |
| Is one phenotype clustered overall? | `RipleysK()` (target =) | K, L vs radius |
| At *what distance* is it clustered? | `PairCorrelation()` | g(r) |
| Is a marker spatially autocorrelated? | `MoransI()` | I, z, p |
| Quick "is it random?" check | `QuadratAnalysis()` | chi-sq, VMR |
| How near is A to B? | `CrossNNDistance()` | distances |
| Which cells touch? | `DelaunayNetwork()` | edge list |
| Marker-defined communities | `ExpressionClusters()` | cluster labels |
All of these run on a kd-tree search engine, so they scale to large sections;
a 20,000-cell tissue completes each statistic in seconds.
## Next steps
- `?NeighbourhoodEnrichment`, `?RipleysK`, … for full parameter documentation
- Combine with `InteractionMatrix()` / `InteractionPlot()` for the co-localisation picture
- **[Visualisation Gallery](phenoscapR-04-visualisation.html)** for every plot
