spicyR 1.17.3
if (!require("BiocManager")) {
install.packages("BiocManager")
}
BiocManager::install("spicyR")
# load required packages
library(spicyR)
library(ggplot2)
library(SingleCellExperiment)
library(SpatialExperiment)
library(imcRtools)
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.
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 SingleCellExperiment
object containing single-cell data of 160 images
from 8 subjects, with 20 images per subject.
data("diabetesData")
diabetesData
#> 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).
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.
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. However, given that there are 3 conditions, we can specify the desired
conditions by setting the order of our condition
factor. spicy()
will choose the first level of the factor as
the base condition and the second level as the comparison condition. spicy()
will also naturally coerce the condition
column into a
factor if not already a factor.
spicyTestPair <- spicy(
diabetesData,
condition = "stage",
subject = "case",
from = "beta",
to = "delta"
)
topPairs(spicyTestPair)
#> intercept coefficient p.value adj.pvalue from to
#> beta__delta 182.2128 -34.00085 0.008724782 0.008724782 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.
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"
)
topPairs(spicyTest)
#> intercept coefficient p.value adj.pvalue
#> beta__delta 1.815457e+02 -48.736479 0.0005032491 0.07169646
#> delta__beta 1.817943e+02 -48.166076 0.0005601286 0.07169646
#> B__unknown 0.000000e+00 11.770769 0.0052340673 0.42052603
#> delta__delta 2.089550e+02 -52.061196 0.0125129204 0.42052603
#> unknown__macrophage 1.007336e+01 -15.826911 0.0207411810 0.42052603
#> unknown__B 1.474842e-15 12.142843 0.0225856401 0.42052603
#> macrophage__unknown 1.004422e+01 -14.471640 0.0244668217 0.42052603
#> B__Th -1.571468e-15 26.847613 0.0245045799 0.42052603
#> otherimmune__naiveTc -9.292507e+00 33.584768 0.0255811877 0.42052603
#> ductal__ductal 1.481580e+01 -8.632552 0.0266937176 0.42052603
#> 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
We can also examine the L-function metrics of individual images by using the
convenient bind()
function on our spicyTest results object.
bind(spicyTest)[1:5, 1:5]
#> imageID condition subject ductal__ductal stromal__ductal
#> 1 A09 Onset 6362 3.5882456 2.221666
#> 2 A11 Onset 6362 3.5803460 -5.070033
#> 3 A16 Onset 6362 29.8249524 7.859502
#> 4 A17 Onset 6362 30.8885146 20.136107
#> 5 A25 Onset 6362 0.9820586 -3.173230
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)
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,
colData = colData(diabetesData))
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
#> A09 2.590909
#> A11 4.802817
#> A16 1.769231
#> A17 0.500000
#> A25 0.000000
#> A39 1.354839
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
)
topPairs(spicyTestColPairs)
#> intercept coefficient p.value adj.pvalue
#> otherimmune__acinar 8.013025e+00 -3.92656804 0.002070014 0.1470671
#> macrophage__macrophage 1.524442e-03 0.01186938 0.003008538 0.1470671
#> acinar__delta -2.026981e-17 0.19829545 0.004005448 0.1470671
#> gamma__endothelial 3.206950e-01 0.63099515 0.004186254 0.1470671
#> ductal__stromal 6.531116e+00 -2.57189740 0.004531110 0.1470671
#> ductal__otherimmune 1.200000e+00 5.09868056 0.004819382 0.1470671
#> Th__otherimmune 1.984338e-01 0.61055124 0.005893992 0.1470671
#> unknown__alpha 3.560011e-01 0.65139563 0.006024443 0.1470671
#> Th__stromal 9.233590e+00 -3.41854439 0.006371488 0.1470671
#> endothelial__ductal 8.673736e-04 0.01574398 0.006775857 0.1470671
#> from to
#> otherimmune__acinar otherimmune acinar
#> macrophage__macrophage macrophage macrophage
#> acinar__delta acinar delta
#> gamma__endothelial gamma endothelial
#> ductal__stromal ductal stromal
#> ductal__otherimmune ductal otherimmune
#> Th__otherimmune Th otherimmune
#> unknown__alpha unknown alpha
#> Th__stromal Th stromal
#> endothelial__ductal endothelial ductal
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"
)
)
sessionInfo()
#> R version 4.4.1 (2024-06-14)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 22.04.4 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.20-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.11.0 SpatialExperiment_1.15.1
#> [3] SingleCellExperiment_1.27.2 SummarizedExperiment_1.35.1
#> [5] Biobase_2.65.0 GenomicRanges_1.57.1
#> [7] GenomeInfoDb_1.41.1 IRanges_2.39.2
#> [9] S4Vectors_0.43.2 BiocGenerics_0.51.0
#> [11] MatrixGenerics_1.17.0 matrixStats_1.3.0
#> [13] ggplot2_3.5.1 spicyR_1.17.3
#> [15] BiocStyle_2.33.1
#>
#> loaded via a namespace (and not attached):
#> [1] splines_4.4.1 later_1.3.2
#> [3] bitops_1.0-8 svgPanZoom_0.3.4
#> [5] tibble_3.2.1 polyclip_1.10-7
#> [7] lifecycle_1.0.4 sf_1.0-16
#> [9] rstatix_0.7.2 vroom_1.6.5
#> [11] lattice_0.22-6 MASS_7.3-61
#> [13] MultiAssayExperiment_1.31.4 backports_1.5.0
#> [15] magrittr_2.0.3 sass_0.4.9
#> [17] rmarkdown_2.27 jquerylib_0.1.4
#> [19] yaml_2.3.10 httpuv_1.6.15
#> [21] ClassifyR_3.9.3 sp_2.1-4
#> [23] spatstat.sparse_3.1-0 DBI_1.2.3
#> [25] minqa_1.2.7 RColorBrewer_1.1-3
#> [27] abind_1.4-5 zlibbioc_1.51.1
#> [29] purrr_1.0.2 ggraph_2.2.1
#> [31] RCurl_1.98-1.16 tweenr_2.0.3
#> [33] GenomeInfoDbData_1.2.12 ggrepel_0.9.5
#> [35] RTriangle_1.6-0.13 spatstat.utils_3.0-5
#> [37] terra_1.7-78 pheatmap_1.0.12
#> [39] units_0.8-5 goftest_1.2-3
#> [41] spatstat.random_3.3-1 svglite_2.1.3
#> [43] codetools_0.2-20 DelayedArray_0.31.11
#> [45] scuttle_1.15.2 DT_0.33
#> [47] ggforce_0.4.2 tidyselect_1.2.1
#> [49] raster_3.6-26 UCSC.utils_1.1.0
#> [51] farver_2.1.2 lme4_1.1-35.5
#> [53] viridis_0.6.5 spatstat.explore_3.3-1
#> [55] jsonlite_1.8.8 BiocNeighbors_1.23.0
#> [57] e1071_1.7-14 tidygraph_1.3.1
#> [59] survival_3.7-0 systemfonts_1.1.0
#> [61] tools_4.4.1 Rcpp_1.0.13
#> [63] glue_1.7.0 gridExtra_2.3
#> [65] SparseArray_1.5.31 xfun_0.46
#> [67] mgcv_1.9-1 EBImage_4.47.0
#> [69] dplyr_1.1.4 HDF5Array_1.33.5
#> [71] scam_1.2-17 shinydashboard_0.7.2
#> [73] withr_3.0.1 numDeriv_2016.8-1.1
#> [75] BiocManager_1.30.23 fastmap_1.2.0
#> [77] boot_1.3-30 rhdf5filters_1.17.0
#> [79] fansi_1.0.6 digest_0.6.36
#> [81] R6_2.5.1 mime_0.12
#> [83] colorspace_2.1-1 tensor_1.5
#> [85] jpeg_0.1-10 spatstat.data_3.1-2
#> [87] utf8_1.2.4 tidyr_1.3.1
#> [89] generics_0.1.3 data.table_1.15.4
#> [91] class_7.3-22 graphlayouts_1.1.1
#> [93] httr_1.4.7 htmlwidgets_1.6.4
#> [95] S4Arrays_1.5.7 pkgconfig_2.0.3
#> [97] gtable_0.3.5 XVector_0.45.0
#> [99] htmltools_0.5.8.1 carData_3.0-5
#> [101] bookdown_0.40 fftwtools_0.9-11
#> [103] scales_1.3.0 ggupset_0.4.0
#> [105] png_0.1-8 spatstat.univar_3.0-0
#> [107] knitr_1.48 tzdb_0.4.0
#> [109] reshape2_1.4.4 rjson_0.2.21
#> [111] nlme_3.1-165 nloptr_2.1.1
#> [113] proxy_0.4-27 cachem_1.1.0
#> [115] rhdf5_2.49.0 stringr_1.5.1
#> [117] KernSmooth_2.23-24 parallel_4.4.1
#> [119] vipor_0.4.7 concaveman_1.1.0
#> [121] pillar_1.9.0 grid_4.4.1
#> [123] vctrs_0.6.5 promises_1.3.0
#> [125] ggpubr_0.6.0 car_3.1-2
#> [127] beachmat_2.21.5 distances_0.1.11
#> [129] xtable_1.8-4 beeswarm_0.4.0
#> [131] evaluate_0.24.0 tinytex_0.52
#> [133] readr_2.1.5 magick_2.8.4
#> [135] cli_3.6.3 locfit_1.5-9.10
#> [137] compiler_4.4.1 rlang_1.1.4
#> [139] crayon_1.5.3 ggsignif_0.6.4
#> [141] labeling_0.4.3 classInt_0.4-10
#> [143] plyr_1.8.9 ggbeeswarm_0.7.2
#> [145] stringi_1.8.4 viridisLite_0.4.2
#> [147] deldir_2.0-4 BiocParallel_1.39.0
#> [149] nnls_1.5 cytomapper_1.17.0
#> [151] lmerTest_3.1-3 munsell_0.5.1
#> [153] tiff_0.1-12 spatstat.geom_3.3-2
#> [155] Matrix_1.7-0 hms_1.1.3
#> [157] bit64_4.0.5 Rhdf5lib_1.27.0
#> [159] shiny_1.9.1 highr_0.11
#> [161] igraph_2.0.3 broom_1.0.6
#> [163] memoise_2.0.1 bslib_0.8.0
#> [165] bit_4.0.5