---
title: "Reconstruction and analysis of pancreatic islets from IMC data"
author: 
  - name: "Samuel Gunz"
    affiliation:
      - &DMLS Department of Molecular Life Sciences, University of Zurich, Switzerland
      - &SIB SIB Swiss Institute of Bioinformatics, University of Zurich, Switzerland
    email: "samuel.gunz@uzh.ch"
  - name: Mark D. Robinson
    affiliation:
      - *DMLS
      - *SIB
package: "`r BiocStyle::Biocpkg('sosta')`"
output:
  BiocStyle::html_document
abstract: >
  Reconstruction and analysis of pancreatic islets from IMC data using the `sosta` package
vignette: >
  %\VignetteIndexEntry{Reconstruction and analysis of pancreatic islets from IMC data}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: sosta.bib
---

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

# Installation

`r BiocStyle::Biocpkg('sosta')` can be loaded from Bioconductor and installed as follows

```{r installation, include = TRUE, eval=FALSE}
if (!requireNamespace("BiocManager")) {
    install.packages("BiocManager")
}
BiocManager::install("sosta")
```

# Setup

```{r setup, message=FALSE}
library("sosta")
library("dplyr")
library("tidyr")
library("ggplot2")
library("ggspavis")
library("imcdatasets")
```

# Introduction

In this vignette we will use data from the package `r BiocStyle::Biocpkg('imcdatasets')`. The dataset contains imaging mass cytometry (IMC) data of pancreatic islets of human donors at different stages of type 1 diabetes (T1D) and healthy controls [@damondMapHumanType2019]. Note that we will only use a subset of the images `full_dataset = FALSE`.

```{r loading, echo=FALSE, message=FALSE}
# load the data
spe <- imcdatasets::Damond_2019_Pancreas("spe", full_dataset = FALSE)
``` 

First we plot the data for illustration. As we have multiple images per patient we will subset one patient and a few slides.

```{r}
plotSpots(
    spe[, spe[["patient_id"]] == 6126 &
        spe[["image_name"]] %in% c("E02", "E03", "E04")],
    annotate = "cell_category",
    sample_id = "image_number",
    in_tissue = NULL,
    y_reverse = FALSE
) + facet_wrap(~image_name)
```

The goal is to reconstruct / segment and quantify the pancreatic islets.

# Reconstruction of pancreatic islets

## Reconstruction of pancreatic islets for one image

In this example we will reconstruct the islets based on the point pattern density of the islet cells. We will start with estimating the parameters that we use for reconstruction afterwards. For one image this can be illustrated as follows.

```{r}
shapeIntensityImage(
    spe,
    marks = "cell_category",
    image_col = "image_name",
    image_id = "E04",
    mark_select = "islet"
)
```

We see the density (pixel) image on the left and a histogram of the intensity values on the right. The estimated threshold is roughly the mean between the two modes of the (truncated) pixel intensity distribution.

This was done for one image. The function `estimateReconstructionParametersSPE` returns two plots with the estimated parameters for each image. The parameter `bndw` is the bandwidth parameter that is used for estimating the intensity profile of the point pattern. The parameter `thres` is the estimated parameter for the density threshold for reconstruction. We subset 20 images to speed up computation.

```{r}
n <- estimateReconstructionParametersSPE(
    spe,
    marks = "cell_category",
    image_col = "image_name",
    mark_select = "islet",
    nimages = 20,
    plot_hist = TRUE
)
```
We will use the mean of the two estimated vectors as our parameters.

```{r}
(thres_spe <- mean(n$thres))
(bndw_spe <- mean(n$bndw))
```

We cam now use the function `reconstructShapeDensity` to segment the islet of one image. The result is a `r BiocStyle::CRANpkg('sf')` polygon [@pebesmaSimpleFeaturesStandardized2018]. 

```{r}
islet <- reconstructShapeDensityImage(
    spe,
    marks = "cell_category",
    image_col = "image_name",
    image_id = "E04",
    mark_select = "islet",
    bndw = bndw_spe,
    dim = 500,
    thres = thres_spe
)
```

We can plot both the points and the estimated islets polygon.

```{r}
plotSpots(
    spe[, spe[["image_name"]] %in% c("E04")],
    annotate = "cell_category",
    sample_id = "image_number",
    in_tissue = NULL,
    y_reverse = FALSE,
) +
    geom_sf(
        data = islet,
        fill = NA,
        color = "darkblue",
        inherit.aes = FALSE, # this is important
        linewidth = 0.75
    )
```

If no arguments are given, the function `reconstructShapeDensityImage` estimates the reconstruction parameters internally.

```{r}
islet_2 <- reconstructShapeDensityImage(
    spe,
    marks = "cell_category",
    image_col = "image_name",
    image_id = "E04",
    mark_select = "islet",
    dim = 500
)
```


```{r}
plotSpots(
    spe[, spe[["image_name"]] %in% c("E04")],
    annotate = "cell_category",
    sample_id = "image_number",
    in_tissue = NULL,
    y_reverse = FALSE,
) +
    geom_sf(
        data = islet_2,
        fill = NA,
        color = "darkblue",
        inherit.aes = FALSE,
        linewidth = 0.75
    )
```

## Reconstruction of pancreatic islets for all images

The function `reconstructShapeDensitySPE` reconstructs the islets for all images in the `spe` object. We use the estimated parameters from above.

```{r, eval=TRUE}
all_islets <- reconstructShapeDensitySPE(
    spe,
    marks = "cell_category",
    image_col = "image_name",
    mark_select = "islet",
    bndw = bndw_spe,
    thres = thres_spe,
    ncores = 2
)
```

# Calculation of structure metrics

As we have islets for all images, we now use the function `totalShapeMetrics` to calculate a set of metrics related to the shape of the islets.

```{r}
islet_shape_metrics <- totalShapeMetrics(all_islets)
```

The result is a simple feature collection with polygons. We will add some patient level information to the simple feature collection for plotting afterwards.

```{r}
patient_data <- colData(spe) |>
    as_tibble() |>
    group_by(image_name) |>
    select(all_of(
        c(
            "patient_stage",
            "tissue_slide",
            "tissue_region",
            "patient_id",
            "patient_disease_duration",
            "patient_age",
            "patient_gender",
            "patient_ethnicity",
            "patient_BMI",
            "sample_id"
        )
    )) |>
    unique()
```


```{r}
all_islets <- dplyr::left_join(all_islets, patient_data, by = "image_name")
all_islets <- cbind(all_islets, t(islet_shape_metrics))
```

## Plot structure metrics

We use PCA to get an overview of the different features. Each dot represent one structure.

```{r}
library(ggfortify)

autoplot(
    prcomp(t(islet_shape_metrics), scale. = TRUE),
    x = 1,
    y = 2,
    data = all_islets,
    color = "patient_stage",
    size = 2,
    # shape = 'type',
    loadings = TRUE,
    loadings.colour = "steelblue3",
    loadings.label = TRUE,
    loadings.label.size = 3,
    loadings.label.repel = TRUE,
    loadings.label.colour = "black"
) +
    scale_color_brewer(palette = "Dark2") +
    theme_bw() +
    coord_fixed()
```

We can use boxplots to investigate differences of shape metrics between patient stages.

```{r, fig.width=8, fig.height=6}
all_islets |>
    sf::st_drop_geometry() |>
    select(patient_stage, rownames(islet_shape_metrics)) |>
    pivot_longer(-patient_stage) |>
    ggplot(aes(x = patient_stage, y = value, fill = patient_stage)) +
    geom_boxplot() +
    facet_wrap(~name, scales = "free") +
    scale_fill_brewer(palette = "Dark2") +
    scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
    guides(fill = "none") +
    theme_bw()
```

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

# References