Spatial Mixed-Effects Modelling with spicy
27 October 2020
Package
spicyR 1.2.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.
diabetesDataSub 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("diabetesDataSub")
diabetesDataSub
cellSummary(diabetesDataSub)
#> DataFrame with 49888 rows and 6 columns
#> imageID cellID imageCellID x y cellType
#> <factor> <character> <character> <numeric> <numeric> <factor>
#> 1 1 cell_1 cell_1 449.5263 1.15789 alpha
#> 2 1 cell_2 cell_2 551.0476 1.47619 alpha
#> 3 1 cell_3 cell_3 38.5067 2.84000 alpha
#> 4 1 cell_4 cell_4 515.8500 1.25000 alpha
#> 5 1 cell_5 cell_5 622.0000 1.41176 alpha
#> ... ... ... ... ... ... ...
#> 49884 564 cell_49884 cell_49884 139.9459 249.439 alpha
#> 49885 564 cell_49885 cell_49885 341.2500 255.100 alpha
#> 49886 564 cell_49886 cell_49886 29.8947 257.553 beta
#> 49887 564 cell_49887 cell_49887 188.4615 258.538 alpha
#> 49888 564 cell_49888 cell_49888 29.2353 272.000 alpha
imagePheno(diabetesDataSub)
#> # A tibble: 160 x 3
#> # Groups: subject [8]
#> imageID condition subject
#> <chr> <fct> <fct>
#> 1 1 Onset 6362
#> 2 2 Onset 6362
#> 3 4 Onset 6362
#> 4 5 Onset 6362
#> 5 9 Onset 6362
#> 6 11 Onset 6362
#> 7 16 Onset 6362
#> 8 17 Onset 6362
#> 9 19 Onset 6362
#> 10 25 Onset 6362
#> # … with 150 more rowsIn 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.
4.1 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 Th cells (to) tend to be found around Tc cells (from).
spicyTestPair <- spicy(diabetesDataSub,
condition = "condition",
subject = "subject",
from = "Tc",
to = "Th")
spicyTestPair
#> [1] 1
topPairs(spicyTestPair)
#> intercept coefficient p.value adj.pvalue from to
#> Tc_Th -8.560209 8.769617 0.02594204 0.02594204 Tc ThWe 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.2 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().
data("spicyTest")spicyTest <- spicy(diabetesDataSub,
condition = "condition",
subject = "subject")spicyTest
#>
#> Number of cell type pairs: 100
#> Number of differentially localised cell type pairs:
#> conditionOnset
#> 3
topPairs(spicyTest)
#> intercept coefficient p.value adj.pvalue from
#> Tc_Tc -5.374072 8.655123 0.0004025104 0.03074952 Tc
#> beta_Th -11.129781 9.584197 0.0006149904 0.03074952 beta
#> Th_beta -10.540529 9.067895 0.0012961801 0.04320600 Th
#> beta_delta 15.044381 -8.371413 0.0035887811 0.07359819 beta
#> delta_beta 14.971690 -8.207731 0.0036799093 0.07359819 delta
#> macrophage_beta -7.548617 6.042619 0.0056501803 0.08544570 macrophage
#> beta_macrophage -7.794805 5.864594 0.0060796672 0.08544570 beta
#> Th_Tc -11.995738 13.286021 0.0068356562 0.08544570 Th
#> Tc_Th -11.722786 13.029044 0.0087099236 0.09677693 Tc
#> alpha_alpha 14.709334 4.047072 0.0185066494 0.18506649 alpha
#> to
#> Tc_Tc Tc
#> beta_Th Th
#> Th_beta beta
#> beta_delta delta
#> delta_beta beta
#> macrophage_beta beta
#> beta_macrophage macrophage
#> Th_Tc Tc
#> Tc_Th Th
#> alpha_alpha alphaAgain, 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.3 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(diabetesDataSub,
condition = "condition",
subject = "subject",
from = "Tc",
to = "Th",
nsim = 1000)spicyTestBootstrap
#>
#> Number of cell type pairs: 81
#> Number of differentially localised cell type pairs:
#> conditionResponders
#> 18
topPairs(spicyTestBootstrap)
#> intercept coefficient p.value adj.pvalue
#> CD8+PD1+PDL1+_SOX10+ -1.849790 -2.277733 0 0
#> CD8+PD1-PDL1-_CD8+PD1+PDL1- 11.049319 7.124038 0 0
#> CD8-PD1+PDL1+_CD8+PD1+PDL1- 21.403625 -14.143257 0 0
#> CD8+PD1+PDL1+_CD8-PD1+PDL1- 6.579050 -2.957959 0 0
#> CD8-PD1+PDL1+_CD8-PD1+PDL1- 11.621590 -6.942920 0 0
#> CD8+PD1-PDL1+_CD8-PD1+PDL1- 4.821599 -6.028484 0 0
#> CD8-PD1-PDL1+_CD8-PD1+PDL1- 4.438305 -2.245794 0 0
#> CD8-PD1-PDL1-_CD8+PD1-PDL1- 9.067690 -2.765159 0 0
#> CD8+PD1+PDL1-_CD8+PD1-PDL1- 11.041912 6.860201 0 0
#> CD8+PD1-PDL1-_CD8+PD1-PDL1- 20.934811 -11.403936 0 0
#> from to
#> CD8+PD1+PDL1+_SOX10+ CD8+PD1+PDL1+ SOX10+
#> CD8+PD1-PDL1-_CD8+PD1+PDL1- CD8+PD1-PDL1- CD8+PD1+PDL1-
#> CD8-PD1+PDL1+_CD8+PD1+PDL1- CD8-PD1+PDL1+ CD8+PD1+PDL1-
#> CD8+PD1+PDL1+_CD8-PD1+PDL1- CD8+PD1+PDL1+ CD8-PD1+PDL1-
#> CD8-PD1+PDL1+_CD8-PD1+PDL1- CD8-PD1+PDL1+ CD8-PD1+PDL1-
#> CD8+PD1-PDL1+_CD8-PD1+PDL1- CD8+PD1-PDL1+ CD8-PD1+PDL1-
#> CD8-PD1-PDL1+_CD8-PD1+PDL1- CD8-PD1-PDL1+ CD8-PD1+PDL1-
#> CD8-PD1-PDL1-_CD8+PD1-PDL1- CD8-PD1-PDL1- CD8+PD1-PDL1-
#> CD8+PD1+PDL1-_CD8+PD1-PDL1- CD8+PD1+PDL1- CD8+PD1-PDL1-
#> CD8+PD1-PDL1-_CD8+PD1-PDL1- CD8+PD1-PDL1- CD8+PD1-PDL1-Indeed, we get improved statistical power compared to the previous method.
5 sessionInfo()
sessionInfo()
#> R version 4.0.3 (2020-10-10)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 18.04.5 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.12-bioc/R/lib/libRblas.so
#> LAPACK: /home/biocbuild/bbs-3.12-bioc/R/lib/libRlapack.so
#>
#> 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
#>
#> attached base packages:
#> [1] parallel stats4 stats graphics grDevices utils datasets
#> [8] methods base
#>
#> other attached packages:
#> [1] S4Vectors_0.28.0 BiocGenerics_0.36.0 ggplot2_3.3.2
#> [4] spicyR_1.2.0 BiocStyle_2.18.0
#>
#> loaded via a namespace (and not attached):
#> [1] Rcpp_1.0.5 lattice_0.20-41 deldir_0.1-29
#> [4] class_7.3-17 utf8_1.1.4 assertthat_0.2.1
#> [7] digest_0.6.27 V8_3.3.1 R6_2.4.1
#> [10] evaluate_0.14 tensor_1.5 pillar_1.4.6
#> [13] rlang_0.4.8 curl_4.3 minqa_1.2.4
#> [16] nloptr_1.2.2.2 magick_2.5.0 rpart_4.1-15
#> [19] Matrix_1.2-18 goftest_1.2-2 rmarkdown_2.5
#> [22] labeling_0.4.2 splines_4.0.3 lme4_1.1-25
#> [25] BiocParallel_1.24.0 statmod_1.4.35 stringr_1.4.0
#> [28] pheatmap_1.0.12 polyclip_1.10-0 munsell_0.5.0
#> [31] spatstat.data_1.4-3 compiler_4.0.3 numDeriv_2016.8-1.1
#> [34] xfun_0.18 pkgconfig_2.0.3 lmerTest_3.1-3
#> [37] mgcv_1.8-33 htmltools_0.5.0 tidyselect_1.1.0
#> [40] tibble_3.0.4 bookdown_0.21 IRanges_2.24.0
#> [43] fansi_0.4.1 crayon_1.3.4 dplyr_1.0.2
#> [46] withr_2.3.0 MASS_7.3-53 grid_4.0.3
#> [49] jsonlite_1.7.1 nlme_3.1-150 gtable_0.3.0
#> [52] lifecycle_0.2.0 magrittr_1.5 concaveman_1.1.0
#> [55] scales_1.1.1 cli_2.1.0 stringi_1.5.3
#> [58] farver_2.0.3 spatstat_1.64-1 ellipsis_0.3.1
#> [61] spatstat.utils_1.17-0 generics_0.0.2 vctrs_0.3.4
#> [64] boot_1.3-25 RColorBrewer_1.1-2 tools_4.0.3
#> [67] glue_1.4.2 purrr_0.3.4 abind_1.4-5
#> [70] yaml_2.2.1 colorspace_1.4-1 BiocManager_1.30.10
#> [73] knitr_1.30