Clustering local indicators of spatial association (LISA) functions
is a methodology for identifying consistent spatial organisation of
multiple cell-types in an unsupervised way. This can be used to enable
the characterization of interactions between multiple cell-types
simultaneously and can complement traditional pairwise analysis. In our
implementation our LISA curves are a localised summary of an L-function
from a Poisson point process model. Our framework lisaClust
can be used to provide a high-level summary of cell-type colocalization
in high-parameter spatial cytometry data, facilitating the
identification of distinct tissue compartments or identification of
complex cellular microenvironments.
TO illustrate our lisaClust
framework, here we consider
a very simple toy example where two cell-types are completely separated
spatially. We simulate data for two different images.
set.seed(51773)
x <- round(c(
runif(200), runif(200) + 1, runif(200) + 2, runif(200) + 3,
runif(200) + 3, runif(200) + 2, runif(200) + 1, runif(200)
), 4) * 100
y <- round(c(
runif(200), runif(200) + 1, runif(200) + 2, runif(200) + 3,
runif(200), runif(200) + 1, runif(200) + 2, runif(200) + 3
), 4) * 100
cellType <- factor(paste("c", rep(rep(c(1:2), rep(200, 2)), 4), sep = ""))
imageID <- rep(c("s1", "s2"), c(800, 800))
cells <- data.frame(x, y, cellType, imageID)
ggplot(cells, aes(x, y, colour = cellType)) +
geom_point() +
facet_wrap(~imageID) +
theme_minimal()
First we store our data in a SingleCellExperiment
object.
We can then use the convenience function lisaClust
to
simultaneously calculate local indicators of spatial association (LISA)
functions using the lisa
function and perform k-means
clustering. The number of clusters can be specified with the
k =
parameter. In the example below, we’ve chosen
k = 2
, resulting in a total of 2 clusters.
These clusters are stored in colData
of the
SingleCellExperiment
object, as a new column with the
column name regions
.
SCE <- lisaClust(SCE, k = 2)
colData(SCE) |> head()
## DataFrame with 6 rows and 5 columns
## x y cellType imageID region
## <numeric> <numeric> <factor> <character> <character>
## 1 36.72 38.58 c1 s1 region_2
## 2 61.38 41.29 c1 s1 region_2
## 3 33.59 80.98 c1 s1 region_2
## 4 50.17 64.91 c1 s1 region_2
## 5 82.93 35.60 c1 s1 region_2
## 6 83.13 2.69 c1 s1 region_2
lisaClust
also provides the convenient
hatchingPlot
function to visualise the different regions
that have been demarcated by the clustering. hatchingPlot
outputs a ggplot
object where the regions are marked by
different hatching patterns. In a real biological dataset, this allows
us to plot both regions and cell-types on the same visualization.
In the example below, we can visualise our stimulated data where our
2 cell types have been separated neatly into 2 distinct regions based on
which cell type each region is dominated by. region_2
is
dominated by the red cell type c1
, and
region_1
is dominated by the blue cell type
c2
.
## Using other clustering methods.
While the lisaClust
function is convenient, we have not
implemented an exhaustive suite of clustering methods as it is very easy
to do this yourself. There are just two simple steps.
We can calculate local indicators of spatial association (LISA)
functions using the lisa
function. Here the LISA curves are
a localised summary of an L-function from a Poisson point process model.
The radii that will be calculated over can be set with
Rs
.
The LISA curves can then be used to cluster the cells. Here we use
k-means clustering, other clustering methods like SOM could be used. We
can store these cell clusters or cell “regions” in our
SingleCellExperiment
object.
# Custom clustering algorithm
kM <- kmeans(lisaCurves, 2)
# Storing clusters into colData
colData(SCE)$custom_region <- paste("region", kM$cluster, sep = "_")
colData(SCE) |> head()
## DataFrame with 6 rows and 6 columns
## x y cellType imageID region custom_region
## <numeric> <numeric> <factor> <character> <character> <character>
## 1 36.72 38.58 c1 s1 region_2 region_1
## 2 61.38 41.29 c1 s1 region_2 region_1
## 3 33.59 80.98 c1 s1 region_2 region_1
## 4 50.17 64.91 c1 s1 region_2 region_1
## 5 82.93 35.60 c1 s1 region_2 region_1
## 6 83.13 2.69 c1 s1 region_2 region_1
Next, we apply our lisaClust
framework to three images
of pancreatic islets from A Map of Human Type 1 Diabetes Progression
by Imaging Mass Cytometry by Damond et al. (2019).
We will start by reading in the data and storing it as a
SingleCellExperiment
object. Here the data is in a format
consistent with that outputted by CellProfiler.
isletFile <- system.file("extdata", "isletCells.txt.gz", package = "spicyR")
cells <- read.table(isletFile, header = TRUE)
damonSCE <- SingleCellExperiment(
assay = list(intensities = t(cells[, grepl(names(cells), pattern = "Intensity_")])),
colData = cells[, !grepl(names(cells), pattern = "Intensity_")]
)
This data does not include annotation of the cell-types of each cell.
Here we extract the marker intensities from the
SingleCellExperiment
object using assay()
. We
then perform k-means clustering with 10 clusters and store these
cell-type clusters in our SingleCellExperiment
object using
colData() <-
.
markers <- t(assay(damonSCE, "intensities"))
kM <- kmeans(markers, 10)
colData(damonSCE)$cluster <- paste("cluster", kM$cluster, sep = "")
colData(damonSCE)[, c("ImageNumber", "cluster")] |> head()
## DataFrame with 6 rows and 2 columns
## ImageNumber cluster
## <integer> <character>
## 1 1 cluster10
## 2 1 cluster10
## 3 1 cluster10
## 4 1 cluster3
## 5 1 cluster3
## 6 1 cluster1
As before, we can perform k-means clustering on the local indicators
of spatial association (LISA) functions using the lisaClust
function, remembering to specify the imageID
,
cellType
, and spatialCoords
columns in
colData
.
damonSCE <- lisaClust(damonSCE,
k = 2,
Rs = c(10, 20, 50),
imageID = "ImageNumber",
cellType = "cluster",
spatialCoords = c("Location_Center_X", "Location_Center_Y")
)
These regions are stored in colData
and can be
extracted.
colData(damonSCE)[, c("ImageNumber", "region")] |>
head(20)
## DataFrame with 20 rows and 2 columns
## ImageNumber region
## <integer> <character>
## 1 1 region_2
## 2 1 region_2
## 3 1 region_2
## 4 1 region_2
## 5 1 region_2
## ... ... ...
## 16 1 region_2
## 17 1 region_2
## 18 1 region_2
## 19 1 region_2
## 20 1 region_2
lisaClust
also provides a convenient function,
regionMap
, for examining which cell types are located in
which regions. In this example, we use this to check which cell types
appear more frequently in each region than expected by chance.
Here, we clearly see that clusters 2, 5, 1, and 8 are highly concentrated in region 1, whilst all other clusters are thinly spread out across region 2.
We can further segregate these cells by increasing the number of
clusters, ie. increasing the parameter k =
in the
lisaClust()
function, but for the purposes of
demonstration, let’s take a look at the hatchingPlot
of
these regions.
regionMap(damonSCE,
imageID = "ImageNumber",
cellType = "cluster",
spatialCoords = c("Location_Center_X", "Location_Center_Y"),
type = "bubble"
)
Finally, we can use hatchingPlot
to construct a
ggplot
object where the regions are marked by different
hatching patterns. This allows us to visualize the two regions and ten
cell-types simultaneously.
hatchingPlot(damonSCE,
imageID = "ImageNumber",
cellType = "cluster",
spatialCoords = c("Location_Center_X", "Location_Center_Y")
)
sessionInfo()
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04 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.26.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=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: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] SingleCellExperiment_1.27.2 SummarizedExperiment_1.35.1
## [3] Biobase_2.65.0 GenomicRanges_1.57.1
## [5] GenomeInfoDb_1.41.1 IRanges_2.39.1
## [7] S4Vectors_0.43.1 BiocGenerics_0.51.0
## [9] MatrixGenerics_1.17.0 matrixStats_1.3.0
## [11] ggplot2_3.5.1 spicyR_1.17.1
## [13] lisaClust_1.13.1 BiocStyle_2.33.1
##
## loaded via a namespace (and not attached):
## [1] deldir_2.0-4 rlang_1.1.4
## [3] magrittr_2.0.3 compiler_4.4.1
## [5] spatstat.geom_3.2-9 mgcv_1.9-1
## [7] vctrs_0.6.5 reshape2_1.4.4
## [9] stringr_1.5.1 pkgconfig_2.0.3
## [11] SpatialExperiment_1.15.1 crayon_1.5.3
## [13] fastmap_1.2.0 backports_1.5.0
## [15] magick_2.8.3 XVector_0.45.0
## [17] labeling_0.4.3 utf8_1.2.4
## [19] rmarkdown_2.27 nloptr_2.1.1
## [21] UCSC.utils_1.1.0 purrr_1.0.2
## [23] xfun_0.45 MultiAssayExperiment_1.31.4
## [25] zlibbioc_1.51.1 cachem_1.1.0
## [27] jsonlite_1.8.8 goftest_1.2-3
## [29] highr_0.11 DelayedArray_0.31.7
## [31] tweenr_2.0.3 spatstat.utils_3.0-5
## [33] BiocParallel_1.39.0 broom_1.0.6
## [35] parallel_4.4.1 R6_2.5.1
## [37] stringi_1.8.4 bslib_0.7.0
## [39] RColorBrewer_1.1-3 spatstat.data_3.1-2
## [41] boot_1.3-30 car_3.1-2
## [43] numDeriv_2016.8-1.1 jquerylib_0.1.4
## [45] Rcpp_1.0.12 knitr_1.48
## [47] tensor_1.5 splines_4.4.1
## [49] Matrix_1.7-0 tidyselect_1.2.1
## [51] abind_1.4-5 yaml_2.3.9
## [53] codetools_0.2-20 spatstat.random_3.2-3
## [55] spatstat.explore_3.2-7 curl_5.2.1
## [57] lmerTest_3.1-3 lattice_0.22-6
## [59] tibble_3.2.1 plyr_1.8.9
## [61] withr_3.0.0 evaluate_0.24.0
## [63] survival_3.7-0 polyclip_1.10-6
## [65] ggupset_0.4.0 pillar_1.9.0
## [67] BiocManager_1.30.23 ggpubr_0.6.0
## [69] carData_3.0-5 ClassifyR_3.9.2
## [71] generics_0.1.3 munsell_0.5.1
## [73] scales_1.3.0 minqa_1.2.7
## [75] class_7.3-22 glue_1.7.0
## [77] pheatmap_1.0.12 maketools_1.3.0
## [79] tools_4.4.1 sys_3.4.2
## [81] data.table_1.15.4 lme4_1.1-35.5
## [83] ggsignif_0.6.4 buildtools_1.0.0
## [85] grid_4.4.1 tidyr_1.3.1
## [87] colorspace_2.1-0 nlme_3.1-165
## [89] GenomeInfoDbData_1.2.12 ggforce_0.4.2
## [91] cli_3.6.3 spatstat.sparse_3.1-0
## [93] fansi_1.0.6 S4Arrays_1.5.4
## [95] scam_1.2-17 dplyr_1.1.4
## [97] concaveman_1.1.0 V8_4.4.2
## [99] gtable_0.3.5 rstatix_0.7.2
## [101] sass_0.4.9 digest_0.6.36
## [103] SparseArray_1.5.18 farver_2.1.2
## [105] rjson_0.2.21 htmltools_0.5.8.1
## [107] lifecycle_1.0.4 httr_1.4.7
## [109] MASS_7.3-61