Spatial Mixed-Effects Modelling with spicy
30 August 2023
Package
spicyR 1.12.2
Contents
# load required packages
library(spicyR)
library(ggplot2)
library(SingleCellExperiment)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_SCE is a SingleCellExperiment object containing single-cell data of 160 images
from 8 subjects, with 20 images per subjects.
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 deltaWe 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 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 ductalAgain, 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 signifPlot() function by
providing it the spicy object obtained.
signifPlot(
spicyTest,
breaks = c(-3, 3, 1),
marksToPlot = c(
"alpha", "beta", "gamma", "delta",
"B", "naiveTc", "Th", "Tc", "neutrophil", "macrophage"
)
)
5 sessionInfo()
sessionInfo()
#> R version 4.3.1 (2023-06-16)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 22.04.3 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.17-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] SingleCellExperiment_1.22.0 SummarizedExperiment_1.30.2
#> [3] Biobase_2.60.0 GenomicRanges_1.52.0
#> [5] GenomeInfoDb_1.36.2 IRanges_2.34.1
#> [7] S4Vectors_0.38.1 BiocGenerics_0.46.0
#> [9] MatrixGenerics_1.12.3 matrixStats_1.0.0
#> [11] ggplot2_3.4.3 spicyR_1.12.2
#> [13] BiocStyle_2.28.0
#>
#> loaded via a namespace (and not attached):
#> [1] RColorBrewer_1.1-3 jsonlite_1.8.7
#> [3] MultiAssayExperiment_1.26.0 magrittr_2.0.3
#> [5] spatstat.utils_3.0-3 magick_2.7.5
#> [7] farver_2.1.1 nloptr_2.0.3
#> [9] rmarkdown_2.24 zlibbioc_1.46.0
#> [11] vctrs_0.6.3 minqa_1.2.5
#> [13] spatstat.explore_3.2-1 DelayedMatrixStats_1.22.6
#> [15] RCurl_1.98-1.12 rstatix_0.7.2
#> [17] htmltools_0.5.6 S4Arrays_1.0.6
#> [19] broom_1.0.5 Rhdf5lib_1.22.0
#> [21] rhdf5_2.44.0 sass_0.4.7
#> [23] bslib_0.5.1 plyr_1.8.8
#> [25] cachem_1.0.8 lifecycle_1.0.3
#> [27] pkgconfig_2.0.3 Matrix_1.6-1
#> [29] R6_2.5.1 fastmap_1.1.1
#> [31] GenomeInfoDbData_1.2.10 digest_0.6.33
#> [33] numDeriv_2016.8-1.1 colorspace_2.1-0
#> [35] tensor_1.5 dqrng_0.3.0
#> [37] ggpubr_0.6.0 beachmat_2.16.0
#> [39] labeling_0.4.3 fansi_1.0.4
#> [41] spatstat.sparse_3.0-2 polyclip_1.10-4
#> [43] abind_1.4-5 mgcv_1.9-0
#> [45] compiler_4.3.1 withr_2.5.0
#> [47] backports_1.4.1 BiocParallel_1.34.2
#> [49] carData_3.0-5 highr_0.10
#> [51] HDF5Array_1.28.1 ggforce_0.4.1
#> [53] R.utils_2.12.2 ggsignif_0.6.4
#> [55] MASS_7.3-60 concaveman_1.1.0
#> [57] DelayedArray_0.26.7 rjson_0.2.21
#> [59] tools_4.3.1 goftest_1.2-3
#> [61] R.oo_1.25.0 glue_1.6.2
#> [63] nlme_3.1-163 rhdf5filters_1.12.1
#> [65] grid_4.3.1 ClassifyR_3.4.9
#> [67] reshape2_1.4.4 generics_0.1.3
#> [69] gtable_0.3.4 spatstat.data_3.0-1
#> [71] R.methodsS3_1.8.2 tidyr_1.3.0
#> [73] data.table_1.14.8 car_3.1-2
#> [75] utf8_1.2.3 XVector_0.40.0
#> [77] spatstat.geom_3.2-4 pillar_1.9.0
#> [79] stringr_1.5.0 limma_3.56.2
#> [81] splines_4.3.1 dplyr_1.1.2
#> [83] tweenr_2.0.2 lattice_0.21-8
#> [85] survival_3.5-7 deldir_1.0-9
#> [87] tidyselect_1.2.0 locfit_1.5-9.8
#> [89] scuttle_1.10.2 knitr_1.43
#> [91] bookdown_0.35 edgeR_3.42.4
#> [93] xfun_0.40 DropletUtils_1.20.0
#> [95] pheatmap_1.0.12 scam_1.2-14
#> [97] stringi_1.7.12 yaml_2.3.7
#> [99] boot_1.3-28.1 evaluate_0.21
#> [101] codetools_0.2-19 tibble_3.2.1
#> [103] BiocManager_1.30.22 cli_3.6.1
#> [105] munsell_0.5.0 jquerylib_0.1.4
#> [107] Rcpp_1.0.11 spatstat.random_3.1-5
#> [109] parallel_4.3.1 sparseMatrixStats_1.12.2
#> [111] bitops_1.0-7 lme4_1.1-34
#> [113] SpatialExperiment_1.10.0 lmerTest_3.1-3
#> [115] scales_1.2.1 purrr_1.0.2
#> [117] crayon_1.5.2 rlang_1.1.1