Introduction to lisaClust
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.
Installation
Tissue domains in triple-negative breast cancer
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 674Finding 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 7151k-means starts from random centres, so set a seed to make the regions reproducible.
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.
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 stroma T cell tumour
#> 3.6 31.4 11.9 53.0The cell types in each domain:
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
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.
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.
Cellular niches in colorectal cancer
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] 28Finding 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))
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_5For 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 4With 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.4The 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)
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 160Choosing 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.
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.
Reporting results
A methods sentence might read: “We used lisaClust (version 1.21.8) 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},
#> }Session information
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [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: Etc/UTC
#> 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.8 BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] RColorBrewer_1.1-3 sys_3.4.3 jsonlite_2.0.0
#> [4] magrittr_2.0.5 spatstat.utils_3.2-5 magick_2.9.1
#> [7] farver_2.1.2 rmarkdown_2.32 vctrs_0.7.3
#> [10] memoise_2.0.1 spatstat.explore_3.8-3 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 buildtools_1.0.0 mime_0.13
#> [25] lifecycle_1.0.5 pkgconfig_2.0.3 Matrix_1.7-6
#> [28] R6_2.6.1 fastmap_1.2.0 shiny_1.14.0
#> [31] digest_0.6.39 AnnotationDbi_1.75.2 tensor_1.5.1
#> [34] RSQLite_3.53.3 filelock_1.0.3 labeling_0.4.3
#> [37] fansi_1.0.7 spatstat.sparse_3.2-0 httr_1.4.9
#> [40] polyclip_1.10-7 abind_1.4-8 compiler_4.6.1
#> [43] bit64_4.8.6 withr_3.0.3 S7_0.2.2
#> [46] BiocParallel_1.47.0 viridis_0.6.5 DBI_1.3.0
#> [49] rappdirs_0.3.4 DelayedArray_0.39.8 rjson_0.2.23
#> [52] tools_4.6.1 otel_0.2.0 httpuv_1.6.17
#> [55] goftest_1.2-3 glue_1.8.1 nlme_3.1-171
#> [58] promises_1.5.0 grid_4.6.1 gtable_0.3.6
#> [61] spatstat.data_3.1-9 data.table_1.18.6.1 utf8_1.2.6
#> [64] XVector_0.53.0 spatstat.geom_3.8-3 BiocVersion_3.24.0
#> [67] pillar_1.11.1 stringr_1.6.0 later_1.4.8
#> [70] splines_4.6.1 lattice_0.23-1 survival_3.8-12
#> [73] bit_4.6.0 deldir_2.0-4 tidyselect_1.2.1
#> [76] maketools_1.3.2 Biostrings_2.81.9 knitr_1.52
#> [79] gridExtra_2.3.1 xfun_0.61 pheatmap_1.0.13
#> [82] stringi_1.8.9 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] Rcpp_1.1.2 spatstat.random_3.5-2 png_0.1-9
#> [94] spatstat.univar_3.2-0 parallel_4.6.1 tidygate_1.0.19
#> [97] ellipsis_0.3.3 blob_1.3.0 viridisLite_0.4.3
#> [100] scales_1.4.0 purrr_1.2.2 crayon_1.5.3
#> [103] rlang_1.3.0 KEGGREST_1.53.6