---
title: "Spatial Mixed-Effects Modelling with spicy"
date: "`r BiocStyle::doc_date()`"
params:
  test: FALSE
author:
- name: Nicolas Canete
  affiliation:  
  - &WIMR Westmead Institute for Medical Research, University of Sydney, Australia
  email: nicolas.canete@sydney.edu.au
- name: Ellis Patrick
  affiliation:
  - &WIMR Westmead Institute for Medical Research, University of Sydney, Australia
  - School of Mathematics and Statistics, University of Sydney, Australia
package: "`r BiocStyle::pkg_ver('spicyR')`"
vignette: >
  %\VignetteIndexEntry{"Introduction to spicy"}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
output: 
  BiocStyle::html_document
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  warning = FALSE
)
library(BiocStyle)
```


# Installation

```{r, eval = FALSE}
if (!require("BiocManager")) {
  install.packages("BiocManager")
}
BiocManager::install("spicyR")
```

```{r warning=FALSE, message=FALSE}
# load required packages
library(spicyR)
library(ggplot2)
library(SingleCellExperiment)
library(SpatialExperiment)
library(imcRtools)
```

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

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


```{r}
data("diabetesData_SCE")
diabetesData_SCE
```

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

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

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

```{r message=FALSE}
spicyTestPair <- spicy(
  diabetesData_SCE,
  condition = "stage",
  subject = "case",
  from = "beta",
  to = "delta"
)

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 delta cells cells are significantly less likely to be found near beta cells
in the onset diabetes group compared to the non-diabetic control.

## 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()`.

```{r message=FALSE}
spicyTest <- spicy(
  diabetesData_SCE,
  condition = "stage",
  subject = "case"
)

topPairs(spicyTest)
```
We can also examine the L-function metrics of individual images by using the
convenient `bind()` function on our spicyTest results object.
```{r}
bind(spicyTest)
```


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.
```{r}
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`.

```{r}
spicyBoxPlot(results = spicyTest,
             # from = "beta",
             # to = "delta",
             rank = 1)
```

## 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. 
```{r}
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"))

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

head(pairAbundances["delta__delta"])
```

`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.
```{r}
spicyTestColPairs <- spicy(
  diabetesData_SPE,
  condition = "stage",
  subject = "case",
  alternateResult = pairAbundances,
  weights = FALSE
)

topPairs(spicyTestColPairs)
```

Again, we can present this `spicy` object as a bubble plot using the 
`signifPlot()` function by providing it with the `spicy` object.
```{r}
signifPlot(
  spicyTestColPairs,
  marksToPlot = c(
    "alpha", "acinar", "ductal", "naiveTc", "neutrophil", "Tc",
    "Th", "otherimmune"
  )
)
```

# sessionInfo()

```{r}
sessionInfo()
```