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 melanoma image dataset with two conditions: Responders and Non-Responders. With this data set, we want to see if there are differences in how cell types are localised with respect to each other.

cells is a SegmentedCells object containing single-cell data of 135 images from 27 subjects, with 5 images per subjects. There are 9 Non-Responder subjects and 18 Responder 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 DataFrame object providing the corresponding subject (subject) and condition (condition) for each image.

data("melanomaResponders")
melanomaResponders
cellSummary(melanomaResponders)
#> DataFrame with 296523 rows and 6 columns
#>         imageID      cellID imageCellID         x         y      cellType
#>        <factor> <character> <character> <numeric> <numeric>      <factor>
#> 1             1      cell_1      cell_1      49.5      1.25 CD8-PD1-PDL1-
#> 2             1      cell_2      cell_2      85.0      1.00 CD8-PD1-PDL1-
#> 3             1      cell_3      cell_3     169.0      1.25 SOX10+       
#> 4             1      cell_4      cell_4     245.0      1.00 CD8-PD1-PDL1-
#> 5             1      cell_5      cell_5     291.0      1.00 SOX10+       
#> ...         ...         ...         ...       ...       ...           ...
#> 296519      525 cell_351307 cell_351307      4.75    248.75 CD8-PD1-PDL1-
#> 296520      525 cell_351308 cell_351308     80.50    249.00 SOX10+       
#> 296521      525 cell_351309 cell_351309    255.75    249.00 CD8-PD1-PDL1-
#> 296522      525 cell_351310 cell_351310    291.50    248.75 SOX10+       
#> 296523      525 cell_351311 cell_351311    310.25    248.75 SOX10+
imagePheno(melanomaResponders)
#> DataFrame with 135 rows and 3 columns
#>      imageID      condition    subject
#>     <factor>       <factor>   <factor>
#> 1          1 Non-Responders SP12-18569
#> 2          2 Non-Responders SP12-18569
#> 3          3 Non-Responders SP12-18569
#> 4          4 Non-Responders SP12-18569
#> 5          5 Non-Responders SP12-18569
#> ...      ...            ...        ...
#> 131      521     Responders SW14-18104
#> 132      522     Responders SW14-18104
#> 133      523     Responders SW14-18104
#> 134      524     Responders SW14-18104
#> 135      525     Responders SW14-18104

In this data set, cellType can be some combination of the expression of CD8, PD1 or PDL1, or is SOX10+.

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 CD8-PD1+PDL1+ cells (to) thend to be found around CD8+PD1+PDL1- cells (from).

spicyTestPair <- spicy(melanomaResponders, 
                   condition = "condition", 
                   subject = "subject", 
                   from = "CD8+PD1+PDL1-", 
                   to = "CD8-PD1+PDL1+")
#> Calculating pairwise spatial associations
#> Testing for spatial differences across conditions accounting for multiple images per subject
spicyTestPair
#> 
#> Number of cell type pairs: 1
#> Number of differentially localised cell type pairs:
#> [1] 0
topPairs(spicyTestPair)
#>                             intercept coefficient    p.value adj.pvalue
#> CD8+PD1+PDL1-_CD8-PD1+PDL1+  12.30061   -6.676495 0.07270545 0.07270545
#>                                      from            to
#> CD8+PD1+PDL1-_CD8-PD1+PDL1+ CD8+PD1+PDL1- CD8-PD1+PDL1+

We obtain a spicy object which details the results of the mixed effects modelling performed. As the coefficient in spicyTest is negative, we find that CD8-PD1+PDL1+ cells are more likely to be found around CD8+PD1+PDL1- cells in Non-Responders.

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().

spicyTest <- spicy(melanomaResponders, 
                   condition = "condition", 
                   subject = "subject")
spicyTest
#> 
#> Number of cell type pairs: 81
#> Number of differentially localised cell type pairs:
#> conditionResponders 
#>                   0
topPairs(spicyTest)  
#>                             intercept coefficient     p.value adj.pvalue
#> CD8-PD1+PDL1+_CD8+PD1+PDL1- 21.403625  -14.143257 0.006696371  0.3184453
#> CD8+PD1+PDL1-_CD8-PD1+PDL1+ 21.300751  -13.854313 0.007862848  0.3184453
#> CD8-PD1+PDL1-_CD8+PD1-PDL1+  6.348714   -7.144476 0.076475111  0.9174987
#> CD8-PD1+PDL1-_CD8-PD1+PDL1+ 12.485009   -7.625312 0.081119183  0.9174987
#> CD8+PD1-PDL1-_CD8+PD1-PDL1- 20.934811  -11.403936 0.093658867  0.9174987
#> CD8-PD1+PDL1+_CD8-PD1+PDL1- 11.621590   -6.942920 0.112733933  0.9174987
#> CD8+PD1-PDL1+_CD8-PD1+PDL1-  4.821599   -6.028484 0.119793885  0.9174987
#> CD8+PD1-PDL1-_CD8+PD1+PDL1- 11.049319    7.124038 0.123391930  0.9174987
#> CD8-PD1-PDL1-_CD8+PD1-PDL1-  9.067690   -2.765159 0.141393584  0.9174987
#> CD8+PD1+PDL1-_CD8+PD1-PDL1- 11.041912    6.860201 0.143290192  0.9174987
#>                                      from            to
#> 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-
#> CD8+PD1+PDL1-_CD8+PD1-PDL1- CD8+PD1+PDL1- CD8+PD1-PDL1-

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))

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().

spicyTestBootstrap <- spicy(melanomaResponders, 
                   condition = "condition", 
                   subject = "subject", 
                   nsim = 199)
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-

signifPlot(spicyTestBootstrap, breaks=c(-3, 3, 0.5))

Indeed, we get improved statistical power compared to the previous method.

5 sessionInfo()

sessionInfo()
#> R version 4.0.2 (2020-06-22)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 18.04.4 LTS
#> 
#> Matrix products: default
#> BLAS:   /home/biocbuild/bbs-3.11-bioc/R/lib/libRblas.so
#> LAPACK: /home/biocbuild/bbs-3.11-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.26.1    BiocGenerics_0.34.0 ggplot2_3.3.2      
#> [4] spicyR_1.0.1        BiocStyle_2.16.0   
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.5            lattice_0.20-41       deldir_0.1-25        
#>  [4] class_7.3-17          digest_0.6.25         V8_3.2.0             
#>  [7] R6_2.4.1              evaluate_0.14         tensor_1.5           
#> [10] pillar_1.4.6          rlang_0.4.7           curl_4.3             
#> [13] minqa_1.2.4           nloptr_1.2.2.2        magick_2.4.0         
#> [16] rpart_4.1-15          Matrix_1.2-18         goftest_1.2-2        
#> [19] rmarkdown_2.3         labeling_0.3          splines_4.0.2        
#> [22] lme4_1.1-23           BiocParallel_1.22.0   statmod_1.4.34       
#> [25] stringr_1.4.0         pheatmap_1.0.12       polyclip_1.10-0      
#> [28] munsell_0.5.0         spatstat.data_1.4-3   compiler_4.0.2       
#> [31] numDeriv_2016.8-1.1   xfun_0.15             pkgconfig_2.0.3      
#> [34] lmerTest_3.1-2        mgcv_1.8-31           htmltools_0.5.0      
#> [37] tidyselect_1.1.0      tibble_3.0.3          bookdown_0.20        
#> [40] IRanges_2.22.2        crayon_1.3.4          dplyr_1.0.0          
#> [43] withr_2.2.0           MASS_7.3-51.6         grid_4.0.2           
#> [46] jsonlite_1.7.0        nlme_3.1-148          gtable_0.3.0         
#> [49] lifecycle_0.2.0       magrittr_1.5          concaveman_1.1.0     
#> [52] scales_1.1.1          stringi_1.4.6         farver_2.0.3         
#> [55] spatstat_1.64-1       ellipsis_0.3.1        spatstat.utils_1.17-0
#> [58] generics_0.0.2        vctrs_0.3.1           boot_1.3-25          
#> [61] RColorBrewer_1.1-2    tools_4.0.2           glue_1.4.1           
#> [64] purrr_0.3.4           abind_1.4-5           yaml_2.2.1           
#> [67] colorspace_1.4-1      BiocManager_1.30.10   knitr_1.29