SIB Spatial Omics Data Analysis Course 2026 (Part 2): Lattice Data Analysis Vignette

Author

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

Introduction

In contrast to point pattern based methods, we can view the location of cells or spots as a fixed lattice and measure the corresponding marker expression at each location. HTS-based spatial transcriptomic technologies often produce data on a regular lattice (i.e., approximately evenly spaced spots or beads of uniform size and shape), whereas imaging-based technologies yield irregular lattice structures (i.e., with variable cell sizes and shapes, and non-uniform spacing).

For this representation of the cells, we will rely on the SpatialFeatureExperiment package (see below for more details). For preprocessing of the dataset, we refer the reader to the vignette of the Voyager package (Moses et al. 2023). The Voyager package also provides wrapper functions around the package spdep (Pebesma and Bivand 2023) that work directly on the SpatialFeatureExperiment object. The package spdep is designed for the analysis of spatial data with lattice structure.

Dependencies

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

Sys.setenv(PROJ_DATA = "/app/.pixi/envs/default/share/proj") # Required for Voyager to work in RStudio within container

suppressPackageStartupMessages({
  library(Voyager)
  library(SpatialFeatureExperiment)
  library(SFEData)
  library(sf)
  library(scater)
})

theme_set(theme_light())

Setup and Preprocessing

We will load a dataset generated by (McKellar et al. 2021) using the Visium technology (Ståhl et al. 2016). The data shows a sample taken from the tibialis anterior muscle of a mouse.

# Load the dataset

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

if(!file.exists(file.path(basedir, "data","muscle.rds"))){
  
  sfe <- SFEData::McKellarMuscleData(dataset = "full")
  
  dir_path <- paste0(basedir, "/data")
  file_path <- file.path(dir_path, "muscle.rds")
  
  # Check if the directory exists, and create it if it doesn't
  if (!dir.exists(dir_path)) {
    dir.create(dir_path)
  }
  
  saveRDS(sfe, file.path(basedir, "data","muscle.rds"))
}
see ?SFEData and browseVignettes('SFEData') for documentation
downloading 1 resources
retrieving 1 resource
loading from cache
sfe = readRDS(file.path(basedir, "data","muscle.rds"))
# Take spots that are covered with tissue
sfe_tissue <- sfe[, colData(sfe)$in_tissue]
rm(sfe)
invisible(gc())
# Filter out genes with no counts
sfe_tissue <- sfe_tissue[rowSums(counts(sfe_tissue)) > 0, ]
# Convert counts to log-transformed normalized expression values
sfe_tissue <- scater::logNormCounts(sfe_tissue)

SpatialFeatureExperiment (Moses et al. 2023) objects are an extension of the SpatialExperiment object (Righelli et al. 2022). In particular, the objects contain geometric annotations that are encoded as simple features of the sf library (Pebesma and Bivand 2023).

sfe_tissue
class: SpatialFeatureExperiment 
dim: 15043 932 
metadata(0):
assays(2): counts logcounts
rownames(15043): ENSMUSG00000025902 ENSMUSG00000096126 ...
  ENSMUSG00000064368 ENSMUSG00000064370
rowData names(6): Ensembl symbol ... vars cv2
colnames(932): AAACATTTCCCGGATT AAACCTAAGCAGCCGG ... TTGTGTTTCCCGAAAG
  TTGTTGTGTGTCAAGA
colData names(13): barcode col ... in_tissue sizeFactor
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):
spatialCoords names(2) : imageX imageY
imgData names(1): sample_id

unit: full_res_image_pixels
Geometries:
colGeometries: spotPoly (POLYGON) 
annotGeometries: tissueBoundary (POLYGON), myofiber_full (POLYGON), myofiber_simplified (POLYGON), nuclei (POLYGON), nuclei_centroid (POINT) 

Graphs:
Vis5A:

For example, the spots of the Visium dataset are stored as a simple feature collection.

colGeometry(sfe_tissue, "spotPoly")
Simple feature collection with 932 features and 0 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 4903.08 ymin: 3644.16 xmax: 12150.06 ymax: 14931.74
CRS:           NA
First 10 features:
                         geometry
1  POLYGON ((6037.973 5430.826...
2  POLYGON ((7977.344 4461.467...
3  POLYGON ((8547.633 11699.1,...
4  POLYGON ((6867.944 4463.856...
5  POLYGON ((9779.598 4216.29,...
6  POLYGON ((10206.53 9282.577...
7  POLYGON ((6462.832 9531.935...
8  POLYGON ((6186.521 10015.12...
9  POLYGON ((9928.147 8800.585...
10 POLYGON ((8538.28 7355.805,...

Lattice Data

Definition

Lattice data refers to spatial data collected at locations arranged in a regular or irregular grid (lattice). Each location has a defined spatial unit, and the sampling locations are fixed rather than random. This approach contrasts with point pattern analysis, where we assume that the locations were generated by a stochastic process (Zuur et al. 2007; Pebesma and Bivand 2023).

As can be seen below, the Visium dataset is a regular (hexagonal) lattice. The color shows the number of counts detected in each Visium spot. In contrast, the outlines of individual cells after segmentation (e.g., from a higher resolution spatial omics assay, or an H&E stain) could be thought of and represented as an irregular lattice. As this dataset also contains (manual) segmentations of myofibers (muscle cells stored as myofiber_simplified), we will illustrate calculations based on irregular lattice on the myofiber segmentations.

# A plot using the plotSpatialFeature from Voyager
p <- plotSpatialFeature(sfe_tissue,
  "nCounts",
  colGeometryName = "spotPoly"
)

# Extract the segmented myofiber geometries
myofibers <- annotGeometry(sfe_tissue, "myofiber_simplified") |>
  st_geometry()

# We can also use ggplot and geom_sf to plot sf objects
q <- ggplot() +
  geom_sf(data = myofibers, fill = NA) +
  theme_void()

# Using `patchwork` to combine the plots
p | q

Spatial weight matrix

In lattice data analysis, a key concept is the spatial weight matrix, which models the spatial relationships between units in the lattice (i.e., spots or cells). Various methods exist to construct this matrix, such as contiguity-based (direct neighbors), graph-based, or distance-based methods. (Getis 2009; Zuur et al. 2007; Pebesma and Bivand 2023). The documentation of the package spdep gives an overview of the different methods.

For Visium, the most straightforward way is to take the direct neighbours of each spot. This is done using the function findVisiumGraph.

colGraph(sfe_tissue, "visium") <- findVisiumGraph(sfe_tissue)
plotColGraph(sfe_tissue,
  colGraphName = "visium",
  colGeometryName = "spotPoly"
) + theme_void()

In an irregular lattice, the task of finding a spatial weight matrix is more complex, as different options exist. One option is to base the neighbourhood graph on neighbours that are in direct contact with each other (contiguous), as implemented in the poly2nb method.

annotGraph(sfe_tissue, "myofiber_poly2nb") <-
  findSpatialNeighbors(sfe_tissue,
    type = "myofiber_simplified",
    MARGIN = 3,
    method = "poly2nb", # wraps spdep function with same name
    zero.policy = TRUE
  )
p1 <- plotAnnotGraph(sfe_tissue,
  annotGraphName = "myofiber_poly2nb",
  annotGeometryName = "myofiber_simplified"
) + theme_void()

Alternatively, we could take the five nearest neighbours of each myofiber.

annotGraph(sfe_tissue, "knn5") <-
  findSpatialNeighbors(sfe_tissue,
    type = "myofiber_simplified",
    MARGIN = 3, # to use the annotation geometry
    method = "knearneigh", # wraps the spdep function with the same name
    k = 5,
    zero.policy = TRUE
  )
p2 <- plotAnnotGraph(sfe_tissue,
  annotGraphName = "knn5",
  annotGeometryName = "myofiber_simplified"
) + theme_void()

As we can see below, the graphs look quite distinct. On the left side, in the contiguous neighbour graph (neighbours in direct contact) we can notice the formation of patches, while in the kNN graph isolated patches are interconnected.

p1 + p2

Spatial autocorrelation

Spatial autocorrelation measures the association between spatial units while recognizing that the spatial units are not independent measurements. These measures can be global (summarizing the entire field) or local (providing statistics for individual locations) (Pebesma and Bivand 2023).

Global Measures

Global methods give us an overview over the entire field-of-view and summarize the spatial autocorrelation metric to a single value. The metrics are a function of the weight matrix and the variables of interest. The variables of interest can be gene expression, intensity of a marker or the area of the cell. The global measures can be seen as a weighted average of the local metric, as explained below.

In general, a global spatial autocorrelation measure has the form of a double sum over all locations \(i,j\)

\[\sum_i \sum_j f(x_i,x_j) w_{ij}\]

where \(f(x_i,x_j)\) is the measure of association between features of interest and \(w_{ij}\) scales the relationship by a spatial weight as defined in the weight matrix \(W\) (Zuur et al. 2007; Pebesma and Bivand 2023).

Moran’s I

Moran’s I is the most prominent measure of spatial autocorrelation. The values are generally bounded by \(-1\) and \(1\). The expected value is close to \(0\) for large \(n\), the exact value is given by \(\mathbb{E}(I) = -1/(n-1)\). A value higher than the expected value indicates spatial auto-correlation. Negative values indicate negative auto-correlation. Spatial auto-correlation means that similar values tend to be found together in the tissue. In fact, Moran’s I can be interpreted as the Pearson correlation between the value at location \(i\) and the average value of the neighbours of \(i\), (neighbours as defined in the weight matrix \(W\)) (Moran 1950).

In the first example, we will calculate Moran’s I for the number of counts and genes measured in the Visium dataset. First we have a look at the distribution by eye.

plotSpatialFeature(sfe_tissue,
  features = c("nCounts", "nGenes"),
  colGeometryName = "spotPoly"
)

Based on these values we can calculate Moran’s I

calculateUnivariate(t(colData(sfe_tissue)[, c("nCounts", "nGenes")]),
  type = "moran",
  listw = colGraph(sfe_tissue, "visium")
)
DataFrame with 2 rows and 2 columns
            moran         K
        <numeric> <numeric>
nCounts  0.528705   3.00082
nGenes   0.384028   3.88036

We can further use permutation testing to assess statistical significance of our estimates.

sfe_tissue <- colDataUnivariate(sfe_tissue,
  features = c("nCounts", "nGenes"),
  colGraphName = "visium", nsim = 1000,
  type = "moran.mc"
)
plotMoranMC(sfe_tissue, c("nCounts", "nGenes"),
  linewidth = 2
)

Here we can see the Moran’s I values for the number of counts and genes that were identified in the Visium dataset and a permutation-based null distribution. As we can see, the Moran’s I values both indicate positive spatial autocorrelation and are highly significant. This means that regions with similar count and number of measured gene values tend to cluster together. This is an interesting finding and was observed in other spatial transcriptomic datasets, see (Bhuva et al. 2024).

Local Measures for Univariate Data

Often a global measure is not enough. One number determining, e.g., the spatial autocorrelation over an entire tissue slice might not be reflective of tissue heterogeneity. Therefore, local indicators of spatial associations have been developed (Pebesma and Bivand 2023).

Local Moran’s I

Local Moran’s I provides a measure of spatial autocorrelation at each location, highlighting local clusters of similarity or dissimilarity (Anselin 1995). It is defined as:

\[I_i = \frac{x_i - \bar{x}}{\sum_{k=1}^n(x_k-\bar{x})^2/(n-1)} \sum_{j=1}^n w_{ij}(x_j - \bar{x})\]

where the index \(i\) refers to the location for which the measure is calculated. The interpretation is analogous to the global Moran’s I where a value of \(I_i\) higher than \(\mathbb{E}(I) = -w_i/(n-1)\) indicates spatial auto-correlation; smaller values indicate negative auto-correlation. It is important to note that, as for the global counterpart, the value of local Moran’s I could be a result from both the high or low end of the values. Here we will calculate the local Moran’s I value for the measurement of muscle fiber marker gene Myh2:

sfe_tissue <- runUnivariate(sfe_tissue,
  type = "localmoran",
  features = "Myh2",
  colGraphName = "visium",
  swap_rownames = "symbol"
)

# plot the expression values
pExp <- plotSpatialFeature(sfe_tissue,
  features = c("Myh2"), colGeometryName = "spotPoly",
  swap_rownames = "symbol"
)

# plot the local Moran's I values
pLi <- plotLocalResult(sfe_tissue, "localmoran",
  features = c("Myh2"),
  colGeometryName = "spotPoly", swap_rownames = "symbol",
  divergent = TRUE, diverge_center = 0
) + labs(fill = "li(Myh2)") # specify legend


# plot the local Moran's I p-values
pPval <- plotLocalResult(sfe_tissue, "localmoran",
  features = c("Myh2"), "-log10p_adj",
  colGeometryName = "spotPoly", swap_rownames = "symbol",
  divergent = TRUE, diverge_center = -log10(0.05)
) + labs(fill = "-log10(p.adj)") # specify legend

pExp + pLi + pPval

In the local version of Moran’s I, the interpretation is the same as the global version. When interpreting local autocorrelation measures, it is important to consider both the effect size estimates and the significance level. Since the significance level is calculated for each spot separately, it is recommended to adjust for multiple testing. By default, runUnivariate uses the Benjamini & Hochberg correction. The local Moran’s I statistics reveal locations in the tissue that have similar values to their neighbours (c.f., the lower part of the tissue).

TipEXERCISE

Try to recalculate the local Moran’s I with gene Myh1. What do you notice?

sfe_tissue <- runUnivariate(sfe_tissue,
  type = "localmoran",
  features = "Myh1",
  colGraphName = "visium",
  swap_rownames = "symbol"
)

# plot the expression values
pExp <- plotSpatialFeature(sfe_tissue,
  features = c("Myh1"), colGeometryName = "spotPoly",
  swap_rownames = "symbol"
)

# plot the local Moran's I values
pLi <- plotLocalResult(sfe_tissue, "localmoran",
  features = c("Myh1"),
  colGeometryName = "spotPoly", swap_rownames = "symbol",
  divergent = TRUE, diverge_center = 0
) + labs(fill = "li(Myh1)") # specify legend


# plot the local Moran's I p-values
pPval <- plotLocalResult(sfe_tissue, "localmoran",
  features = c("Myh1"), "-log10p_adj",
  colGeometryName = "spotPoly", swap_rownames = "symbol",
  divergent = TRUE, diverge_center = -log10(0.05)
) + labs(fill = "-log10(p.adj)") # specify legend

pExp + pLi + pPval

Advanced – Bivariate Moran’s I

There exists a bivariate version of Moran’s I as well. The package spatialDM (Li et al. 2023) uses an adapted version of bivariate Moran’s I to identify ligand-receptor pairs.

For two continuous observations, the global bivariate Moran’s I is defined as (Wartenberg 1985; Bivand 2022) \[I_B =\frac{\sum_i x_i \sum_j w_{ij} y_j}{\sum_i x_i^2}.\]

Here, \(x_i\) is the value of variable \(x\) at location \(i\), while \(\sum_j w_{ij}y_j\) is the spatial lag of variable \(y\), i.e. a weighted summary of \(y\) in the neighbours of location \(i\) (also called the spatial lag of \(y\)).

There exists a local version that we will explore here. The local version follows that same idea; it compares \(x_i\) with the spatial lag of a second variable \(y\). Note that bivariate Moran’s I is not symmetric.

sfe_tissue <- runBivariate(sfe_tissue, 
                           "localmoran_bv", # wraps the method from spdep
                           c("Myh1", "Myh2"),
                           swap_rownames = "symbol", nsim = 1000)

# this command gives all results of localmoran_bv
localResultFeatures(sfe_tissue, "localmoran_bv")
[1] "Myh1__Myh1" "Myh2__Myh1" "Myh1__Myh2" "Myh2__Myh2"
plotLocalResult(sfe_tissue, "localmoran_bv", 
                c("Myh1__Myh2", "Myh2__Myh1"),
                colGeometryName = "spotPoly",
                divergent = TRUE, diverge_center = 0)

TipEXERCISE

Compare the results with the univariate local Moran’s I plots for genes Myh1 and Myh2. Also compare Myh1–Myh2 with Myh2–Myh1. What do you notice?

The univariate statistic compares each gene with neighbouring values of itself, whereas the bivariate statistic compares one gene with neighbouring values of the other gene.

The Myh1–Myh2 and Myh2–Myh1 maps are similar, but not identical. This is expected because bivariate Moran’s I is asymmetric.

Please note that the result might overestimate the spatial autocorrelation of the variables due to the inherent (non-spatial) correlation of \(x\) and \(y\) (Bivand 2022).

Impact of neighbourhood on autocorrelation measures

Let’s compare the impact of different spatial weight matrices on local autocorrelation analysis. In this Visium dataset we do not have gene expression information for each cell because of the resolution of the spots. We will instead calculate the autocorrelation measures on the area of the segmented myofibers. This allows us to ask whether similarly sized myofibers tend to cluster spatially.

First, we base the spatial weight matrix on contiguous neighbours.

sfe_tissue <- annotGeometryUnivariate(sfe_tissue, 
  "localmoran", "area",
  annotGeometryName = "myofiber_simplified",
  annotGraphName = "myofiber_poly2nb",
  include_self = FALSE, zero.policy = TRUE,
  name = "myofiber_poly2nb"
)

pPoly <- plotLocalResult(sfe_tissue, "myofiber_poly2nb", "area",
  annotGeometryName = "myofiber_simplified",
  annotGraphName = "myofiber_poly2nb",
  divergent = TRUE, diverge_center = 0
) + labs(title = "Local Moran's I (li)\nContiguous Neighbours", 
         fill = "li(area)")

As an alternative, we also base the spatial weight matrix on the 5 nearest neighbours.

sfe_tissue <- annotGeometryUnivariate(sfe_tissue, "localmoran", "area",
  annotGeometryName = "myofiber_simplified",
  annotGraphName = "knn5",
  include_self = FALSE, zero.policy = TRUE,
  name = "knn5"
)

pKnn5 <- plotLocalResult(sfe_tissue, "knn5", "area",
  annotGeometryName = "myofiber_simplified",
  annotGraphName = "knn5",
  divergent = TRUE, diverge_center = 0
) + labs(title = "Local Moran's I (li)\n5 Nearest Neighbours", 
         fill = "li(area)")
pPoly + pKnn5

TipEXERCISE

What happens if we set the nearest neighbourhood graph to 10 nearest neighbours? What do you notice?

annotGraph(sfe_tissue, "knn10") <-
  findSpatialNeighbors(
    sfe_tissue,
    type = "myofiber_simplified",
    MARGIN = 3,
    method = "knearneigh",
    k = 10,
    zero.policy = TRUE
  )

sfe_tissue <- annotGeometryUnivariate(
  sfe_tissue,
  "localmoran",
  "area",
  annotGeometryName = "myofiber_simplified",
  annotGraphName = "knn10",
  include_self = FALSE,
  zero.policy = TRUE,
  name = "knn10"
)

pKnn10 <- plotLocalResult(
  sfe_tissue,
  "knn10",
  "area",
  annotGeometryName = "myofiber_simplified",
  annotGraphName = "knn10",
  divergent = TRUE,
  diverge_center = 0
) +
  labs(
    title = "Local Moran's I (li)\n10 Nearest Neighbours",
    fill = "li(area)"
  )

pKnn5 + pKnn10

The overall spatial pattern remains similar, but there are local differences. With 10 nearest neighbours, each myofiber is compared to a broader neighbourhood, resulting in a slightly smoother local Moran’s I pattern. Note that the same color in the different plots does not correspond to the same value between the plots.

Summary and Considerations

  • Lattice data refers to spatial data collected at fixed locations arranged in regular or irregular grids, contrasting with stochastic point pattern analysis.
  • Regular lattices have uniform, evenly spaced spots, while irregular lattices have variable sizes and shapes with non-uniform spacing.
  • The spatial weight matrix models spatial relationships between lattice units. Different methods exist to define the spatial weight matrix.
  • Global spatial autocorrelation measures, like Moran’s I, summarize spatial correlation over the entire field, while local measures identify local clusters of similarity or dissimilarity.
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] scater_1.38.1                   ggplot2_4.0.3                  
 [3] scuttle_1.20.0                  SingleCellExperiment_1.32.0    
 [5] SummarizedExperiment_1.40.0     Biobase_2.70.0                 
 [7] GenomicRanges_1.62.1            Seqinfo_1.0.0                  
 [9] IRanges_2.44.0                  S4Vectors_0.48.0               
[11] BiocGenerics_0.56.0             generics_0.1.4                 
[13] MatrixGenerics_1.22.0           matrixStats_1.5.0              
[15] sf_1.1-0                        SFEData_1.12.0                 
[17] Voyager_1.12.0                  SpatialFeatureExperiment_1.12.1

loaded via a namespace (and not attached):
  [1] splines_4.5.3             bitops_1.0-9             
  [3] filelock_1.0.3            tibble_3.3.1             
  [5] R.oo_1.27.1               lifecycle_1.0.5          
  [7] httr2_1.2.2               edgeR_4.8.2              
  [9] lattice_0.22-9            MASS_7.3-65              
 [11] backports_1.5.1           magrittr_2.0.5           
 [13] limma_3.66.0              rmarkdown_2.32           
 [15] yaml_2.3.12               sp_2.2-1                 
 [17] DBI_1.3.0                 RColorBrewer_1.1-3       
 [19] multcomp_1.4-32           abind_1.4-8              
 [21] spatialreg_1.4-3          purrr_1.2.2              
 [23] R.utils_2.13.0            RCurl_1.98-1.19          
 [25] TH.data_1.1-5             rappdirs_0.3.4           
 [27] sandwich_3.1-3            ggrepel_0.9.8            
 [29] irlba_2.3.7               terra_1.9-25             
 [31] BiocStyle_2.38.0          units_1.0-1              
 [33] RSpectra_0.16-2           dqrng_0.3.2              
 [35] DelayedMatrixStats_1.32.0 codetools_0.2-20         
 [37] DropletUtils_1.30.0       DelayedArray_0.36.1      
 [39] tidyselect_1.2.1          memuse_4.2-3             
 [41] farver_2.1.2              ScaledMatrix_1.18.0      
 [43] viridis_0.6.5             BiocFileCache_3.0.0      
 [45] jsonlite_2.0.0            BiocNeighbors_2.4.0      
 [47] e1071_1.7-17              survival_3.8-6           
 [49] tools_4.5.3               ggnewscale_0.5.2         
 [51] Rcpp_1.1.1-1.1            glue_1.8.1               
 [53] gridExtra_2.3.1           SparseArray_1.10.8       
 [55] xfun_0.57                 EBImage_4.52.0           
 [57] dplyr_1.2.1               HDF5Array_1.38.0         
 [59] withr_3.0.3               BiocManager_1.30.27      
 [61] fastmap_1.2.0             boot_1.3-32              
 [63] rhdf5filters_1.22.0       bluster_1.20.0           
 [65] spData_2.3.5              digest_0.6.39            
 [67] rsvd_1.0.5                R6_2.6.1                 
 [69] wk_0.9.5                  LearnBayes_2.15.2        
 [71] jpeg_0.1-11               RSQLite_3.53.1           
 [73] R.methodsS3_1.8.2         h5mread_1.2.1            
 [75] data.table_1.17.8         class_7.3-23             
 [77] httr_1.4.9                htmlwidgets_1.6.4        
 [79] S4Arrays_1.10.1           spdep_1.4-2              
 [81] pkgconfig_2.0.3           scico_1.5.0              
 [83] gtable_0.3.6              blob_1.3.0               
 [85] S7_0.2.2                  XVector_0.50.0           
 [87] htmltools_0.5.9           fftwtools_0.9-11         
 [89] scales_1.4.0              png_0.1-9                
 [91] SpatialExperiment_1.20.0  knitr_1.51               
 [93] rjson_0.2.23              coda_0.19-4.1            
 [95] nlme_3.1-169              curl_7.1.0               
 [97] proxy_0.4-29              cachem_1.1.0             
 [99] zoo_1.8-15                rhdf5_2.54.1             
[101] BiocVersion_3.22.0        KernSmooth_2.23-26       
[103] parallel_4.5.3            vipor_0.4.7              
[105] AnnotationDbi_1.72.0      s2_1.1.9                 
[107] pillar_1.11.1             grid_4.5.3               
[109] vctrs_0.7.3               BiocSingular_1.26.1      
[111] dbplyr_2.6.0              beachmat_2.26.0          
[113] sfheaders_0.4.5           cluster_2.1.8.2          
[115] beeswarm_0.4.0            evaluate_1.0.5           
[117] zeallot_0.2.0             magick_2.9.1             
[119] mvtnorm_1.4-1             cli_3.6.6                
[121] locfit_1.5-9.12           compiler_4.5.3           
[123] rlang_1.2.0               crayon_1.5.3             
[125] labeling_0.4.3            classInt_0.4-11          
[127] ggbeeswarm_0.7.3          viridisLite_0.4.3        
[129] deldir_2.0-4              BiocParallel_1.44.0      
[131] Biostrings_2.78.0         tiff_0.1-12              
[133] marginaleffects_1.0.0     Matrix_1.7-5             
[135] ExperimentHub_3.0.0       patchwork_1.3.2          
[137] sparseMatrixStats_1.22.0  bit64_4.8.2              
[139] Rhdf5lib_1.32.0           KEGGREST_1.50.0          
[141] statmod_1.5.2             AnnotationHub_4.0.0      
[143] igraph_2.3.2              memoise_2.0.1            
[145] bit_4.6.0

This tutorial is partially based on the following Voyager vignette.

References

Anselin, Luc. 1995. “Local Indicators of Spatial Association—LISA.” Geographical Analysis 27 (2): 93–115. https://doi.org/10.1111/j.1538-4632.1995.tb00338.x.
Bhuva, Dharmesh D., Chin Wee Tan, Agus Salim, et al. 2024. “Library Size Confounds Biology in Spatial Transcriptomics Data.” Genome Biology 25 (1): 99. https://doi.org/10.1186/s13059-024-03241-7.
Bivand, Roger. 2022. “R Packages for Analyzing Spatial Data: A Comparative Case Study with Areal Data.” Geographical Analysis 54 (3): 488–518. https://doi.org/10.1111/gean.12319.
Getis, Arthur. 2009. “Spatial Weights Matrices.” Geographical Analysis 41 (4): 404–10. https://doi.org/10.1111/j.1538-4632.2009.00768.x.
Li, Zhuoxuan, Tianjie Wang, Pentao Liu, and Yuanhua Huang. 2023. “SpatialDM for Rapid Identification of Spatially Co-Expressed Ligand–Receptor and Revealing Cell–Cell Communication Patterns.” Nature Communications 14 (1): 3995. https://doi.org/10.1038/s41467-023-39608-w.
McKellar, David W., Lauren D. Walter, Leo T. Song, et al. 2021. “Large-Scale Integration of Single-Cell Transcriptomic Data Captures Transitional Progenitor States in Mouse Skeletal Muscle Regeneration.” Communications Biology 4 (1): 1–12. https://doi.org/10.1038/s42003-021-02810-x.
Moran, P. A. P. 1950. “Notes on Continuous Stochastic Phenomena.” Biometrika 37 (1/2): 17–23. https://doi.org/10.2307/2332142.
Moses, Lambda, Pétur Helgi Einarsson, Kayla C. Jackson, et al. 2023. Voyager: Exploratory Single-Cell Genomics Data Analysis with Geospatial Statistics. bioRxiv. https://doi.org/10.1101/2023.07.20.549945.
Pebesma, Edzer, and Roger Bivand. 2023. Spatial Data Science: With Applications in R. 1st ed. Chapman and Hall/CRC. https://doi.org/10.1201/9780429459016.
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.
Ståhl, Patrik L., Fredrik Salmén, Sanja Vickovic, et al. 2016. “Visualization and Analysis of Gene Expression in Tissue Sections by Spatial Transcriptomics.” Science 353 (6294): 78–82. https://doi.org/10.1126/science.aaf2403.
Wartenberg, Daniel. 1985. “Multivariate Spatial Correlation: A Method for Exploratory Geographical Analysis.” Geographical Analysis 17 (4): 263–83. https://doi.org/10.1111/j.1538-4632.1985.tb00849.x.
Zuur, Alain F., Elena N. Ieno, and Graham M. Smith. 2007. Analysing Ecological Data. Statistics for Biology and Health. Springer.