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 rows

In 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 Th

We 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          alpha

Again, 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