SIB Spatial Omics Data Analysis Course 2026: Point Pattern Analysis Vignette (Part 1)

Author

Martin Emons, Samuel Gunz, Mark D. Robinson, Peiying Cai

Introduction

Here, we assume that cells (or transcripts) can be approximated as points given their location (ignoring for instance that cells can have a shape).

`spatstat` package and `ppp` object

The powerful package to analyse point patterns in R is called spatstat (Baddeley and Turner 2005). The main data object to compute on is called a ppp object. ppp objects describe point patterns in two dimensional space, ppx objects create multidimensional point patterns. A ppp object is made up of three specifications (Baddeley and Turner 2005):

The locations of the points in question (\(x\),\(y\) and, optionally, \(z\) coordinates)
The observation window
“Marks” that are associated to each point in the pattern

On this object, various spatstat metrics can be calculated.

`SpatialExperiment` Object

Often, the starting point in spatial omics data analysis in R is a SpatialExperiment (or similar) object. This is a central data structure in the Bioconductor framework to store spatial omics data. The data we consider here is a MERFISH assay of a mouse preoptic hypothalamus (Chen et al. 2015; Moffitt et al. 2018).

if(!file.exists("/home/rstudio/.cache/R/ExperimentHub")){
  dir.create("/home/rstudio/.cache/R/ExperimentHub", recursive = TRUE)
}

suppressPackageStartupMessages({
  
  library(SpatialExperiment)
  library(spatstat.explore)
  library(spatialFDA)

  library(ggplot2)
  library(dplyr)
  library(tidyr)
  library(magrittr)
  library(purrr)
  
})

# load the data from ExperimentHub
#basedir <- here::here("day_3/practical_5/")

basedir <- "/home/rstudio/work/practical_7/"

if(!file.exists(file.path(basedir, "data","brain.rds"))){
  
  load_and_save = function(basedir){
    
    library(ExperimentHub)
    
    #### Code from Helena Crowell ####
    eh <- ExperimentHub()
    # query(eh, "MERFISH")
    df <- eh[["EH7546"]]
    
    # extract cell metadata
    i <- seq_len(9)
    cd <- data.frame(df[, i], row.names = 1)
    
    # set sample identifiers
    names(cd)[grep("Bregma", names(cd))] <- "sample_id"
    
    # rename spatial coordinates
    names(cd)[grep("Centroid", names(cd))] <- c("x", "y")
    
    # simplify annotations
    cd$cluster_id <- cd$Cell_class
    for (. in c("Endothelial", "OD Mature", "OD Immature"))
      cd$cluster_id[grep(., cd$cluster_id)] <- .
    
    # extract & sparsify assay data
    y <- data.frame(df[, -i], row.names = df[, 1])
    y <- as(t(as.matrix(y)), "dgCMatrix")
    
    # construct SPE
    (spe <- SpatialExperiment(
      assays = list(exprs  = y),
      spatialCoordsNames = c("x","y"),
      colData = cd))
    
    # Define the directory and file paths
    dir_path <- paste0(basedir, "/data")
    file_path <- file.path(dir_path, "brain.rds")
    
    # Check if the directory exists, and create it if it doesn't
    if (!dir.exists(dir_path)) {
      dir.create(dir_path)
    }
    
    #save data as rds file
    saveRDS(spe,  file_path)
    
  }
  
  load_and_save(basedir)
  
}

theme_set(theme_light())
#spe <- readRDS("data/brain.rds")
spe <- readRDS(file.path(basedir, "data","brain.rds"))
spe

class: SpatialExperiment 
dim: 161 73655 
metadata(0):
assays(1): exprs
rownames(161): Ace2 Adora2a ... Ucn3 Vgf
rowData names(0):
colnames(73655): 6749ccb4-2ed1-4029-968f-820a287f43c8
  6cac74bd-4ea7-4701-8701-42563cc65eb8 ...
  6b666f81-7b73-4100-9e02-b5381b39f0f3
  fdcddd97-7701-462a-b48f-979111245bd5
colData names(7): Animal_ID Animal_sex ... Neuron_cluster_ID cluster_id
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):
spatialCoords names(2) : x y
imgData names(0):

invisible(gc())

We see that we have an object of class SpatialExperiment with \(161\) genes (rows) and \(73655\) cells. This object is very similar to a SingleCellExperiment object except it has the added spatialCoords slot. One slot in the colData is called sample_id which defines the so called z-slices. The three dimensional tissue is cut in the z-axis into consecutive two dimensional slices (Righelli et al. 2022).

Just to remind users of typical accessor functions for SpatialExperiment objects, here are a few useful accessor functions:

dim(spe)

[1]   161 73655

spatialCoords(spe) %>% head

                                             x        y
6749ccb4-2ed1-4029-968f-820a287f43c8 -3211.562 2608.541
6cac74bd-4ea7-4701-8701-42563cc65eb8 -3207.923 2621.795
9f29bd57-16a5-4b26-b9f5-37598809da9e -3209.578 2633.153
d7eb4e0b-276e-47e3-a55c-0b033180a2fe -3203.854 2756.046
54434f3a-eba9-4aec-af35-c9d317ffa1d5 -3202.683 2608.804
6dcfb0ba-763f-4330-8ff6-aa55fd5d6ba6 -3203.466 2743.163

colData(spe) %>% head(2)

DataFrame with 2 rows and 7 columns
                                     Animal_ID  Animal_sex    Behavior
                                     <numeric> <character> <character>
6749ccb4-2ed1-4029-968f-820a287f43c8         1      Female       Naive
6cac74bd-4ea7-4701-8701-42563cc65eb8         1      Female       Naive
                                       sample_id  Cell_class Neuron_cluster_ID
                                     <character> <character>       <character>
6749ccb4-2ed1-4029-968f-820a287f43c8        0.26   Astrocyte                NA
6cac74bd-4ea7-4701-8701-42563cc65eb8        0.26  Inhibitory               I-5
                                      cluster_id
                                     <character>
6749ccb4-2ed1-4029-968f-820a287f43c8   Astrocyte
6cac74bd-4ea7-4701-8701-42563cc65eb8  Inhibitory

colnames(colData(spe))

[1] "Animal_ID"         "Animal_sex"        "Behavior"         
[4] "sample_id"         "Cell_class"        "Neuron_cluster_ID"
[7] "cluster_id"

rowData(spe)

DataFrame with 161 rows and 0 columns

assays(spe)

List of length 1
names(1): exprs

dim(assay(spe,1))

[1]   161 73655

assay(spe,1) %>% class

[1] "dgCMatrix"
attr(,"package")
[1] "Matrix"

assay(spe,1)[,1:3] %>% head(3)

3 x 3 sparse Matrix of class "dgCMatrix"
        6749ccb4-2ed1-4029-968f-820a287f43c8
Ace2                                .       
Adora2a                             1.638275
Aldh1l1                            21.299750
        6cac74bd-4ea7-4701-8701-42563cc65eb8
Ace2                                .       
Adora2a                             .       
Aldh1l1                             1.578873
        9f29bd57-16a5-4b26-b9f5-37598809da9e
Ace2                                .       
Adora2a                             .       
Aldh1l1                             2.701349

spe$sample_id %>% head

[1] "0.26" "0.26" "0.26" "0.26" "0.26" "0.26"

table(spe$sample_id)


-0.04 -0.09 -0.14 -0.19 -0.24 -0.29  0.01  0.06  0.11  0.16  0.21  0.26 
 6154  6185  6605  6507  6412  6509  6111  6144  5799  6067  5578  5584

EXERCISE

How many ‘z’ slices does the dataset have?
How many cells are annotated as ‘OD Mature’ in each of the 3 slices used below (-0.09, 0.01, 0.21)?

Solution

# Number of z-slices 
length(unique(spe$sample_id))

[1] 12

# Number of OD Mature cells in each selected slice
selected_slices <- c("-0.09", "0.01", "0.21")
cell_counts <- table(spe$sample_id, spe$cluster_id)
cell_counts[selected_slices, "OD Mature"]

-0.09  0.01  0.21 
  240   457  1136

Next, we want to extract three 2D slices of this SpatialExperiment object and convert them into ppp objects. The spatialFDA package contains a function for this conversion: spatialFDA::.dfToppp.

# define the Z-stacks that you want to compare
zstack_list <- list("-0.09", "0.01", "0.21")
# small helper function to extract the z-slices and convert them to `ppp` objects
selectZstacks <- function(zstack, spe) {
  df <- spatialFDA::.speToDf(spe[, spe$sample_id == zstack])
  spatialFDA::.dfToppp(df, marks = "cluster_id")
}
(pp_ls <- lapply(zstack_list, selectZstacks, spe) |>
  setNames(zstack_list))

$`-0.09`
Marked planar point pattern: 6185 points
Multitype, with levels = 
   Ambiguous Astrocyte Endothelial Ependymal Excitatory Inhibitory Microglia OD 
Immature OD Mature Pericytes
window: rectangle = [1306.7716, 3098.8543] x [2113.368, 3904.386] units

$`0.01`
Marked planar point pattern: 6111 points
Multitype, with levels = 
   Ambiguous Astrocyte Endothelial Ependymal Excitatory Inhibitory Microglia OD 
Immature OD Mature Pericytes
window: rectangle = [1221.5635, 3013.4248] x [-3994.535, -2201.755] units

$`0.21`
Marked planar point pattern: 5578 points
Multitype, with levels = 
   Ambiguous Astrocyte Endothelial Ependymal Excitatory Inhibitory Microglia OD 
Immature OD Mature Pericytes
window: rectangle = [1327.6103, 3125.2492] x [2198.257, 3994.499] units

We see that we obtain a list of three ppp objects for the three z-slices \(-0.09, 0.01, 0.21\).

We can plot one of these slices, e.g. slice \(-0.09\) with ggplot

# create a dataframe from ppp and plot with ggplot
slice_df <- as.data.frame(pp_ls[["-0.09"]])
ggplot(slice_df, aes(x, y, colour = marks)) +
  geom_point(size = .5) +
  coord_equal()

EXERCISE

Modify the ggplot to highlight OD Mature cells

Solution

ggplot(slice_df, aes(x, y, colour = marks)) +
  geom_point(size = .5) +
  geom_point(
    data = subset(slice_df, marks == "OD Mature"),
    aes(x, y),
    shape = 1,
    colour = "black"
  ) +
  coord_equal()

Windows

One important aspect of a point pattern is the observation window, which represents the region in which a pattern is observed or, e.g., a survey was conducted (Baddeley et al. 2015, 85). In most microscopy use cases, we encounter so-called window sampling. Window sampling describes the case where we don’t observe the entire point pattern in a window but just a sample (Baddeley et al. 2015, 143–45).

The window of a point pattern does not need to be rectangular; we can receive round biopsies or calculate convex hulls around our sample (Baddeley et al. 2015, 143–45).

Let’s investigate the observation window for the slice \(-0.09\).

# subset point pattern list
pp_sub <- pp_ls$`-0.09`
# base R plot of all marks
pp_sub |> plot()

Window(pp_sub)

window: rectangle = [1306.7716, 3098.8543] x [2113.368, 3904.386] units

Here, we have a rectangular window around all points.

Let’s investigate what a round window would look like:

pp_sub_round <- pp_sub
# calculate circle with radius 850 µm and a center at the centroid of the window would look like
w <- disc(r = 850, centroid.owin(Window(pp_sub)))
Window(pp_sub_round) <- w
pp_sub_round |> plot()

Correctly assigning windows is important. The window should represent the space where points are expected. This means, in window sampling, one should not restrict the window. This would lead to a false underestimation of the area where the points can be potentially observed. This problem of where we can observe points and where not (beyond the boundary of the window) leads to a range of problems collectively called edge effects (Baddeley et al. 2015, 143–45). A short discussion of this is given below.

Marks

The next concept that defines a point pattern is that marks can be associated with the points. Point patterns can even be unmarked (e.g., unmark(pp_sub)). In the context of cell biology, we can distinguish between discrete marks (e.g., cell types) or continuous marks (e.g., gene expression).

Discrete Marks

In our example, we have a multitype point pattern, meaning there are different cell types that serve as marks for the point pattern. Multitype means that we consider all marks together. An alternative is multivariate, where we consider the marks independently (Baddeley et al. 2015, 564 ff.).

Compare the multitype case:

pp_sub |> plot()

with the multivariate view on the same pattern:

pp_sub |>
  split() |>
  plot()

Continuous Marks

Marks can also be continuous as in the case of gene expression. We choose some genes from the original paper and look at their distribution (Baddeley et al. 2015, 637; Moffitt et al. 2018).

# subset the SpatialExperiment to our example slice -0.09
sub <- spe[, spe$sample_id == "-0.09"]

#  Genes from Fig. 6 of Moffitt et al. (2018)
genes <- c("Slc18a2", "Esr1", "Pgr")

gex <- assay(sub)[genes, ] %>%
  t() %>%
  as.matrix() %>%
  data.frame() %>%
  set_rownames(NULL)

# gene expression to marks
marks(pp_sub) <- gex

Now that we have points with multivariate continuous marks

# create a dataframe in long format for plotting
pp_df <- pp_sub %>%
  as.data.frame() %>%
  pivot_longer(cols = 3:5)

ggplot(pp_df, aes(x, y, colour = log(value + 1))) +
  geom_point(size = .5) +
  facet_wrap(~name) +
  coord_equal() +
  scale_color_continuous(type = "viridis")

EXERCISE

Adjust the code blocks above to also visualize the expression of ‘Prlr’

Solution

Make the following two changes:

genes <- c("Slc18a2", "Esr1", "Pgr", "Prlr")

pivot_longer(cols = all_of(genes))

Then re-run the modified code blocks.

Within Mark Comparison

We can compare patterns between marks of the same type, which is referred to as a within-mark comparison in our vignette. We can compare discrete marks, i.e., the distribution of one single mark (e.g., a cell type).

Here, we plot the distribution of mature oligodendrocytes across three slices of one 3D brain sample.

# create a dataframe from the point pattern
pp_df_discrete <- lapply(zstack_list, function(x) {
  df <- pp_ls[[x]] %>% as.data.frame()
  df$stack <- x
  return(df)
}) %>% bind_rows()

# select OD Mature cells
pp_df_odmature <- pp_df_discrete[pp_df_discrete$marks == "OD Mature", ]

ggplot(pp_df_odmature, aes(x, y, colour = marks)) +
  geom_point(size = 0.5) +
  facet_wrap(~stack, scales = "free") +
  theme(aspect.ratio = 1)

Continuous marks can be compared as well, e.g. the expression of a gene across slices of a tissue.

Correlation

Correlation is a second order quantity that measures the dependence between points (Baddeley et al. 2015, 199). A famous way to measure this is with Ripley’s \(K\), which is a cumulative function that quantifies the “number of \(r\)-neighbours of a typical random point” (Baddeley et al. 2015, 204; Ripley 1976).

Global Measures

Global correlation measures quantify the correlation in the entire window. Global Ripley’s \(K\) is defined as:

\[\hat{K}(r) = \frac{|W|}{n(n-1)}\sum_{i=1}^n\sum_{\substack{j=1 \\ j \neq i}}^n\mathbb{1}\{d_{ij}\leq r\} e_{ij}(r)\]

In the formula above we note a few things:

The function is normalised by the number of points \(n\) and the window size \(|W|\)
the term \(e_{ij}(r)\) is an edge correction, which is covered briefly below.

Ripley’s \(K\) function can be variance-stabilized, which is referred to as Besag’s \(L\) (Canete et al. 2022; Besag 1977). The idea behind variance stabilisation is to “uncouple” the relationship between mean and variance. By taking the square root of the function in question, the variance is nearly constant across the function (Bartlett 1947).

\[L(r) = \sqrt{\frac{K(r)}{\pi}}\]

m <- calcMetricPerFov(spe, selection = 'OD Mature',
                      subsetby = 'sample_id',
                      fun = 'Kest',
                      marks = 'cluster_id',
                      by = c('sample_id'))

head(m, 2)

p <- plotMetricPerFov(m %>% subset(sample_id %in% c('-0.09', 
                                                  '0.01', '0.21')), 
                      correction = "border", 
                      theo = TRUE, x = "r", 
                      imageId = 'sample_id', 
                      legend.position = "right",
                      linewidth = 3) +
  labs(colour = "Slice")
p

The strongest estimate of association between oligodendrocytes is found for slices \(0.01\). Slice \(0.21\) does not show such a high degree of association at radii \(\leq300\) as the other two slices. This means that the apparent clustering we see in the distribution of points is mainly due to an overall higher number of cells in slice \(0.21\) and not a higher degree of association per se. The black line indicates the expected \(K\) function under completely spatially random Poisson process (Baddeley et al. 2015, 132 ff.).

EXERCISE

Choose another summary function from spatstat.explore (see the Summary statistics for a point pattern section in ?spatstat.explore).

For example, use Besag’s \(L\) (Lest) to compare the spatial arrangement of OD Mature cells across the same three sections. How do the resulting curves compare with the Ripley’s \(K\) curves above?

Solution

calcMetricPerFov(
  spe, selection = "OD Mature",
  subsetby = "sample_id",
  fun = "Lest",
  marks = "cluster_id",
  by = "sample_id"
) |>
  dplyr::filter(sample_id %in% c("-0.09", "0.01", "0.21")) |>
  plotMetricPerFov(
    correction = "border",
    theo = TRUE, x = "r",
    imageId = "sample_id",
    linewidth = 3
  )

A similar analysis can be performed for continuous marks using a mark-weighted correlation function (markcorr). For a continuous mark such as gene expression, it evaluates whether pairs of cells separated by distance (r) tend to have jointly high or low mark values. The mark weighted correlation function is defined as:

\[k_f(r) = \frac{\mathbb{E}[f(m(u),m(v))|u,v \in X]}{\mathbb{E}[f(M,M')]}\]

where the numerator is the conditional expectation of the marks at location \(u,v\) separated by a radius \(r\) and \(f\) can be any function linking the two marks. The denominator is the expectation of two random marks \(M,M'\) (Baddeley et al. 2015, 603).

EXERCISE

Calculate the mark correlation function for the expression of Slc18a2, Esr1, and Pgr in section -0.09

Hint: modify subsetby, fun, and marks, and set continuous = TRUE. Use ?markcorr for more information.

Solution

genes <- c('Slc18a2', 'Esr1', 'Pgr')

calcMetricPerFov(
  spe[,spe$sample_id == "-0.09"],
  selection = "all",
  subsetby = genes,
  fun = 'markcorr',
  marks = genes,
  by = 'sample_id',
  continuous = TRUE
) |> 
  plotMetricPerFov(
    correction = "iso",
    theo = TRUE,
    x = "r",
    imageId = 'gene', linewidth = 3
    ) + 
  theme(legend.position = "bottom")

invisible(gc())

At short distances, all three genes show positive spatial association, with Slc18a2 showing the strongest association. This indicates that nearby cells tend to have jointly high expression values.

Local Measures

Next to observation window metrics, we can calculate point level statistics as well. One such option is the local indicators of spatial association (LISA). This gives one curve per point in the field of view (Baddeley et al. 2015, 247–48).

pp <- subset(pp_ls[["0.01"]], marks %in% "OD Mature")
df <- localL(pp) %>% as.data.frame

# data.frame manipulation for plotting
dfm <- pivot_longer(df,  cols = colnames(df)[colnames(df) != "r"], names_to = "variable")
get_sel <- dfm %>%
  filter(r > 200.5630 & r < 201.4388, variable != "theo") %>%
  mutate(sel = value) %>%
  select(variable, sel)
dfm <- dfm %>% left_join(get_sel)

ggplot(dfm, aes(x = r, y = value, 
                     group = variable, colour = sel)) +
  geom_line(linewidth = 1) +
  scale_color_continuous(type = "viridis") +
  geom_vline(xintercept = 200) +
  theme(legend.position = "bottom") +
  ggtitle("LISA curves of slice 0.01")

It may be useful to compress these curves into something easier to interpret, such as plotting the first functional principal components from a functional PCA (Baddeley et al. 2015, 247–48; Ramsay and Silverman 2005).

# extract the functional response matrix
mat <- df %>% select(!c(r, theo)) %>% t()

dat <- data.frame(ID = rownames(mat))
dat$Y <- mat
dat$sel <- get_sel$sel

# perform functional PCA
res <- functionalPCA(dat, unique(df$r))

ggplot(res$scores %>% as.data.frame(),
       aes(V1, V2, colour = (dat[['sel']]))) +
  geom_point() +
  coord_equal() +
  theme_light() +
  scale_color_continuous(type = "viridis") +
  xlab('PC1') +
  ylab('PC2')

EXERCISE

Play with the scale of the curves (e.g., although we’ve used L, which is a variance-stabilized transform of the K function, it may not have fully resolved the mean-variance relationship); for example, consider sqrt() or asinh(. / 50) or similar. Recompute and plot the functional PC embeddings.
Do the exploratory analysis the other way around; for each point, plot the value of the functional PC1 in the spatial context (above, we coloured by the value of the local L-curve at r=~200).

Solution

# Exercise 1
# transform the local L curves
df_trans <- df %>%
  mutate(across(!c(r, theo), ~ asinh(.x / 50)))

# Exercise 2
pp_df <- as.data.frame(pp)
pp_df$fPC1 <- res$scores[, 1]

ggplot(pp_df, aes(x, y, colour = fPC1)) +
  geom_point(size = 2) +
  coord_equal() +
  scale_color_continuous(type = "viridis") +
  theme_bw() +
  labs(title = "Functional PC1 in spatial context",
       colour = "fPC1")

invisible(gc())

Cross Mark Comparison

Similar to above, analyses can be performed between two cell types. The corresponding functions are called cross functions (Baddeley et al. 2015, 594 ff.). As an exercise, try to implement (similar to the analyses above) a cross comparison between two cell types of interest. With the provided functions this is possible, just give a function you like as input and a vector with two cell types you wish to compare. You can look at functions in the documentation of spatstat.explore. We provide an example below:

calcCrossMetricPerFov(
  spe,
  selection = c("OD Mature", "Microglia"),
  subsetby = 'sample_id',
  fun = 'Lcross.inhom',
  marks = 'cluster_id',
  rSeq = seq(0, 500, length.out = 100),
  by = c('sample_id')
) |>
plotCrossMetricPerFov(
  theo = TRUE,
  correction = "iso",
  x = "r",
  imageId = 'sample_id', 
  linewidth = 3
) |>
  pluck(2)

We note that there is not a very strong co-localisation indicated by the \(L\) curves between mature oligodendrocytes and microglia cells. If we look at their spatial distribution that makes sense since microglia cells are distributed more or less homogeneously in the respective slices.

EXERCISE

Compare the colocalisation of OD Mature and Ependymal cells in section 0.01.

Choose an appropriate cross-type summary function from spatstat.explore. See the Summary statistics for a multitype point pattern section in ?spatstat.explore for available functions. calcMetricPerFov is a helper (wrapper) function that is useful to operate on SpatialExperiment objects.

Solution

Here, we use the inhomogeneous cross-type (L)-function:

calcMetricPerFov(
  spe[, spe$sample_id == "0.01"],
  selection = c("OD Mature", "Ependymal"),
  subsetby = "sample_id",
  fun = "Lcross.inhom",
  marks = "cluster_id",
  rSeq = seq(0, 500, length.out = 100),
  by = "sample_id",
  correction = "isotropic"
) |>
  plotMetricPerFov(
    theo = TRUE,
    correction = "iso",
    x = "r",
    imageId = "sample_id",
    linewidth = 3
  )

EXERCISE (advanced!)

Above, you computed functional PCA of a set of local L curves (i.e., one curve for each cell). Perhaps more interesting is a curve for each field of view (FOV) across a larger set of FOVs and samples. Find a dataset of your choice (e.g., imcdatasets), compute a spatial statistics curve (of choice) for each FOV across samples. Explore the spatial heterogeneity across your dataset.

Edge effects and their corrections for spatial metrics

Edge effects describe the phenomenon that not the entire point process is observed, but rather only the part within the window \(W\). This means the value of various statistics could be biased along the edges (Baddeley et al. 2015, 213).

There are many corrections for edge effects that are briefly listed here (Baddeley et al. 2015, 214–19):

Border correction:

In border correction the summation of data points is restricted to \(x_i\) for which \(b(x_i,r)\) is completely in the window \(W\).

Isotropic correction:

We can regard edge effect as a sampling bias. Larger distances (e.g. close to the edges) are less likely to be observed. This can be corrected for.

Translation correction:

A stationary point process \(X\) is invariant to translations. So the entire point process can be shifted by a vector \(s\) to be at the position \(X+s\).

Summary and Considerations

Point patterns are realisations of a point process. In the analysis, we make inferences about the point process.
A point process assumes stochasticity. Therefore, low-resolution HTS-based approaches are not suitable for point pattern analysis.
There are global metrics for the comparison within a cell type or between cell types.
There are corresponding metrics for single cells and their interactions.
Point pattern analysis allows for the analysis of continuous gene expression marks as well.

sessionInfo()

R version 4.5.3 (2026-03-11)
Platform: x86_64-conda-linux-gnu
Running under: Ubuntu 26.04 LTS

Matrix products: default
BLAS/LAPACK: /app/.pixi/envs/default/lib/libopenblasp-r0.3.33.so;  LAPACK version 3.12.0

locale:
 [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
 [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
 [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
[10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   

time zone: Etc/UTC
tzcode source: system (glibc)

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] purrr_1.2.2                 magrittr_2.0.5             
 [3] tidyr_1.3.2                 dplyr_1.2.1                
 [5] ggplot2_4.0.3               spatialFDA_1.2.1           
 [7] spatstat.explore_3.8-0      nlme_3.1-169               
 [9] spatstat.random_3.4-5       spatstat.geom_3.8-1        
[11] spatstat.univar_3.2-0       spatstat.data_3.1-9        
[13] SpatialExperiment_1.20.0    SingleCellExperiment_1.32.0
[15] SummarizedExperiment_1.40.0 Biobase_2.70.0             
[17] GenomicRanges_1.62.1        Seqinfo_1.0.0              
[19] IRanges_2.44.0              S4Vectors_0.48.0           
[21] BiocGenerics_0.56.0         generics_0.1.4             
[23] MatrixGenerics_1.22.0       matrixStats_1.5.0          

loaded via a namespace (and not attached):
  [1] Rdpack_2.6.6          bitops_1.0-9          DBI_1.3.0            
  [4] deldir_2.0-4          httr2_1.2.2           rlang_1.2.0          
  [7] compiler_4.5.3        RSQLite_3.53.1        mgcv_1.9-4           
 [10] png_0.1-9             vctrs_0.7.3           pkgconfig_2.0.3      
 [13] crayon_1.5.3          fastmap_1.2.0         dbplyr_2.6.0         
 [16] magick_2.9.1          XVector_0.50.0        magic_1.6-1          
 [19] labeling_0.4.3        deSolve_1.42          rmarkdown_2.32       
 [22] pracma_2.4.6          nloptr_2.2.1          bit_4.6.0            
 [25] xfun_0.57             cachem_1.1.0          refund_0.1-40        
 [28] jsonlite_2.0.0        goftest_1.2-3         blob_1.3.0           
 [31] DelayedArray_0.36.1   spatstat.utils_3.2-3  cluster_2.1.8.2      
 [34] parallel_4.5.3        R6_2.6.1              RColorBrewer_1.1-3   
 [37] boot_1.3-32           rainbow_3.8           Rcpp_1.1.1-1.1       
 [40] knitr_1.51            tensor_1.5.1          hdrcde_3.5.0         
 [43] Matrix_1.7-5          splines_4.5.3         tidyselect_1.2.1     
 [46] abind_1.4-8           yaml_2.3.12           RLRsim_3.1-9         
 [49] curl_7.1.0            lattice_0.22-9        tibble_3.3.1         
 [52] withr_3.0.3           ks_1.15.2             KEGGREST_1.50.0      
 [55] S7_0.2.2              pbs_1.1               evaluate_1.0.5       
 [58] grpreg_3.6.0          polyclip_1.10-7       BiocFileCache_3.0.0  
 [61] ExperimentHub_3.0.0   mclust_6.1.2          Biostrings_2.78.0    
 [64] pillar_1.11.1         BiocManager_1.30.27   fda_6.3.0            
 [67] filelock_1.0.3        KernSmooth_2.23-26    reformulas_0.4.4     
 [70] pcaPP_2.0-5           RCurl_1.98-1.19       BiocVersion_3.22.0   
 [73] scales_1.4.0          minqa_1.2.8           BiocStyle_2.38.0     
 [76] gamm4_0.3-0           glue_1.8.1            tools_4.5.3          
 [79] AnnotationHub_4.0.0   fds_1.9               lme4_2.0-1           
 [82] mvtnorm_1.4-1         grid_4.5.3            rbibutils_2.4.1      
 [85] colorspace_2.1-2      AnnotationDbi_1.72.0  patchwork_1.3.2      
 [88] cli_3.6.6             spatstat.sparse_3.2-0 rappdirs_0.3.4       
 [91] viridisLite_0.4.3     S4Arrays_1.10.1       gtable_0.3.6         
 [94] digest_0.6.39         SparseArray_1.10.8    rjson_0.2.23         
 [97] htmlwidgets_1.6.4     farver_2.1.2          memoise_2.0.1        
[100] htmltools_0.5.9       lifecycle_1.0.5       httr_1.4.9           
[103] bit64_4.8.2           MASS_7.3-65

References

Baddeley, Adrian, Ege Rubak, and Rolf Turner. 2015. Spatial Point Patterns. 1st ed. CRC Interdisciplinary Statistics Series. CRC Press, Taylor & Francis Group. https://doi.org/10.1201/b19708.

Baddeley, Adrian, and Rolf Turner. 2005. “Spatstat: An R Package for Analyzing Spatial Point Patterns.” Journal of Statistical Software 12 (January): 1–42. https://doi.org/10.18637/jss.v012.i06.

Bartlett, M. S. 1947. “The Use of Transformations.” Biometrics 3 (1): 39–52. https://doi.org/10.2307/3001536.

Besag, Julian. 1977. “Contribution to the Discussion on Dr Ripley’s Paper.” JR Stat Soc B 39: 193–95.

Canete, Nicolas P, Sourish S Iyengar, John T Ormerod, Heeva Baharlou, Andrew N Harman, and Ellis Patrick. 2022. “spicyR: Spatial Analysis of in Situ Cytometry Data in R.” Bioinformatics 38 (11): 3099–105. https://doi.org/10.1093/bioinformatics/btac268.

Chen, Kok Hao, Alistair N. Boettiger, Jeffrey R. Moffitt, Siyuan Wang, and Xiaowei Zhuang. 2015. “Spatially Resolved, Highly Multiplexed RNA Profiling in Single Cells.” Science 348 (6233): aaa6090. https://doi.org/10.1126/science.aaa6090.

Moffitt, Jeffrey R., Dhananjay Bambah-Mukku, Stephen W. Eichhorn, et al. 2018. “Molecular, Spatial, and Functional Single-Cell Profiling of the Hypothalamic Preoptic Region.” Science 362 (6416): eaau5324. https://doi.org/10.1126/science.aau5324.

Ramsay, JO, and BW Silverman. 2005. “Principal Components Analysis for Functional Data.” Functional Data Analysis, 147–72. https://doi.org/10.1007/0-387-22751-2_8.

Righelli, Dario, Lukas M Weber, Helena L Crowell, et al. 2022. “SpatialExperiment: Infrastructure for Spatially-Resolved Transcriptomics Data in R Using Bioconductor.” Bioinformatics 38 (11): 3128–31. https://doi.org/10.1093/bioinformatics/btac299.

Ripley, B. D. 1976. “The Second-Order Analysis of Stationary Point Processes.” Journal of Applied Probability 13 (2): 255–66. https://doi.org/10.2307/3212829.

Introduction

spatstat package and ppp object

SpatialExperiment Object

Windows

Marks

Discrete Marks

Continuous Marks

Within Mark Comparison

Correlation

Global Measures

Local Measures

Cross Mark Comparison

Edge effects and their corrections for spatial metrics

Summary and Considerations

References

`spatstat` package and `ppp` object

`SpatialExperiment` Object