knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) library(BiocStyle)

# load required packages library(spicyR) library(ggplot2)

if (!require("BiocManager")) install.packages("BiocManager") BiocManager::install("spicyR")

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 `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 cellSummary(diabetesData) imagePheno(diabetesData)

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

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 topPairs(spicyTestPair)

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.

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
topPairs(spicyTest)
```

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

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
topPairs(spicyTestBootstrap)
```

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

```
sessionInfo()
```

**Any scripts or data that you put into this service are public.**

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.