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

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. 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")
spicyTestPair
#> [1] 2
topPairs(spicyTestPair)
#>            intercept coefficient      p.value   adj.pvalue from    to
#> beta_delta  97.19869   -31.24584 0.0001555308 0.0001555308 beta delta

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

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.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 Tc

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

5 sessionInfo()

sessionInfo()
#> R version 4.0.4 (2021-02-15)
#> 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] ggplot2_3.3.3       S4Vectors_0.28.1    BiocGenerics_0.36.0
#> [4] spicyR_1.2.1        BiocStyle_2.18.1   
#> 
#> loaded via a namespace (and not attached):
#>  [1] sass_0.3.1            tidyr_1.1.3           jsonlite_1.7.2       
#>  [4] splines_4.0.4         bslib_0.2.4           assertthat_0.2.1     
#>  [7] statmod_1.4.35        highr_0.8             BiocManager_1.30.10  
#> [10] spatstat.geom_1.65-5  yaml_2.2.1            numDeriv_2016.8-1.1  
#> [13] pillar_1.5.1          lattice_0.20-41       glue_1.4.2           
#> [16] digest_0.6.27         RColorBrewer_1.1-2    polyclip_1.10-0      
#> [19] minqa_1.2.4           colorspace_2.0-0      htmltools_0.5.1.1    
#> [22] Matrix_1.3-2          spatstat.sparse_2.0-0 pkgconfig_2.0.3      
#> [25] pheatmap_1.0.12       magick_2.7.0          bookdown_0.21        
#> [28] fftwtools_0.9-11      purrr_0.3.4           spatstat.core_1.65-5 
#> [31] scales_1.1.1          tensor_1.5            spatstat.utils_2.1-0 
#> [34] BiocParallel_1.24.1   lme4_1.1-26           tibble_3.1.0         
#> [37] mgcv_1.8-34           farver_2.1.0          generics_0.1.0       
#> [40] IRanges_2.24.1        ellipsis_0.3.1        withr_2.4.1          
#> [43] magrittr_2.0.1        crayon_1.4.1          deldir_0.2-10        
#> [46] evaluate_0.14         fansi_0.4.2           nlme_3.1-152         
#> [49] MASS_7.3-53.1         tools_4.0.4           data.table_1.14.0    
#> [52] lifecycle_1.0.0       stringr_1.4.0         munsell_0.5.0        
#> [55] compiler_4.0.4        jquerylib_0.1.3       concaveman_1.1.0     
#> [58] rlang_0.4.10          debugme_1.1.0         grid_4.0.4           
#> [61] nloptr_1.2.2.2        goftest_1.2-2         labeling_0.4.2       
#> [64] rmarkdown_2.7         boot_1.3-27           gtable_0.3.0         
#> [67] lmerTest_3.1-3        abind_1.4-5           DBI_1.1.1            
#> [70] R6_2.5.0              knitr_1.31            dplyr_1.0.5          
#> [73] utf8_1.2.1            stringi_1.5.3         spatstat.data_2.0-0  
#> [76] Rcpp_1.0.6            vctrs_0.3.6           rpart_4.1-15         
#> [79] tidyselect_1.1.0      xfun_0.22