---
title: "01 -- Overview of sosta"
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: >
  In this vignette we show different functionalities of `sosta`. First we perform a density-based reconstruction of structures based on spatial coordinates. Next, we will calcuate different metrics on the structure and cell level.
vignette: >
  %\VignetteIndexEntry{Overview of sosta}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
bibliography: sosta.bib
editor_options: 
  chunk_output_type: console
---

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

# Installation

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

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

# Setup

For this vignette, we will need several additional packages:

```{r setup, message=FALSE, warning=FALSE}
library("sosta")
library("dplyr")
library("tidyr")
library("ggplot2")
library("sf")
library("SpatialExperiment")

theme_set(theme_bw())
```

# Introduction

As an example, we will load an simulated dataset with three images and three cell types A, B and C, stored as a `r BiocStyle::Biocpkg('SpatialExperiment')` object:

```{r loading, echo=TRUE, message=FALSE}
# load the data
data("sostaSPE")
sostaSPE
``` 

```{r}
cbind(colData(sostaSPE), spatialCoords(sostaSPE)) |>
    as.data.frame() |>
    ggplot(aes(x = x, y = y, color = cellType)) +
    geom_point(size = 0.25) +
    facet_wrap(~imageName) +
    coord_equal()
```

The goal is to reconstruct, or segment, the cellular "structures" given by cell type A.

# Structure reconstruction

## Reconstruction of structures for one image

In this example, we will reconstruct cellular *structures* based on the point pattern density of the cells of type A. We will start with estimating parameters that are used for reconstruction afterwards. For one image, this can be illustrated as follows:

```{r}
shapeIntensityImage(
    sostaSPE,
    marks = "cellType",
    imageCol = "imageName",
    imageId = "image1",
    markSelect = "A"
)
```

We see the pixel-wise density 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. The dimensions of the pixel image are specified on the bottom left. The dimensions correspond to the density image, setting a higher value for the `dim` parameter will result in a higher resolution density image but will not change how many points are within a pixel. This can be changed by varying the smoothing parameter (`bndw`).

This was done for one image. The function `estimateReconstructionParametersSPE` returns two plots with the estimated parameters for each image in the dataset. 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.

```{r}
n <- estimateReconstructionParametersSPE(
    sostaSPE,
    marks = "cellType",
    imageCol = "imageName",
    markSelect = "A",
    plotHist = TRUE
)
```
We will use the mean of the two estimated vectors as our parameters.

```{r}
(thresSPE <- mean(n$thres))
(bndwSPE <- mean(n$bndw))
```

We now use the function `reconstructShapeDensity` to segment the cell-type-A structure into regions. The result is a collection of `r BiocStyle::CRANpkg('sf')` polygons [@pebesmaSimpleFeaturesStandardized2018]. 

```{r}
(struct <- reconstructShapeDensityImage(
    sostaSPE,
    marks = "cellType",
    imageCol = "imageName",
    imageId = "image1",
    markSelect = "A",
    bndw = bndwSPE,
    dim = 500,
    thres = thresSPE
))
```

Let's plot both the points and the segmented polygons.

```{r}
cbind(
    colData(sostaSPE[, sostaSPE$imageName == "image1"]),
    spatialCoords(sostaSPE[, sostaSPE$imageName == "image1"])
) |>
    as.data.frame() |>
    ggplot(aes(x = x, y = y, color = cellType)) +
    geom_point(size = 0.5) +
    facet_wrap(~imageName) +
    coord_equal() +
    geom_sf(
        data = struct,
        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. The bandwidth is estimated using cross validation using the function `bw.diggle` of the `r BiocStyle::CRANpkg('spatstat.explore')` package. The threshold is estimated by taking the mean between the two modes of the pixel intensity distribution as illustrated above.

```{r}
struct2 <- reconstructShapeDensityImage(
    sostaSPE,
    marks = "cellType",
    imageCol = "imageName",
    imageId = "image1",
    markSelect = "A",
    dim = 500
)
```


```{r}
cbind(
    colData(sostaSPE[, sostaSPE$imageName == "image1"]),
    spatialCoords(sostaSPE[, sostaSPE$imageName == "image1"])
) |>
    as.data.frame() |>
    ggplot(aes(x = x, y = y, color = cellType)) +
    geom_point(size = 0.5) +
    facet_wrap(~imageName) +
    coord_equal() +
    geom_sf(
        data = struct2,
        fill = NA,
        color = "darkblue",
        inherit.aes = FALSE, # this is important
        linewidth = 0.75
    )
```

## Reconstruction of structures for all images

The function `reconstructShapeDensitySPE` reconstructs the cell-type-A structure for all images in the `spe` object. We use the estimated parameters from above.

```{r, eval=TRUE}
allStructs <- reconstructShapeDensitySPE(
    sostaSPE,
    marks = "cellType",
    imageCol = "imageName",
    markSelect = "A",
    bndw = bndwSPE,
    thres = thresSPE,
    nCores = 1
)
```

## Intersection with cells

Next, we assign cells in the `SpatialExperiment` object to their specific structure. This is useful for later calculations. Such that the function can match the cells to the correct structures, we need the `colnames` to be defined in the `SpatialExperiment` object.

```{r}
# Define colnames by numbering the cells
colnames(sostaSPE) <- paste0("cell_", c(1:dim(sostaSPE)[2]))
```


```{r}
assign <- assingCellsToStructures(sostaSPE, allStructs,
    imageCol = "imageName", nCores = 1
)
# Assign using the correct order of the columns in the spe object
sostaSPE$structAssign <- assign[colnames(sostaSPE)]
```


```{r}
cbind(colData(sostaSPE), spatialCoords(sostaSPE)) |>
    as.data.frame() |>
    ggplot(aes(x = x, y = y, color = structAssign)) +
    geom_point(size = 0.25) +
    facet_wrap(~imageName) +
    coord_equal()
```

# Structure level metrics

## Proportion of cell types within structures

Using the function `cellTypeProportions`, we can estimate the proportion of cell types within each individual structure.

```{r}
cellTypeProportions(sostaSPE, "structAssign", "cellType")
```

## Shape Metrics

The function `totalShapeMetrics` calculates a set of geometric metrics related to the shape of the structures.

```{r}
shapeMetrics <- totalShapeMetrics(allStructs)
head(shapeMetrics)
```


```{r}
cbind(allStructs, t(shapeMetrics)) |>
    ggplot() +
    geom_sf(aes(fill = Area)) +
    facet_wrap(~imageName)
```

# Cell level metrics

## Distance to structure border

Using the function `minBoundaryDistances`, we can compute the distance to the border structure. Negative values indicate that the points lie inside the structure.

```{r}
sostaSPE$minDist <- minBoundaryDistances(
    spe = sostaSPE, imageCol = "imageName",
    structColumn = "structAssign", allStructs = allStructs
)

cbind(colData(sostaSPE), spatialCoords(sostaSPE)) |>
    as.data.frame() |>
    ggplot(aes(x = x, y = y, color = minDist)) +
    geom_point(size = 0.25) +
    facet_wrap(~imageName) +
    coord_equal() +
    scale_colour_gradient2() +
    geom_sf(
        data = allStructs,
        fill = NA,
        inherit.aes = FALSE
    ) +
    facet_wrap(~imageName)
```

This information can be used to define border cells by thresholding to a range of positive and negative values.

```{r}
sostaSPE$border <- ifelse(abs(sostaSPE$minDist) < 3, TRUE, FALSE)


cbind(colData(sostaSPE), spatialCoords(sostaSPE)) |>
    as.data.frame() |>
    ggplot(aes(x = x, y = y, color = border)) +
    geom_point(size = 0.25) +
    facet_wrap(~imageName) +
    coord_equal() +
    geom_sf(
        data = allStructs,
        fill = NA,
        inherit.aes = FALSE
    ) +
    facet_wrap(~imageName)
```

Alternatively, borders can be defined using `st_difference` and `st_buffer` on the structure object. In this case we will have an `sf` polygon that correspond to the border region.

```{r}
borders <- lapply(
    st_geometry(allStructs),
    \(x) st_difference(st_buffer(x, 3), st_buffer(x, -3))
) |>
    st_as_sfc() |>
    st_as_sf() # both functions needed for proper conversion

borders$imageName <- allStructs$imageName
borders$borderID <- paste0("border_", allStructs$structID)

borderAssign <- assingCellsToStructures(sostaSPE,
    borders,
    imageCol = "imageName",
    uniqueId = "borderID",
    nCores = 1
)
# Assign using the correct order of the columns in the spe object
sostaSPE$borderSf <- borderAssign[colnames(sostaSPE)]
```


```{r}
cbind(colData(sostaSPE), spatialCoords(sostaSPE)) |>
    as.data.frame() |>
    ggplot(aes(x = x, y = y, color = borderSf)) +
    geom_point(size = 0.25) +
    facet_wrap(~imageName) +
    coord_equal() +
    geom_sf(
        data = borders,
        fill = NA,
        inherit.aes = FALSE
    ) +
    facet_wrap(~imageName)
```

## Structure boundary vs FOV boundary

When defining border cells, the field of view (FOV) can cut through the reconstructed polygon(s). This can cause cells that are actually inside the structure to be classified as border cells. To avoid this, we compute the distance to the FOV boundary and remove cells that are close to it.

The FOV in this example ranges from coordinates 0 to 128 in the $x$ and $y$ direction. We therefore define polygons that span the entire FOV. Since there are multiple samples, we replicate these polygons and label them accordingly.

```{r}
# create the fov bounding box using st_bbox
fov_bbox <- st_bbox(c(xmin = 0, ymin = 0, xmax = 128, ymax = 128))

# convert to polygon
fov <- st_as_sfc(fov_bbox)

# one fov per image
fovBorder <- rep(fov, length(unique(sostaSPE$imageName))) |>
  st_as_sf()

# name the polygons accordingly
fovBorder$imageName <- unique(sostaSPE$imageName)

fovBorder
```

Now we use the function `minBoundaryDistances` to get the distances to the FOV border. 

```{r}
sostaSPE$fovBorderDist <- minBoundaryDistances(
  spe = sostaSPE,
  imageCol = "imageName",
  structColumn = "imageName", # here any variable that is true for all cells
  allStructs = fovBorder
) |> 
  abs()
```

```{r}
cbind(colData(sostaSPE), spatialCoords(sostaSPE)) |>
    as.data.frame() |>
    ggplot(aes(x = x, y = y, color = fovBorderDist)) +
    geom_point(size = 0.25) +
    facet_wrap(~imageName) +
    coord_equal() +
    facet_wrap(~imageName)
```

We can use the distance to the FOV border to correct the border assignments.

```{r}
sostaSPE$borderCorrected <- 
  ifelse(sostaSPE$fovBorderDist > 5, sostaSPE$borderSf, NA)
```

```{r}
cbind(colData(sostaSPE), spatialCoords(sostaSPE)) |>
    as.data.frame() |>
    ggplot(aes(x = x, y = y, color = borderCorrected)) +
    geom_point(size = 0.25) +
    facet_wrap(~imageName) +
    coord_equal() +
    facet_wrap(~imageName)
```

# Session Info

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

# References
