1 Introduction

For each cell, lisaClust (Patrick et al. 2023) counts the cells of every type within a few radii of it and compares each count with what would be expected if that type were spread evenly over the image. These local indicators of spatial association (LISA) describe the cell’s surroundings: a tumour cell deep in a tumour nest has many tumour neighbours and few immune neighbours, while a tumour cell at the invasive margin has both. Clustering the cells by their LISA groups cells with similar surroundings, whatever their own type, into regions.

The number of regions sets the scale of the question. With a few regions, they are broad tissue domains, such as tumour, stroma and immune infiltrate. With 20 or more, they are cellular niches: recurring local combinations of cell types, such as T cells mixed with macrophages at a tumour border.

lisaClust needs only the type and position of every cell, so it applies to any technology that gives segmented, typed cells: spatial proteomics such as imaging mass cytometry, MIBI and CODEX, and spatial transcriptomics such as Xenium, CosMx, MERSCOPE or high-definition spatial transcriptomics (Patrick et al. 2023). It is not designed for spot-based data such as Visium, where a spot holds several cells. It works alongside other SydneyBioX packages: spicyR tests whether pairs of cell types co-localise differently between groups of patients, and Statial tests co-localisation within tissue compartments. Here we use spicyR’s colTest() to compare regions between groups of patients, and tidySpatialExperiment to summarise the cells of each image with dplyr verbs.

2 Installation

if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("lisaClust", "spicyR", "tidySpatialExperiment", "SpatialDatasets"))
library(lisaClust)
library(spicyR)
library(tidySpatialExperiment)
library(SpatialDatasets)
library(ggplot2)
theme_set(theme_classic())

3 Tissue domains in triple-negative breast cancer

3.1 The data

Keren et al. (2018) imaged tumours from 40 patients with triple-negative breast cancer by MIBI-TOF, one image per patient, and classified each tumour by how its immune and tumour cells are arranged: compartmentalised, where immune cells and tumour cells occupy separate areas; mixed, where they are interspersed; and cold, with few immune cells. The data are a SpatialExperiment in the SpatialDatasets package, with the coordinates in pixels.

kerenSPE <- spe_Keren_2018()
kerenSPE |> distinct(imageID, tumour_type) |> count(tumour_type)
#> # A tibble: 3 × 2
#>   tumour_type           n
#>   <fct>             <int>
#> 1 cold                  6
#> 2 compartmentalised    15
#> 3 mixed                19
sort(table(kerenSPE$cellType), decreasing = TRUE)
#> 
#> Keratin_Tumour    Macrophages     CD8_T_cell     CD4_T_cell         B_cell 
#>          99487          20616          15698          12438           9115 
#>    Mesenchymal   Other_Immune     DC_or_Mono       dn_T_CD3         Tumour 
#>           8170           6891           5049           3848           3167 
#>    Mono_or_Neu    Neutrophils    Endothelial   Unidentified          Tregs 
#>           3110           3018           2086           1725           1341 
#>             DC             NK 
#>           1245            674

3.2 Finding domains

lisaClust() computes the LISA of every cell and clusters them with k-means into k regions, stored in a new colData column (region by default). It looks for columns called imageID and cellType, and takes the coordinates of a SpatialExperiment from spatialCoords(). Here we ask for four domains, describing each cell’s surroundings within 20, 50 and 100 pixels.

set.seed(51773)
kerenSPE <- lisaClust(kerenSPE, k = 4, r = c(20, 50, 100))
table(kerenSPE$region)
#> 
#> region_1 region_2 region_3 region_4 
#>    62148   104863    23516     7151

k-means starts from random centres, so set a seed to make the regions reproducible.

3.3 What is in each domain?

regionMap() shows how much more often each cell type is found in each region than if cell types were spread evenly over the regions: above one, the type is enriched in that region.

regionMap(kerenSPE, type = "bubble")

The regions are numbered arbitrarily. nameRegions() names each one by the marker cell type it is most enriched for, here among tumour cells, CD8 T cells, B cells and mesenchymal cells. Regions most enriched for the same marker get the same name, so they are merged into one domain.

markers <- c(tumour = "Keratin_Tumour", `T cell` = "CD8_T_cell", `B cell` = "B_cell", stroma = "Mesenchymal")
kerenSPE <- nameRegions(kerenSPE, markers, regionName = "domain")
round(100 * prop.table(table(kerenSPE$domain)), 1)
#> 
#> B cell T cell stroma tumour 
#>    3.6   11.9   31.4   53.0

The cell types in each domain:

regionMap(kerenSPE, region = "domain", type = "bubble")

3.4 Domains in two tumours

hatchingPlot() draws the cells coloured by type, with each region outlined and filled with its own hatching, so types and regions can be seen together. Here are the compartmentalised tumour and the mixed tumour with the most cells in the T cell domain. With tidySpatialExperiment, group_by() and summarise() work on the cells of a SpatialExperiment and return a table.

tcellShare <- kerenSPE |>
  group_by(imageID, tumour_type) |>
  summarise(share = mean(domain == "T cell"), .groups = "drop")
examples <- tcellShare |>
  filter(tumour_type != "cold") |>
  group_by(tumour_type) |>
  slice_max(share)
examples
#> # A tibble: 2 × 3
#> # Groups:   tumour_type [2]
#>   imageID tumour_type       share
#>   <chr>   <fct>             <dbl>
#> 1 16      compartmentalised 0.403
#> 2 29      mixed             0.170
hatchingPlot(kerenSPE, useImages = examples$imageID, region = "domain")

3.5 Do the domains differ between kinds of tumour?

Each patient has one image, so the share of each domain in an image is a value per patient. regionBoxPlot() draws these shares by group.

regionBoxPlot(kerenSPE, condition = "tumour_type", region = "domain")

spicyR’s colTest() computes the share of each domain in each image and compares the groups with a Wilcoxon test; here we compare compartmentalised with mixed tumours.

domainTests <- kerenSPE |>
  filter(tumour_type != "cold") |>
  colTest(condition = "tumour_type", feature = "domain", type = "wilcox")
domainTests
#>        tval.W    pval adjPval cluster
#> tumour      6 3.2e-08 1.3e-07  tumour
#> T cell    250 9.1e-05 1.8e-04  T cell
#> B cell    240 7.5e-04 1.0e-03  B cell
#> stroma    230 3.1e-03 3.1e-03  stroma
medians <- with(tcellShare, tapply(share, tumour_type, median))
round(100 * medians, 1)
#>              cold compartmentalised             mixed 
#>               1.3              15.2               4.5
pT <- domainTests["T cell", "pval"]

The T cell domain makes up a median of 15.2% of the cells in a compartmentalised tumour and 4.5% in a mixed one (Wilcoxon p = 9.1e-05). In compartmentalised tumours the immune cells form their own areas, which become a domain; in mixed tumours they sit among the tumour cells and fall into the tumour domain.

Keren et al. (2018) classified the tumours from how their immune and tumour cells are arranged, so this difference is a check that the domains capture that arrangement, not an independent finding.

4 Cellular niches in colorectal cancer

4.1 The data

Schürch et al. (2020) imaged colorectal cancers from 35 patients by CODEX, four tissue microarray cores per patient. The patients had one of two patterns of immune response: a Crohn’s-like reaction (CLR, with lymphoid follicles; 17 patients) or diffuse inflammatory infiltration (DII; 18 patients). We leave out the objects the authors labelled as debris (dirt).

schurchSPE <- spe_Schurch_2020() |>
  filter(cellType != "dirt") |>
  mutate(condition = ifelse(group == 1, "CLR", "DII"))
schurchSPE |> distinct(patients, condition) |> count(condition)
#> # A tibble: 2 × 2
#>   condition     n
#>   <chr>     <int>
#> 1 CLR          17
#> 2 DII          18
length(unique(schurchSPE$cellType))
#> [1] 28

4.2 Finding niches

With 20 regions, the regions are finer combinations of cell types. We store them in a column called niche.

set.seed(51773)
schurchSPE <- lisaClust(schurchSPE, k = 20, r = c(20, 50, 100), regionName = "niche")
regionMap(schurchSPE, region = "niche", type = "bubble") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

4.3 Comparing niches between groups of patients

The question is whether a niche makes up more of the tissue in one group of patients than in the other. The four cores of a patient come from the same tumour, so they are not independent: we pool the cells of each patient’s cores, so that each patient contributes one share of each niche, and compare the two groups of patients. Setting imageID = "patients" in colTest() does this pooling. With 20 niches tested, colTest() also adjusts the p-values for multiple testing (adjPval).

patientTests <- colTest(schurchSPE, condition = "condition", feature = "niche",
                        imageID = "patients", type = "wilcox")
head(patientTests, 5)
#>           tval.W   pval adjPval   cluster
#> region_13    250 0.0016   0.019 region_13
#> region_7     240 0.0019   0.019  region_7
#> region_11     68 0.0043   0.029 region_11
#> region_4      74 0.0083   0.042  region_4
#> region_5     220 0.0200   0.080  region_5

For comparison, the same test treating each core as an independent sample:

coreTests <- colTest(schurchSPE, condition = "condition", feature = "niche",
                     imageID = "imageID", type = "wilcox")
c(patients = sum(patientTests$adjPval < 0.05), cores = sum(coreTests$adjPval < 0.05))
#> patients    cores 
#>        4        4

With the 35 patients as the units, 4 niches differ between CLR and DII patients at a 5% false discovery rate (smallest adjusted p-value 0.019). Treating each of the 140 cores as an independent sample instead, 4 niches differ, with a smallest adjusted p-value of 0.00088.

Cores or images from the same patient are alike, so treating them as independent overstates the evidence and can turn up differences that the patients do not support. Compare niches with patients as the units, and adjust for the number of niches tested.

regionMap(type = "table") gives the enrichment of each cell type in each niche as a matrix, so we can list the cell types the niche with the smallest p-value is most enriched for, and its median share in each group.

top <- patientTests$cluster[1]
enrichment <- regionMap(schurchSPE, region = "niche", type = "table")
topTypes <- sort(enrichment[, top], decreasing = TRUE)
round(topTypes[1:5], 1)
#>            b_cells        cd4_t_cells        cd3_t_cells         cd11c_d_cs 
#>                8.3                2.1                2.1                1.4 
#> cd4_t_cells_cd45ro 
#>                1.4

topMedians <- schurchSPE |>
  group_by(patients, condition) |>
  summarise(share = mean(niche == top), .groups = "drop") |>
  with(tapply(share, condition, median))
round(100 * topMedians, 1)
#> CLR DII 
#> 6.4 3.4

The niche with the smallest p-value, region_13, is most enriched for b_cells, cd4_t_cells, cd3_t_cells. It makes up a median of 6.4% of the cells in CLR patients and 3.4% in DII patients, as expected for the lymphoid follicles that define a Crohn’s-like reaction:

regionBoxPlot(schurchSPE, condition = "condition", region = "niche", imageID = "patients", regions = top)

5 Using your own clustering

lisa() returns the LISA themselves, one row per cell and one column per radius and cell type, so you can cluster them any way you like.

curves <- lisa(kerenSPE[, kerenSPE$imageID %in% examples$imageID], r = c(20, 50, 100))
dim(curves)
#> [1] 13031    51
curves[1:3, 1:4]
#>        20_Keratin_Tumour 20_CD8_T_cell 20_dn_T_CD3 20_CD4_T_cell
#> cell_1         3.7228902    -0.8831332  -0.2041394    -0.8778542
#> cell_2        -0.7323238     5.5545686  -0.2041394    -0.8778542
#> cell_3        -0.7323238     5.4766730  -0.2041394    -0.8778542

# for example, k-means with several random starts
set.seed(51773)
km <- kmeans(curves, centers = 4, nstart = 10, iter.max = 50)
table(km$cluster)
#> 
#>    1    2    3    4 
#> 5466 1362 6043  160

6 Choosing the settings

The number of regions. A few regions describe the large compartments of a tissue; many regions separate finer neighbourhoods but split rarer combinations across several regions and make each one harder to interpret. Look at regionMap() for a few values of k and choose the coarsest one that separates the structures you care about; generic rules for the number of clusters, such as the elbow or silhouette, can merge tissue structures that are known to be distinct (Patrick et al. 2023).

The radii. Choose them from the scales of interaction you expect, in the units of the coordinates: a radius of two or three cell diameters describes the cell’s immediate neighbours, and larger radii describe its wider surroundings. Several radii (r) are used together; run time grows with the number and size of the radii, so a few well-chosen ones are better than many.

Uneven cell density. By default the expected count assumes cells are spread evenly over the image. With sigma, each neighbour is weighted by the inverse of the local density of all cells, estimated with a Gaussian kernel of bandwidth sigma, so that a dense area does not by itself look enriched for every cell type.

7 How it works

For cell i, cell type j and radius r, let nij(r) be the number of cells of type j within r of cell i, counting cell i itself when it is of type j. If the cells of type j were spread evenly over the image window W, with density \(\lambda_j = N_j / |W|\), the expected count would be

\[E_{ij}(r) = \lambda_j \, \pi r^2 e_i(r),\]

where \(e_i(r)\) is the share of the disc of radius r around cell i that lies inside the window, which corrects for cells near the edge. The LISA is the standardised difference

\[\frac{n_{ij}(r) - E_{ij}(r)}{\sqrt{E_{ij}(r)}},\]

a local version of Ripley’s K function for the pair (lisaFunc = "L" gives \(\sqrt{n_{ij}(r)} - \sqrt{E_{ij}(r)}\) instead). The window is the convex hull of the image’s cells by default (window). The neighbour search and the edge correction run in C++, so the LISA of hundreds of thousands of cells take seconds.

8 Reporting results

A methods sentence might read: “We used lisaClust (version 1.21.7) to compute, for each cell, local indicators of spatial association with every cell type at radii of 20, 50 and 100 pixels, and clustered the cells into 20 niches by k-means. The share of each niche in each patient was compared between groups with a Wilcoxon test (spicyR’s colTest), with patients as the units, and p-values were adjusted across niches by the Benjamini–Hochberg method.” Show regionMap() and a hatchingPlot() of example images with the results.

citation("lisaClust")
#> To cite package 'lisaClust' in publications use:
#> 
#>   Patrick E, Canete NP, Iyengar SS, Harman AN, Sutherland GT, Yang P
#>   (2023). "Spatial analysis for highly multiplexed imaging data to
#>   identify tissue microenvironments." _Cytometry Part A_, *103*,
#>   593-599. doi:10.1002/cyto.a.24729
#>   <https://doi.org/10.1002/cyto.a.24729>.
#> 
#> A BibTeX entry for LaTeX users is
#> 
#>   @Article{,
#>     title = {Spatial analysis for highly multiplexed imaging data to identify tissue microenvironments},
#>     author = {Ellis Patrick and Nicolas P. Canete and Sourish S. Iyengar and Andrew N. Harman and Greg T. Sutherland and Pengyi Yang},
#>     journal = {Cytometry Part A},
#>     volume = {103},
#>     pages = {593--599},
#>     year = {2023},
#>     doi = {10.1002/cyto.a.24729},
#>   }

9 Session information

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] SpatialDatasets_1.11.0          ExperimentHub_3.3.2            
#>  [3] AnnotationHub_4.3.2             BiocFileCache_3.3.0            
#>  [5] dbplyr_2.6.0                    tidySpatialExperiment_1.9.0    
#>  [7] ggplot2_4.0.3                   tidyr_1.3.2                    
#>  [9] dplyr_1.2.1                     tidySingleCellExperiment_1.23.0
#> [11] ttservice_0.5.3                 SpatialExperiment_1.23.0       
#> [13] SingleCellExperiment_1.35.2     SummarizedExperiment_1.43.0    
#> [15] Biobase_2.73.2                  GenomicRanges_1.65.4           
#> [17] Seqinfo_1.3.2                   IRanges_2.47.5                 
#> [19] S4Vectors_0.51.10               BiocGenerics_0.59.12           
#> [21] generics_0.1.4                  MatrixGenerics_1.25.0          
#> [23] matrixStats_1.5.0               spicyR_1.99.11                 
#> [25] lisaClust_1.21.7                BiocStyle_2.41.0               
#> 
#> loaded via a namespace (and not attached):
#>   [1] RColorBrewer_1.1-3     jsonlite_2.0.0         magrittr_2.0.5        
#>   [4] spatstat.utils_3.2-5   magick_2.9.1           farver_2.1.2          
#>   [7] rmarkdown_2.32         vctrs_0.7.3            memoise_2.0.1         
#>  [10] spatstat.explore_3.8-3 tinytex_0.61           htmltools_0.5.9       
#>  [13] S4Arrays_1.13.2        BiocBaseUtils_1.15.1   curl_8.0.0            
#>  [16] SparseArray_1.13.4     sass_0.4.10            bslib_0.12.0          
#>  [19] htmlwidgets_1.6.4      httr2_1.3.0            plotly_4.12.1         
#>  [22] cachem_1.1.0           mime_0.13              lifecycle_1.0.5       
#>  [25] pkgconfig_2.0.3        Matrix_1.7-6           R6_2.6.1              
#>  [28] fastmap_1.2.0          shiny_1.14.0           digest_0.6.39         
#>  [31] AnnotationDbi_1.75.2   tensor_1.5.1           RSQLite_3.53.3        
#>  [34] labeling_0.4.3         filelock_1.0.3         fansi_1.0.7           
#>  [37] spatstat.sparse_3.2-0  httr_1.4.9             polyclip_1.10-7       
#>  [40] abind_1.4-8            compiler_4.6.1         bit64_4.8.6           
#>  [43] withr_3.0.3            S7_0.2.2               BiocParallel_1.47.0   
#>  [46] viridis_0.6.5          DBI_1.3.0              rappdirs_0.3.4        
#>  [49] DelayedArray_0.39.8    rjson_0.2.23           tools_4.6.1           
#>  [52] otel_0.2.0             httpuv_1.6.17          goftest_1.2-3         
#>  [55] glue_1.8.1             nlme_3.1-171           promises_1.5.0        
#>  [58] grid_4.6.1             gtable_0.3.6           spatstat.data_3.1-9   
#>  [61] data.table_1.18.6.1    utf8_1.2.6             XVector_0.53.0        
#>  [64] spatstat.geom_3.8-3    BiocVersion_3.24.0     pillar_1.11.1         
#>  [67] stringr_1.6.0          later_1.4.8            splines_4.6.1         
#>  [70] lattice_0.23-1         survival_3.8-12        bit_4.6.0             
#>  [73] deldir_2.0-4           tidyselect_1.2.1       Biostrings_2.81.9     
#>  [76] knitr_1.52             gridExtra_2.3.1        bookdown_0.48         
#>  [79] xfun_0.61              pheatmap_1.0.13        stringi_1.8.9         
#>  [82] fftwtools_0.9-11       yaml_2.3.12            evaluate_1.0.5        
#>  [85] codetools_0.2-20       tibble_3.3.1           BiocManager_1.30.27   
#>  [88] cli_3.6.6              xtable_1.8-8           jquerylib_0.1.4       
#>  [91] dichromat_2.0-1        Rcpp_1.1.2             spatstat.random_3.5-2 
#>  [94] png_0.1-9              spatstat.univar_3.2-0  parallel_4.6.1        
#>  [97] tidygate_1.0.19        ellipsis_0.3.3         blob_1.3.0            
#> [100] viridisLite_0.4.3      scales_1.4.0           purrr_1.2.2           
#> [103] crayon_1.5.3           rlang_1.3.0            KEGGREST_1.53.6

References

Keren, L, M Bosse, D Marquez, R Angoshtari, S Jain, S Varma, S Yang, et al. 2018. “A Structured Tumor-Immune Microenvironment in Triple Negative Breast Cancer Revealed by Multiplexed Ion Beam Imaging.” Cell 174 (6): 1373–1387.e19.
Patrick, Ellis, Nicolas P Canete, Sourish S Iyengar, Andrew N Harman, Greg T Sutherland, and Pengyi Yang. 2023. “Spatial Analysis for Highly Multiplexed Imaging Data to Identify Tissue Microenvironments.” Cytometry Part A 103: 593–99. https://doi.org/10.1002/cyto.a.24729.
Schürch, Christian M et al. 2020. “Coordinated Cellular Neighborhoods Orchestrate Antitumoral Immunity at the Colorectal Cancer Invasive Front.” Cell 182 (5): 1341–59. https://doi.org/10.1016/j.cell.2020.07.005.