Contents

1 Installation

if (!require("BiocManager")) {
  install.packages("BiocManager")
}
BiocManager::install("spicyR")
# load required packages
library(spicyR)
library(ggplot2)
library(SingleCellExperiment)
library(SpatialExperiment)
library(imcRtools)

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_SCE is a SingleCellExperiment object containing single-cell data of 160 images from 8 subjects, with 20 images per subject.

data("diabetesData_SCE")
diabetesData_SCE
#> class: SingleCellExperiment 
#> dim: 0 253777 
#> metadata(0):
#> assays(0):
#> rownames: NULL
#> rowData names(0):
#> colnames: NULL
#> colData names(11): imageID cellID ... group stage
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):

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 localisation between two different cell types, we measure the level of localisation between two cell types by modelling with the L-function. Specifically, the mean difference between the obtained function and the theoretical function is used as a measure for the level of localisation. 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 localisation for a specific pair of cells

Firstly, we can test whether one cell type tends to be more localised with 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_SCE,
  condition = "stage",
  subject = "case",
  from = "beta",
  to = "delta"
)

topPairs(spicyTestPair)
#>             intercept coefficient     p.value  adj.pvalue from    to
#> beta__delta   179.729   -58.24478 0.000109702 0.000109702 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 delta cells cells are significantly less likely to be found near beta cells in the onset diabetes group compared to the non-diabetic control.

4.2 Test for change in localisation 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_SCE,
  condition = "stage",
  subject = "case"
)

spicyTest
#>         conditionOnset conditionLong-duration 
#>                      0                     15
topPairs(spicyTest)
#>                          intercept coefficient      p.value adj.pvalue
#> beta__delta           1.815458e+02  -48.735693 0.0005033247 0.07169649
#> delta__beta           1.817943e+02  -48.166076 0.0005601288 0.07169649
#> B__unknown            4.327556e-15   11.770938 0.0052338392 0.42051606
#> delta__delta          2.089550e+02  -52.061196 0.0125129422 0.42051606
#> unknown__macrophage   1.007337e+01  -15.826919 0.0207410908 0.42051606
#> unknown__B            0.000000e+00   12.142848 0.0225855404 0.42051606
#> macrophage__unknown   1.004424e+01  -14.471666 0.0244668075 0.42051606
#> B__Th                 3.142958e-15   26.847934 0.0245039854 0.42051606
#> otherimmune__naiveTc -9.292508e+00   33.584755 0.0255812944 0.42051606
#> ductal__ductal        1.481580e+01   -8.632569 0.0266935703 0.42051606
#>                             from         to
#> beta__delta                 beta      delta
#> delta__beta                delta       beta
#> B__unknown                     B    unknown
#> delta__delta               delta      delta
#> unknown__macrophage      unknown macrophage
#> unknown__B               unknown          B
#> macrophage__unknown   macrophage    unknown
#> B__Th                          B         Th
#> otherimmune__naiveTc otherimmune    naiveTc
#> ductal__ductal            ductal     ductal

Again, we obtain a spicy object which outlines the result of the mixed effects models performed for each pairwise combination of cell types.

We can represent this as a bubble plot using the signifPlot() function by providing it the spicy object obtained. Here, we can observe that the most significant relationships occur between pancreatic beta and delta cells, suggesting that the 2 cell types are far more localised during diabetes onset compared to non-diabetics.

signifPlot(
  spicyTest,
  breaks = c(-3, 3, 1),
  marksToPlot = c(
    "alpha", "beta", "gamma", "delta",
    "B", "naiveTc", "Th", "Tc", "neutrophil", "macrophage"
  )
)

If we’re interested and wish to examine a specific cell type-cell type relationship in more detail, we can use spicyBoxPlot, specifying the relationship we want to examine.

In the bubble plot above, we can see that the most significant relationship is between beta and delta islet cells in the pancreas. To examine this further, we can specify either from = beta and to = delta or rank = 1 parameters in spicyBoxPlot.

spicyBoxPlot(results = spicyTest,
             # from = "beta",
             # to = "delta",
             rank = 1)

4.3 Mixed effects modelling for custom metrics

spicyR can also be applied to custom distance or abundance metrics. Here, we provide an example where we apply the spicy function to a custom abundance metric.

We first obtain the custom abundance metric by converting the a SpatialExperiment object from the existing SingleCellExperiment object. A KNN interactions graph is then generated with the function buildSpatialGraph from the imcRtools package. This generates a colPairs object inside of the SpatialExperiment object. spicyR provides the function convPairs for converting a colPairs object stored within a SingleCellExperiment object into an abundance matrix by effectively calculating the average number of nearby cells types for every cell type. For example, if there exists on average 5 neutrophils for every macrophage in image 1, the column neutrophil__macrophage would have a value of 5 for image 1.

spicy can take any input of pairwise cell type combinations across multiple images and run a mixed effects model to determine collective differences across conditions.

diabetesData_SPE <- SpatialExperiment(diabetesData_SCE,
                                      colData = colData(diabetesData_SCE)) 

spatialCoords(diabetesData_SPE) <- data.frame(colData(diabetesData_SPE)$x, colData(diabetesData_SPE)$y) |> as.matrix()
spatialCoordsNames(diabetesData_SPE) <- c("x", "y")

diabetesData_SPE <- imcRtools::buildSpatialGraph(diabetesData_SPE, img_id = "imageID", type = "knn", k = 20, coords = c("x", "y"))
#> 'sample_id's are duplicated across 'SpatialExperiment' objects to cbind; appending sample indices.
#> The returned object is ordered by the 'imageID' entry.

pairAbundances <- convPairs(diabetesData_SPE,
                  colPair = "knn_interaction_graph")

head(pairAbundances["delta__delta"])
#>     delta__delta
#> P15     2.863636
#> Q06     2.217391
#> Q20     8.148148
#> P36     3.076923
#> P34     3.166667
#> Q01     1.230769

spicy can take any input of pairwise cell type combinations across multiple images and run a mixed effects model to determine collective differences across conditions. To check out other custom distance metrics which can be used, feel free to check out the Statial package.

spicyTestColPairs <- spicy(
  diabetesData_SPE,
  condition = "stage",
  subject = "case",
  alternateResult = pairAbundances,
  weights = FALSE
)
#> Cell count weighting set to FALSE for alternate results
#> Testing for spatial differences across conditions accounting for multiple images per subject

topPairs(spicyTestColPairs)
#>                             intercept coefficient     p.value adj.pvalue
#> naiveTc__otherimmune     8.016027e+00 -3.94290293 0.001996075  0.1465276
#> gamma__Th                1.524442e-03  0.01185203 0.002963771  0.1465276
#> macrophage__beta        -2.026981e-17  0.19829545 0.004005448  0.1465276
#> Th__delta                3.206950e-01  0.63004541 0.004191793  0.1465276
#> neutrophil__otherimmune  6.532343e+00 -2.57014203 0.004569644  0.1465276
#> B__Th                    1.200000e+00  5.10701389 0.004805464  0.1465276
#> endothelial__stromal     1.988152e-01  0.60945565 0.006007150  0.1465276
#> B__delta                 3.560011e-01  0.65139563 0.006024443  0.1465276
#> alpha__otherimmune       9.231215e+00 -3.41798738 0.006421034  0.1465276
#> otherimmune__Th          8.673736e-04  0.01574398 0.006775857  0.1465276
#>                                from          to
#> naiveTc__otherimmune        naiveTc otherimmune
#> gamma__Th                     gamma          Th
#> macrophage__beta         macrophage        beta
#> Th__delta                        Th       delta
#> neutrophil__otherimmune  neutrophil otherimmune
#> B__Th                             B          Th
#> endothelial__stromal    endothelial     stromal
#> B__delta                          B       delta
#> alpha__otherimmune            alpha otherimmune
#> otherimmune__Th         otherimmune          Th

Again, we can present this spicy object as a bubble plot using the signifPlot() function by providing it with the spicy object.

signifPlot(
  spicyTestColPairs,
  marksToPlot = c(
    "alpha", "acinar", "ductal", "naiveTc", "neutrophil", "Tc",
    "Th", "otherimmune"
  )
)

5 sessionInfo()

sessionInfo()
#> R Under development (unstable) (2023-11-11 r85510)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 22.04.3 LTS
#> 
#> Matrix products: default
#> BLAS:   /home/biocbuild/bbs-3.19-bioc/R/lib/libRblas.so 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
#> 
#> 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       
#> 
#> time zone: America/New_York
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] imcRtools_1.9.0             SpatialExperiment_1.13.0   
#>  [3] SingleCellExperiment_1.25.0 SummarizedExperiment_1.33.0
#>  [5] Biobase_2.63.0              GenomicRanges_1.55.1       
#>  [7] GenomeInfoDb_1.39.1         IRanges_2.37.0             
#>  [9] S4Vectors_0.41.1            BiocGenerics_0.49.1        
#> [11] MatrixGenerics_1.15.0       matrixStats_1.1.0          
#> [13] ggplot2_3.4.4               spicyR_1.15.3              
#> [15] BiocStyle_2.31.0           
#> 
#> loaded via a namespace (and not attached):
#>   [1] splines_4.4.0               later_1.3.1                
#>   [3] bitops_1.0-7                svgPanZoom_0.3.4           
#>   [5] tibble_3.2.1                polyclip_1.10-6            
#>   [7] lifecycle_1.0.4             sf_1.0-14                  
#>   [9] rstatix_0.7.2               vroom_1.6.4                
#>  [11] lattice_0.22-5              MASS_7.3-60.1              
#>  [13] MultiAssayExperiment_1.29.0 backports_1.4.1            
#>  [15] magrittr_2.0.3              sass_0.4.7                 
#>  [17] rmarkdown_2.25              jquerylib_0.1.4            
#>  [19] yaml_2.3.7                  httpuv_1.6.12              
#>  [21] ClassifyR_3.7.1             sp_2.1-1                   
#>  [23] spatstat.sparse_3.0-3       DBI_1.1.3                  
#>  [25] minqa_1.2.6                 RColorBrewer_1.1-3         
#>  [27] abind_1.4-5                 zlibbioc_1.49.0            
#>  [29] purrr_1.0.2                 ggraph_2.1.0               
#>  [31] RCurl_1.98-1.13             tweenr_2.0.2               
#>  [33] GenomeInfoDbData_1.2.11     ggrepel_0.9.4              
#>  [35] RTriangle_1.6-0.12          spatstat.utils_3.0-4       
#>  [37] terra_1.7-55                pheatmap_1.0.12            
#>  [39] units_0.8-4                 goftest_1.2-3              
#>  [41] spatstat.random_3.2-1       DelayedMatrixStats_1.25.1  
#>  [43] svglite_2.1.2               codetools_0.2-19           
#>  [45] DelayedArray_0.29.0         scuttle_1.13.0             
#>  [47] DT_0.30                     ggforce_0.4.1              
#>  [49] tidyselect_1.2.0            raster_3.6-26              
#>  [51] farver_2.1.1                lme4_1.1-35.1              
#>  [53] viridis_0.6.4               spatstat.explore_3.2-5     
#>  [55] jsonlite_1.8.7              BiocNeighbors_1.21.0       
#>  [57] e1071_1.7-13                ellipsis_0.3.2             
#>  [59] tidygraph_1.2.3             survival_3.5-7             
#>  [61] systemfonts_1.0.5           tools_4.4.0                
#>  [63] Rcpp_1.0.11                 glue_1.6.2                 
#>  [65] gridExtra_2.3               SparseArray_1.3.1          
#>  [67] xfun_0.41                   mgcv_1.9-0                 
#>  [69] EBImage_4.45.0              dplyr_1.1.3                
#>  [71] HDF5Array_1.31.0            scam_1.2-14                
#>  [73] shinydashboard_0.7.2        withr_2.5.2                
#>  [75] numDeriv_2016.8-1.1         BiocManager_1.30.22        
#>  [77] fastmap_1.1.1               boot_1.3-28.1              
#>  [79] rhdf5filters_1.15.1         fansi_1.0.5                
#>  [81] digest_0.6.33               R6_2.5.1                   
#>  [83] mime_0.12                   colorspace_2.1-0           
#>  [85] tensor_1.5                  jpeg_0.1-10                
#>  [87] spatstat.data_3.0-3         utf8_1.2.4                 
#>  [89] tidyr_1.3.0                 generics_0.1.3             
#>  [91] data.table_1.14.8           class_7.3-22               
#>  [93] graphlayouts_1.0.2          htmlwidgets_1.6.2          
#>  [95] S4Arrays_1.3.0              pkgconfig_2.0.3            
#>  [97] gtable_0.3.4                XVector_0.43.0             
#>  [99] htmltools_0.5.7             carData_3.0-5              
#> [101] bookdown_0.36               fftwtools_0.9-11           
#> [103] scales_1.2.1                png_0.1-8                  
#> [105] knitr_1.45                  tzdb_0.4.0                 
#> [107] reshape2_1.4.4              rjson_0.2.21               
#> [109] nlme_3.1-163                nloptr_2.0.3               
#> [111] proxy_0.4-27                cachem_1.0.8               
#> [113] rhdf5_2.47.0                stringr_1.5.1              
#> [115] KernSmooth_2.23-22          parallel_4.4.0             
#> [117] vipor_0.4.5                 concaveman_1.1.0           
#> [119] pillar_1.9.0                grid_4.4.0                 
#> [121] vctrs_0.6.4                 promises_1.2.1             
#> [123] ggpubr_0.6.0                car_3.1-2                  
#> [125] beachmat_2.19.0             distances_0.1.9            
#> [127] xtable_1.8-4                beeswarm_0.4.0             
#> [129] evaluate_0.23               readr_2.1.4                
#> [131] magick_2.8.1                cli_3.6.1                  
#> [133] locfit_1.5-9.8              compiler_4.4.0             
#> [135] rlang_1.1.2                 crayon_1.5.2               
#> [137] ggsignif_0.6.4              labeling_0.4.3             
#> [139] classInt_0.4-10             plyr_1.8.9                 
#> [141] ggbeeswarm_0.7.2            stringi_1.8.1              
#> [143] viridisLite_0.4.2           deldir_1.0-9               
#> [145] BiocParallel_1.37.0         nnls_1.5                   
#> [147] cytomapper_1.15.0           lmerTest_3.1-3             
#> [149] munsell_0.5.0               tiff_0.1-11                
#> [151] spatstat.geom_3.2-7         Matrix_1.6-3               
#> [153] hms_1.1.3                   bit64_4.0.5                
#> [155] sparseMatrixStats_1.15.0    Rhdf5lib_1.25.0            
#> [157] shiny_1.7.5.1               highr_0.10                 
#> [159] igraph_1.5.1                broom_1.0.5                
#> [161] bslib_0.5.1                 bit_4.0.5