Spatial Mixed-Effects Modelling with spicy
26 October 2021
Package
spicyR 1.6.0
Contents
# load required packages
library(spicyR)
library(ggplot2)1 Installation
if (!require("BiocManager"))
install.packages("BiocManager")
BiocManager::install("spicyR")2 Overview
This guide will provide a step-by-step guide on how mixed effects models can be applied to multiple segmented and labelled images to identify how the localisation of different cell types can change across different conditions. Here, the subject is modelled as a random effect, and the different conditions are modelled as a fixed effect.
3 Example data
Here, we use a subset of the Damond et al 2019 imaging mass cytometry dataset. We will compare the spatial distributions of cells in the pancreatic islets of individuals with early onset diabetes and healthy controls.
diabetesData is a SegmentedCells object containing single-cell data of 160 images
from 8 subjects, with 20 images per subjects.
cellSummary() returns a DataFrame object providing the location (x and y)
and cell type (cellType) of each cell and the image it belongs to (imageID).
imagePheno() returns a tibble object providing the corresponding subject
(subject) and condition (condition) for each image.
data("diabetesData")
diabetesData
#> A SegmentedCells object with...
#> Number of images:120
#> Number of cells:253777
#> Number of cell types: 16 [ ductal, acinar, ..., beta ]
#> Number of intensities: 0 [ ]
#> Number of morphologies: 0 [ ]
#> Number of image phenotypes: 6 [ imageID, case, ..., stage ]
cellSummary(diabetesData)
#> DataFrame with 253777 rows and 6 columns
#> imageID cellID imageCellID x y cellType
#> <factor> <character> <character> <numeric> <numeric> <factor>
#> 1 P15 cell_1101219 P15_1 219 2 ductal
#> 2 P15 cell_1101220 P15_2 277 0 acinar
#> 3 P15 cell_1101221 P15_3 22 0 ductal
#> 4 P15 cell_1101222 P15_4 30 6 ductal
#> 5 P15 cell_1101223 P15_5 105 0 acinar
#> ... ... ... ... ... ... ...
#> 253773 D25 cell_295062 D25_1016 304 322 ductal
#> 253774 D25 cell_295063 D25_1017 230 325 ductal
#> 253775 D25 cell_295064 D25_1018 328 325 acinar
#> 253776 D25 cell_295065 D25_1019 94 325 Tc
#> 253777 D25 cell_295066 D25_1020 268 325 acinar
imagePheno(diabetesData)
#> DataFrame with 120 rows and 6 columns
#> imageID case slide part group stage
#> <factor> <integer> <character> <character> <integer> <factor>
#> 1 P15 6089 P Tail 3 Long-duration
#> 2 Q06 6089 Q Body 3 Long-duration
#> 3 Q20 6089 Q Body 3 Long-duration
#> 4 P36 6089 P Tail 3 Long-duration
#> 5 P34 6089 P Tail 3 Long-duration
#> ... ... ... ... ... ... ...
#> 116 D03 6418 D Body 1 Long-duration
#> 117 C13 6418 C Tail 1 Long-duration
#> 118 D11 6418 D Body 1 Long-duration
#> 119 D06 6418 D Body 1 Long-duration
#> 120 D25 6418 D Body 1 Long-durationIn this data set, cell types include immune cell types (B cells, naive T cells, T Helper cells, T cytotoxic cells, neutrophils, macrophages) and pancreatic islet cells (alpha, beta, gamma, delta).
4 Mixed Effects Modelling
To investigate changes in colocalisation between two different cell types, we
measure the level of colocalisation between two cell types by modelling with the
Lcross() function in the spatstat package. Specifically, the mean difference
between the obtained function and the theoretical function is used as a measure
for the level of colocalisation. Differences of this statistic between two
conditions is modelled using a weighted mixed effects model, with condition as
the fixed effect and subject as the random effect.
spicyTestBootstrap
## Testing for change in colocalisation for a specific pair of cells
Firstly, we can see whether one cell type tends to be around another cell type
in one condition compared to the other. This can be done using the spicy()
function, where we include condition, and subject. In this example, we want
to see whether or not Delta cells (to) tend to be found around Beta cells (from)
in onset diabetes images compared to non-diabetic images.
spicyTestPair <- spicy(diabetesData,
condition = "stage",
subject = "case",
from = "beta",
to = "delta")
topPairs(spicyTestPair)
#> intercept coefficient p.value adj.pvalue from to
#> beta__delta 184.3982 -2.141889 0.8677429 0.8677429 beta deltaWe obtain a spicy object which details the results of the mixed effects
modelling performed. As the coefficient in spicyTest is positive, we find
that Th cells cells are more likely to be found around beta cells in the onset
diabetes group compared to the non-diabetic control.
4.1 Test for change in colocalisation for all pairwise cell combinations
Here, we can perform what we did above for all pairwise combinations of cell
types by excluding the from and to parameters from spicy().
spicyTest <- spicy(diabetesData,
condition = "stage",
subject = "case")data("spicyTest")spicyTest
#>
#> Number of cell type pairs: 256
#> Number of differentially localised cell type pairs:
#> conditionOnset conditionLong-duration
#> 7 7
topPairs(spicyTest)
#> intercept coefficient p.value adj.pvalue from
#> Th_Th -8.550349 55.505421 0.0002069040 0.03207907 Th
#> beta_ductal -19.461966 10.416870 0.0005559049 0.03207907 beta
#> ductal_beta -19.407845 10.382899 0.0005695262 0.03207907 ductal
#> delta_beta 101.111037 -30.418789 0.0006118914 0.03207907 delta
#> beta_delta 101.052488 -30.370186 0.0006265444 0.03207907 beta
#> unknown_beta -6.704261 8.546938 0.0010115597 0.03891502 unknown
#> beta_unknown -6.686305 8.505666 0.0010640825 0.03891502 beta
#> gamma_neutrophil -3.234976 -19.102452 0.0229478395 0.39078920 gamma
#> neutrophil_gamma -3.228265 -19.065803 0.0230795918 0.39078920 neutrophil
#> ductal_Th -1.454500 -10.508455 0.0268678016 0.39078920 ductal
#> to
#> Th_Th Th
#> beta_ductal ductal
#> ductal_beta beta
#> delta_beta beta
#> beta_delta delta
#> unknown_beta beta
#> beta_unknown unknown
#> gamma_neutrophil neutrophil
#> neutrophil_gamma gamma
#> ductal_Th ThAgain, we obtain a spicy object which outlines the result of the mixed effects
models performed for each pairwise combination if cell types.
We can represent this as a heatmap using the spatialMEMMultiPlot() function by
providing it the spicy object obtained.
signifPlot(spicyTest,
breaks=c(-3, 3, 0.5),
marksToPlot = c("alpha", "beta", "gamma", "delta",
"B", "naiveTc", "Th", "Tc", "neutrophil", "macrophage"))
4.2 Bootstrapping with spicy
There are multiple ways for calculating p-values for mixed effects models. We
have also implemented a bootstrapping approach. All that is needed is a choice
for the number of resamples used in the bootstrap which can be set with the
nsim parameter in spicy().
data("spicyTestBootstrap")spicyTestBootstrap <- spicy(diabetesData,
condition = "stage",
subject = "case",
from = "beta",
to = "Tc",
nsim = 1000)spicyTestBootstrap
#>
#> Number of cell type pairs: 1
#> Number of differentially localised cell type pairs:
#> [1] 2
topPairs(spicyTestBootstrap)
#> intercept coefficient p.value adj.pvalue from to
#> beta_Tc -16.63134 11.28218 0.014 0.014 beta TcIndeed, we get improved statistical power compared to the previous method.
5 sessionInfo()
sessionInfo()
#> R version 4.1.1 (2021-08-10)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 20.04.3 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.14-bioc/R/lib/libRblas.so
#> LAPACK: /home/biocbuild/bbs-3.14-bioc/R/lib/libRlapack.so
#>
#> 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
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] ggplot2_3.3.5 S4Vectors_0.32.0 BiocGenerics_0.40.0
#> [4] spicyR_1.6.0 BiocStyle_2.22.0
#>
#> loaded via a namespace (and not attached):
#> [1] sass_0.4.0 tidyr_1.1.4 jsonlite_1.7.2
#> [4] splines_4.1.1 bslib_0.3.1 assertthat_0.2.1
#> [7] highr_0.9 BiocManager_1.30.16 spatstat.geom_2.3-0
#> [10] yaml_2.2.1 numDeriv_2016.8-1.1 pillar_1.6.4
#> [13] lattice_0.20-45 glue_1.4.2 digest_0.6.28
#> [16] RColorBrewer_1.1-2 polyclip_1.10-0 minqa_1.2.4
#> [19] colorspace_2.0-2 htmltools_0.5.2 Matrix_1.3-4
#> [22] spatstat.sparse_2.0-0 pkgconfig_2.0.3 pheatmap_1.0.12
#> [25] magick_2.7.3 bookdown_0.24 purrr_0.3.4
#> [28] spatstat.core_2.3-0 scales_1.1.1 tensor_1.5
#> [31] spatstat.utils_2.2-0 BiocParallel_1.28.0 lme4_1.1-27.1
#> [34] tibble_3.1.5 mgcv_1.8-38 farver_2.1.0
#> [37] generics_0.1.1 IRanges_2.28.0 ellipsis_0.3.2
#> [40] withr_2.4.2 magrittr_2.0.1 crayon_1.4.1
#> [43] deldir_1.0-6 evaluate_0.14 fansi_0.5.0
#> [46] nlme_3.1-153 MASS_7.3-54 tools_4.1.1
#> [49] data.table_1.14.2 lifecycle_1.0.1 stringr_1.4.0
#> [52] munsell_0.5.0 scam_1.2-12 compiler_4.1.1
#> [55] jquerylib_0.1.4 concaveman_1.1.0 rlang_0.4.12
#> [58] grid_4.1.1 nloptr_1.2.2.2 goftest_1.2-3
#> [61] labeling_0.4.2 rmarkdown_2.11 boot_1.3-28
#> [64] gtable_0.3.0 lmerTest_3.1-3 abind_1.4-5
#> [67] DBI_1.1.1 R6_2.5.1 knitr_1.36
#> [70] dplyr_1.0.7 fastmap_1.1.0 utf8_1.2.2
#> [73] stringi_1.7.5 spatstat.data_2.1-0 parallel_4.1.1
#> [76] Rcpp_1.0.7 vctrs_0.3.8 rpart_4.1-15
#> [79] tidyselect_1.1.1 xfun_0.27