Perform linear and mixed-effects models to assess and visualise changes in cell localisation across disease conditions.
spicyR 1.19.3
if (!require("BiocManager")) {
install.packages("BiocManager")
}
BiocManager::install("spicyR")# load required packages
library(spicyR)
library(ggplot2)
library(SpatialExperiment)
library(SpatialDatasets)
library(imcRtools)
library(dplyr)
library(survival)This guide provides step-by-step instructions on how to apply a linear model to multiple segmented and labelled images to assess how the localisation of different cell types changes across different disease conditions.
We use the Keren et al. (2018) breast cancer dataset to compare the spatial distribution of immune cells in individuals with different levels of tumour infiltration (cold and compartmentalised).
The data is stored as a SpatialExperiment object and contains single-cell spatial data from 41 images.
kerenSPE <- SpatialDatasets::spe_Keren_2018()The cell types in this dataset includes 11 immune cell types (double negative CD3 T cells, CD4 T cells, B cells, monocytes, macrophages, CD8 T cells, neutrophils, natural killer cells, dendritic cells, regulatory T cells), 2 structural cell types (endothelial, mesenchymal), 2 tumour cell types (keratin+ tumour, tumour) and one unidentified category.
To investigate changes in localisation between two different cell types, we measure the level of localisation between two cell types by modelling with the L-function. The L-function is a variance-stabilised K-function given by the equation
\[ \widehat{L_{ij}} (r) = \sqrt{\frac{\widehat{K_{ij}}(r)}{\pi}} \]
with \(\widehat{K_{ij}}\) defined as
\[ \widehat{K_{ij}} (r) = \frac{|W|}{n_i n_j} \sum_{n_i} \sum_{n_j} 1 \{d_{ij} \leq r \} e_{ij} (r) \]
where \(\widehat{K_{ij}}\) summarises the degree of co-localisation of cell type \(j\) with cell type \(i\), \(n_i\) and \(n_j\) are the number of cells of type \(i\) and \(j\), \(|W|\) is the image area, \(d_{ij}\) is the distance between two cells and \(e_{ij} (r)\) is an edge correcting factor.
Specifically, the mean difference between the experimental function and the theoretical function is used as a measure for the level of localisation, defined as
\[ u = \sum_{r' = r_{\text{min}}}^{r_{\text{max}}} \widehat L_{ij, \text{Experimental}} (r') - \widehat L_{ij, \text{Poisson}} (r') \]
where \(u\) is the sum is taken over a discrete range of \(r\) between \(r_{\text{min}}\) and \(r_{\text{max}}\). Differences of the statistic \(u\) between two conditions is modelled using a weighted linear model.
Firstly, we can test whether one cell type tends to be more localised with another cell type in one condition compared to the other. This can be done using the spicy() function, where we specify the condition parameter.
In this example, we want to see whether or not neutrophils (to) tend to be found around CD8 T cells (from) in compartmentalised tumours compared to cold tumours. Given that there are 3 conditions, we can specify the desired conditions by setting the order of our condition factor. spicy will choose the first level of the factor as the base condition and the second level as the comparison condition. spicy will also naturally coerce the condition column into a factor if it is not already a factor. The column containing cell type annotations can be specified using the cellType argument. By default, spicy uses the column named cellType in the SpatialExperiment object.
spicyTestPair <- spicy(
kerenSPE,
condition = "tumour_type",
from = "CD8_T_cell",
to = "Neutrophils"
)
topPairs(spicyTestPair)
#> intercept coefficient p.value adj.pvalue
#> CD8_T_cell__Neutrophils -109.081 112.0185 2.166646e-05 2.166646e-05
#> from to
#> CD8_T_cell__Neutrophils CD8_T_cell NeutrophilsWe obtain a spicy object which details the results of the
modelling performed. The topPairs() function can be used to obtain the associated coefficients and p-value.
As the coefficient in spicyTestPair is positive, we find that neutrophils are significantly more likely to be found near CD8 T cells in the compartmentalised tumours group compared to the cold tumour group.
We can perform what we did above for all pairwise combinations of cell
types by excluding the from and to parameters in spicy().
spicyTest <- spicy(
kerenSPE,
condition = "tumour_type"
)
topPairs(spicyTest)
#> intercept coefficient p.value adj.pvalue
#> Macrophages__dn_T_CD3 56.446064 -50.08474 1.080273e-07 3.035568e-05
#> dn_T_CD3__Macrophages 54.987151 -48.38664 2.194018e-07 3.082595e-05
#> Macrophages__DC_or_Mono 73.239404 -59.90361 5.224660e-06 4.893765e-04
#> DC_or_Mono__Macrophages 71.777087 -58.46833 7.431172e-06 5.220399e-04
#> dn_T_CD3__dn_T_CD3 -63.786032 100.61010 2.878804e-05 1.208706e-03
#> Neutrophils__dn_T_CD3 -63.141840 69.64356 2.891872e-05 1.208706e-03
#> dn_T_CD3__Neutrophils -63.133725 70.15508 3.011012e-05 1.208706e-03
#> DC__Macrophages 96.893239 -92.55112 1.801300e-04 5.758129e-03
#> Macrophages__DC 96.896215 -93.25194 1.844241e-04 5.758129e-03
#> CD4_T_cell__Keratin_Tumour -4.845037 -22.14995 2.834659e-04 7.409016e-03
#> from to
#> Macrophages__dn_T_CD3 Macrophages dn_T_CD3
#> dn_T_CD3__Macrophages dn_T_CD3 Macrophages
#> Macrophages__DC_or_Mono Macrophages DC_or_Mono
#> DC_or_Mono__Macrophages DC_or_Mono Macrophages
#> dn_T_CD3__dn_T_CD3 dn_T_CD3 dn_T_CD3
#> Neutrophils__dn_T_CD3 Neutrophils dn_T_CD3
#> dn_T_CD3__Neutrophils dn_T_CD3 Neutrophils
#> DC__Macrophages DC Macrophages
#> Macrophages__DC Macrophages DC
#> CD4_T_cell__Keratin_Tumour CD4_T_cell Keratin_TumourAgain, we obtain a spicy object which outlines the result of the linear
models performed for each pairwise combination of cell types.
We can also examine the L-function metrics of individual images by using the
convenient bind() function on our spicyTest results object.
bind(spicyTest)[1:5, 1:5]
#> imageID condition Keratin_Tumour__Keratin_Tumour
#> 1 1 mixed -2.300602
#> 2 2 mixed -1.989699
#> 3 3 compartmentalised 11.373530
#> 4 4 compartmentalised 33.931133
#> 5 5 compartmentalised 17.922818
#> dn_T_CD3__Keratin_Tumour B_cell__Keratin_Tumour
#> 1 -5.298543 -20.827279
#> 2 -16.020022 3.025815
#> 3 -21.857447 -24.962913
#> 4 -36.438476 -40.470221
#> 5 -20.816783 -38.138076The results can be represented as a bubble plot using the signifPlot() function.
signifPlot(
spicyTest,
breaks = c(-3, 3, 1),
marksToPlot = c("Macrophages", "DC_or_Mono", "dn_T_CD3", "Neutrophils",
"CD8_T_cell", "Keratin_Tumour")
)