---

title: "Cell/Spot Analysis"

output: 

  html_document:

    toc: true

    toc_float:

      collapsed: false

      smooth_scroll: false

---

<style>

.title{

  display: none;

}

body {

  text-align: justify

}

.center {

  display: block;

  margin-left: auto;

  margin-right: auto;

}

</style>

```{css, echo=FALSE}

.watch-out {

  color: black;

}

```

```{r setup, include=FALSE}

# use rmarkdown::render_site(envir = knitr::knit_global())

knitr::opts_chunk$set(highlight = TRUE, echo = TRUE)

```

<br>

# Xenium Data Analysis

VoltRon is an end-to-end spatial omic analysis package which also supports investigating spatial points in single cell resolution. VoltRon includes essential built-in functions capable of **filtering**, **processing** and **clustering** as well as **visualizing** spatial datasets with a goal of cell type discovery and annotation. 

In this use case, we analyse readouts of the experiments conducted on example tissue sections analysed by the [Xenium In Situ](https://www.10xgenomics.com/platforms/xenium) platform. Two tissue sections of 5 $\mu$m tickness are derived from a single formalin-fixed, paraffin-embedded (FFPE) breast cancer tissue block. More information on the spatial datasets and the study can be also be found on the [BioArxiv preprint](https://www.biorxiv.org/content/10.1101/2022.10.06.510405v1).

You can import these readouts from the [10x Genomics website](https://www.10xgenomics.com/products/xenium-in-situ/preview-dataset-human-breast) (specifically, import **In Situ Replicate 1/2**). Alternatively, you can **download a zipped collection of Xenium readouts** from [here](https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/SpatialDataAlignment/Xenium_vs_Visium/10X_Xenium_Visium.zip). 

<br>

## Building VoltRon objects

VoltRon includes built-in functions for converting readouts of Xenium experiments into VoltRon objects. The **importXenium** function locates all readout documents under the output folder of the Xenium experiment, and forms a VoltRon object. We will import both Xenium replicates separately, and merge them after some image manipulation.

```{r eval = FALSE, class.source="watch-out"}

library(VoltRon)

Xen_R1 <- importXenium("Xenium_R1/outs", sample_name = "XeniumR1", import_molecules = TRUE)

Xen_R2 <- importXenium("Xenium_R2/outs", sample_name = "XeniumR2", import_molecules = TRUE)

```

Before moving on to the downstream analysis of the imaging-based data, we can inspect both Xenium images. We use the **vrImages** function to call and visualize reference images of all VoltRon objects. Observe that the DAPI image of the second Xenium replicate is dim, hence we might need to increase the brightness.  

```{r eval = FALSE, class.source="watch-out"}

vrImages(Xen_R1)

vrImages(Xen_R2)

```

<table>

<tbody>

  <tr style = "vertical-align: center">

  <td style = "width:50%; vertical-align: center"> <img src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/xeniumr1.png" class="center"></td>

  <td style = "width:50%; vertical-align: center"> <img src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/xeniumr2.png" class="center"></td>

  </tr>

</tbody>

</table>

<br>

We can adjust the brightness of the second Xenium replicate using the **modulateImage** function where we can change the brightness and saturation of the reference image of this VoltRon object. This functionality is optional for VoltRon objects and should be used when images require further adjustments. 

```{r eval = FALSE, class.source="watch-out"}

Xen_R2 <- modulateImage(Xen_R2, brightness = 800)

vrImages(Xen_R2)

```

<img width="40%" height="40%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/xeniumr2_new.png" class="center">

<br>

Once both VoltRon objects are created and images are well-tuned, we can merge these two into a single VoltRon object. 

```{r eval = FALSE, class.source="watch-out"}

Xen_list <- list(Xen_R1, Xen_R2)

Xen_data <- merge(Xen_list[[1]], Xen_list[-1])

```

```

VoltRon Object 

XeniumR1: 

  Layers: Section1 

XeniumR2: 

  Layers: Section1 

Assays: Xenium(Main) 

```

<br>

## Spatial Visualization

With **vrSpatialPlot**, we can visualize Xenium experiments in both cellular and subcellular context. Since we have not yet started analyzing raw counts of cells, we can first visualize some transcripts of interest. We first visualize mRNAs of ACTA2, a marker for smooth muscle cell actin, and TCF7, an early exhausted t cell marker. 

We can interactively select a subset of interest within the tissue section and visualize the localization of these transcripts. Here we subset a ductal carcinoma niche, and visualize visualize mRNAs of **(i)** ACTA2, a marker for smooth muscle cell actin, and **(ii)** TCF7, an early exhausted t cell marker.

```{r eval = FALSE, class.source="watch-out"}

Xen_R1_subsetinfo <- subset(Xen_R1, interactive = TRUE)

Xen_R1_subset <- Xen_R1_subsetinfo$subsets[[1]]

vrSpatialPlot(Xen_R1_subset, assay = "Xenium_mol", group.by = "gene",

              group.id = c("ACTA2", "KRT15", "TACSTD2", "CEACAM6"), pt.size = 0.2, legend.pt.size = 5)

```

<img width="70%" height="70%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_transcripts_visualize.png" class="center">

We can also visualize count data of cells in the Xenium replicates. The behaviour of **vrSpatialFeaturePlot** (and most plotting functions in VoltRon) depend on the number of assays associated with the assay type (e.g. Xenium is both cell and subcellular type). Here, we have two assays, and we visualize two features, hence the resulting plot would include four panels. Prior to spatial visualization, we can normalize the counts to correct for count depth of cells by **(i)** dividing counts with total counts in each cell, **(ii)** multiply with some constant (default: 10000), and followed by **(iii)** log transformation of the counts.  

```{r eval = FALSE, class.source="watch-out", fig.align='center'}

Xen_data <- normalizeData(Xen_data, sizefactor = 1000)

vrSpatialFeaturePlot(Xen_data, features = c("ACTA2", "TCF7"), alpha = 1, pt.size = 0.7)

```

<img width="90%" height="90%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_spatialfeature_xenium.png" class="center">

## Processing and Embedding

Some number of cells in both Xenium replicates might have extremely low counts. Although cells are detected at these locations, the low total counts of cells would make it challenging for phenotyping and clustering these cells. Hence, we remove such cells from the VoltRon objects.

```{r eval = FALSE, class.source="watch-out"}

Xen_data <- subset(Xen_data, Count > 5)

```

VoltRon is capable of reducing dimensionality of datasets using both PCA and UMAP which we gonna use to build profile-specific neighborhood graphs and partition the data into cell types. 

```{r eval = FALSE, class.source="watch-out"}

Xen_data <- getPCA(Xen_data, dims = 20)

Xen_data <- getUMAP(Xen_data, dims = 1:20)

```

We can also visualize the normalized expression of these features on embedding spaces (e.g. UMAP) using **vrEmbeddingFeaturePlot** function. 

```{r eval = FALSE, class.source="watch-out", fig.align='center'}

vrEmbeddingFeaturePlot(Xen_data, features = c("LRRC15", "TCF7"), embedding = "umap", 

                       pt.size = 0.4)

```

<img width="100%" height="100%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_featureplot_xenium.png" class="center">

<br>

## Clustering

Next, we build neighborhood graphs with the **shared nearest neighbors (SNN)** of cells which are constructed from dimensionally reduced gene expression profiles. The function **getProfileNeighbors** also has an option of building **k-nearest neighbors (kNN)** graphs. 

```{r eval = FALSE, class.source="watch-out", fig.align='center'}

Xen_data <- getProfileNeighbors(Xen_data, dims = 1:20, method = "SNN")

vrGraphNames(Xen_data)

```

```

[1] "SNN"

```

We can later conduct a clustering of cells using the **leiden's method** from the igraph package, which is utilized with the **getClusters** function.

```{r eval = FALSE, class.source="watch-out", fig.align='center'}

Xen_data <- getClusters(Xen_data, resolution = 1.0, label = "Clusters", graph = "SNN")

```

Now we can label each cell with the associated clustering index and take a look at the clustering accuracy on the embedding space, and we can also visualize these clusters on a spatial context. 

```{r eval = FALSE, class.source="watch-out", fig.align='center'}

vrEmbeddingPlot(Xen_data, group.by = "Clusters", embedding = "umap", 

                pt.size = 0.4, label = TRUE)

```

<img width="60%" height="60%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_embedplot_xenium.png" class="center">

```{r eval = FALSE, class.source="watch-out", fig.align='center'}

vrSpatialPlot(Xen_data, group.by = "Clusters", pt.size = 0.18, background.color = "black")

```

<img width="100%" height="100%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_spatial_xenium.png" class="center">

<br>

## Annotation

We can annotate each of these clusters according to their positive markers across 313 features. One can use the **FindAllMarkers** from the [Seurat](https://satijalab.org/seurat/) package to pinpoint these markers by first utilizing the **as.Seurat** function first on the Xenium assays of the VoltRon object. 

For more information on conversion to other packages, please visit the [Converting VoltRon Objects](conversion.html).

Let us create a new metadata feature from the **Clusters** column, called **CellType**, we can insert this new metadata column directly to the object.

```{r eval = FALSE, class.source="watch-out", fig.align='center'}

clusters <- factor(Xen_data$Clusters, levels = sort(unique(Xen_data$Clusters)))

levels(clusters) <- c("DCIS_1", 

                      "DCIS_2", 

                      "CD4_TCells",

                      "Adipocytes",

                      "PLD4+_LILRA4+_CD4+_Cells",

                      "ACTA2_myoepithelial",

                      "IT_2",

                      "Macrophages",

                      "MastCells",

                      "Bcells",

                      "StromalCells",

                      "CD8_TCells",

                      "CD8_TCells",

                      "EndothelialCells",

                      "StromalCells",

                      "MyelomaCells",

                      "IT_1",

                      "IT_2",

                      "ACTA2_myoepithelial",

                      "DCIS_2",

                      "IT_3",

                      "KRT15_myoepithelial")

Xen_data$CellType <- as.character(clusters)

```

**vrSpatialPlot** function can visualize multiple types of metadata columns, and users can change the location of the legends as well.

```{r eval = FALSE, class.source="watch-out", fig.align='center'}

vrSpatialPlot(Xen_data, group.by = "CellType", pt.size = 0.13, background.color = "black", 

              legend.loc = "top", n.tile = 500)

```

<img width="100%" height="100%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_spatial_xenium_annotated.png" class="center">

<br>

# Visium Data Analysis

Spot-based spatial transcriptomic assays capture spatially-resolved gene expression profiles that are somewhat closer to single cell resolution. However, each spot still include a few number of cells that are likely from a combination of cell types within the tissue of origin. VoltRon analyzes spot level spatial data sets and even allows selecting a highly variable subset of features to cluster spots into meaningful groups of in situ spots for detecting niches of interests

## Import ST Data

For this tutorial we will analyze spot-based transcriptomic assays from Mouse Brain generated by the [Visium](https://www.10xgenomics.com/products/spatial-gene-expression) instrument. 

You can find and download readouts of all four Visium sections [here](https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Cellanalysis/Visium/MouseBrainSerialSections.zip). The **Mouse Brain Serial Section 1** datasets can be downloaded from [here](https://www.10xgenomics.com/resources/datasets?menu%5Bproducts.name%5D=Spatial%20Gene%20Expression&query=&page=1&configure%5BhitsPerPage%5D=50&configure%5BmaxValuesPerFacet%5D=1000) (specifically, please filter for **Species=Mouse**, **AnatomicalEntity=brain**, **Chemistry=v1** and **PipelineVersion=v1.1.0**).

We will now import each of four samples separately and merge them into one VoltRon object. There are four sections in total given two serial anterior and serial posterior sections, hence we have **two tissue blocks each having two layers**. 

```{r eval = FALSE, class.source="watch-out"}

library(VoltRon)

Ant_Sec1 <- importVisium("Sagittal_Anterior/Section1/", sample_name = "Anterior1")

Pos_Sec1 <- importVisium("Sagittal_Posterior/Section1/", sample_name = "Posterior1")

# merge datasets

MBrain_Sec <- merge(Ant_Sec1, Pos_Sec1, samples = c("Anterior", "Posterior"))

MBrain_Sec

```

```

VoltRon Object 

Anterior: 

  Layers: Section1 

Posterior: 

  Layers: Section1

Assays: Visium(Main) 

```

VoltRon maps metadata features on the spatial images, multiple features can be provided for all assays/layers associated with the main assay (Visium). 

```{r eval = FALSE, class.source="watch-out"}

vrSpatialFeaturePlot(MBrain_Sec, features = "Count", crop = TRUE, alpha = 1, ncol = 2)

```

<img width="80%" height="80%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_visium_firstplot.png" class="center">

<br>

## Feature Selection

VoltRon captures the nearly full transcriptome of the Visium data which then can be filtered from a list of features ranked by their variance and importance. We use the **variance stabilization transformation (vst)** on each individual assay using the **getFeatures** function and combine these ranked list to capture features important for all assay of the Visium data later with **getVariableFeatures** function.

```{r eval = FALSE, class.source="watch-out"}

head(vrFeatures(MBrain_Sec))

```

```

[1] "Xkr4"    "Gm1992"  "Gm19938" "Gm37381" "Rp1"     "Sox17"  

```

```{r eval = FALSE, class.source="watch-out"}

length(vrFeatures(MBrain_Sec))

```

```

[1] 33502

```

```{r eval = FALSE, class.source="watch-out"}

MBrain_Sec <- normalizeData(MBrain_Sec)

MBrain_Sec <- getFeatures(MBrain_Sec, n = 3000)

head(vrFeatureData(MBrain_Sec))

```

```

                mean          var    adj_var  rank

Xkr4    0.0248608534 0.0249941807 0.02800216 14114

Gm1992  0.0000000000 0.0000000000 0.00000000     0

Gm19938 0.0285714286 0.0322197476 0.03224908 13889

Gm37381 0.0000000000 0.0000000000 0.00000000     0

Rp1     0.0003710575 0.0003710575 0.00000000     0

Sox17   0.1907235622 0.2219629135 0.23715920 10304

```

```{r eval = FALSE, class.source="watch-out"}

selected_features <- getVariableFeatures(MBrain_Sec)

head(selected_features, 20)

```

```

[1] "Bc1"     "mt-Co1"  "mt-Co3"  "mt-Atp6" "mt-Co2"  "mt-Cytb" "mt-Nd4"  "mt-Nd1"  "mt-Nd2"  

[2] "Fth1"    "Hbb-bs"  "Cst3"    "Gapdh"   "Tmsb4x"  "Mbp"     "Rplp1"   "Ttr"     "Ppia"    

[3] "Ckb"     "mt-Nd3" 

```

## Embedding

Now we can learn and visualize PCA and UMAP embeddings on this smaller number of selected features

```{r eval = FALSE, class.source="watch-out"}

MBrain_Sec <- getPCA(MBrain_Sec, features = selected_features, dims = 30)

MBrain_Sec <- getUMAP(MBrain_Sec, dims = 1:30)

vrEmbeddingPlot(MBrain_Sec, embedding = "umap")

```

<img width="65%" height="65%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_visium_umap.png" class="center">

<br>

## Clustering

```{r eval = FALSE, class.source="watch-out"}

MBrain_Sec <- getProfileNeighbors(MBrain_Sec, dims = 1:30, k = 10, method = "SNN")

vrGraphNames(MBrain_Sec)

```

```

[1] "SNN"

```

```{r eval = FALSE, class.source="watch-out"}

MBrain_Sec <- getClusters(MBrain_Sec, resolution = 0.5, label = "Clusters", graph = "SNN")

vrEmbeddingPlot(MBrain_Sec, embedding = "umap", group.by = "Clusters")

```

<img width="65%" height="65%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_visium_umap_clusters.png" class="center">

<br>

```{r eval = FALSE, class.source="watch-out"}

vrSpatialPlot(MBrain_Sec, group.by = "Clusters")

```

<img width="65%" height="65%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_visium_spatial_clusters1.png" class="center">

<br>

# MELC Data Analysis

VoltRon also provides support for imaging based proteomics assays. In this next use case, we analyze cells characterized by **multi-epitope ligand cartography (MELC)** with a panel of 44 parameters. We use the already segmented cells on which expression of **43 protein features** (excluding DAPI) were mapped to these cells. 

We use the segmented cells over microscopy images collected from **control** and **COVID-19** lung tissues of donors categorized based on disease durations (**control**, **acute**, **chronic** and **prolonged**). Each image is associated with one of few field of views (FOVs) from a single tissue section of a donor. See [GSE190732](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE190732) for more information. You can download the **IFdata.csv** file and the folder with the **DAPI** images [here](https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Cellanalysis/MELC/GSE190732.zip). 

We import the **protein intensities**, **metadata** and **coordinates** associated with segmented cells across FOVs of samples.

```{r eval = FALSE, class.source="watch-out"}

library(VoltRon)

IFdata <- read.csv("IFdata.csv")

data <- IFdata[,c(2:43)]

metadata <- IFdata[,c("disease_state", "object_id", "cluster", "Clusters",

                      "SourceID", "Sample", "FOV", "Section")]

coordinates <- as.matrix(IFdata[,c("posX","posY")], rownames.force = TRUE)

```

<br>

## Importing MELC data

Before analyzing MELC assays across FOVs, we should **build a VoltRon object** for each individual FOV/Section by using the **formVoltron** function. We then merge these sections to respective tissue blocks by defining their samples of origins. We can also define **assay names**, **assay types** and **sample (i.e. block) names** of these objects. 

```{r eval = FALSE, class.source="watch-out", fig.align='center'}

library(dplyr)

library(magick)

vr_list <- list()

sample_metadata <- metadata %>% select(Sample, FOV, Section) %>% distinct()

for(i in 1:nrow(sample_metadata)){

  vrassay <- sample_metadata[i,]

  cells <- rownames(metadata)[metadata$Section == vrassay$Section]

  image <- image_read(paste0("DAPI/", vrassay$Sample, "/DAPI_", vrassay$FOV, ".tif"))

  vr_list[[vrassay$Section]] <- formVoltRon(data = t(data[cells,]), 

                                            metadata = metadata[cells,],

                                            image = image, 

                                            coords = coordinates[cells,],

                                            main.assay = "MELC", 

                                            assay.type = "cell",

                                            sample_name = vrassay$Section)

}

```

Before moving forward with merging FOVs, we should **flip coordinates** of cells and perhaps also then **resize** these images. The main reason for this coordinate flipping is that the y-axis of most digital images are of the opposite direction to the commonly used coordinate spaces. 

```{r eval = FALSE, class.source="watch-out", fig.align='center'}

for(i in 1:nrow(sample_metadata)){

  vrassay <- sample_metadata[i,]

  vr_list[[vrassay$Section]] <- flipCoordinates(vr_list[[vrassay$Section]])

  vr_list[[vrassay$Section]] <- resizeImage(vr_list[[vrassay$Section]], size = 600)

}

```

Finally, we merge these assays into one VoltRon object. The **samples** arguement in the merge function determines which assays are layers of a single tissue sample/block.

```{r eval = FALSE, class.source="watch-out", fig.align='center'}

vr_merged <- merge(vr_list[[1]], vr_list[-1], samples = sample_metadata$Sample)

vr_merged 

```

```

VoltRon Object 

control_case_3: 

  Layers: Section1 Section2 

control_case_2: 

  Layers: Section1 Section2 

control_case_1: 

  Layers: Section1 Section2 Section3 

acute_case_3: 

  Layers: Section1 Section2 

acute_case_1: 

  Layers: Section1 Section2 

... 

There are 13 samples in total 

Assays: MELC(Main) 

```

<br>

The prolonged case 4 has two fields of views (FOVs). By subsetting on the sample of a prolonged case, we can visualize only these two sections, and visualize the protein expression of CD31 and Pancytokeratin which are markers of endothelial and epithelial cells.

```{r eval = FALSE, class.source="watch-out", fig.align='center'}

vr_subset <- subset(vr_merged, samples = "prolonged_case_4")

g1 <- vrSpatialFeaturePlot(vr_subset, features = c("CD31", "Pancytokeratin"), alpha = 1, 

                           pt.size = 0.7, background.color = "black")

```

<img width="90%" height="90%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_spatialfeature.png" class="center">

<br>

## Dimensionality Reduction

We can utilize dimensional reduction of the available protein markers using the getPCA and getUMAP functions, but now with relatively lower numbers of principal components which are enough to capture the information across 44 features.

```{r eval = FALSE, class.source="watch-out", fig.align='center'}

vr_merged <- getPCA(vr_merged, dims = 10)

vr_merged <- getUMAP(vr_merged, dims = 1:10)

vrEmbeddingFeaturePlot(vr_merged, features = c("CD31", "Pancytokeratin"), embedding = "umap")

```

<img width="100%" height="100%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_embedding.png" class="center">

<br>

## Clustering

Now we can visualize the clusters across these sections and perhaps also check for clusters that may reside in only specific disease conditions. 

```{r eval = FALSE, class.source="watch-out", fig.align='center'}

# SNN graph and clusters

vr_merged <- getProfileNeighbors(vr_merged, dims = 1:10, k = 10, method = "SNN")

vrGraphNames(vr_merged)

```

```

[1] "SNN"

```

```{r eval = FALSE, class.source="watch-out", fig.align='center'}

vr_merged <- getClusters(vr_merged, resolution = 0.8, label = "MELC_Clusters", graph = "SNN")

# install patchwork package

if (!requireNamespace("patchwork", quietly = TRUE))

  install.packages("patchwork")

library(patchwork)

# visualize conditions and clusters

vr_merged$Condition <- gsub("_[0-9]$", "", vr_merged$Sample)

g1 <- vrEmbeddingPlot(vr_merged, group.by = c("Condition"), embedding = "umap")

g2 <- vrEmbeddingPlot(vr_merged, group.by = c("MELC_Clusters"), embedding = "umap", 

                      label = TRUE)

g1 | g2

```

<img width="100%" height="100%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_embeddingclusters.png" class="center">

<br>

## Visualization of Markers

VoltRon provides both violin plots (**vrViolinPlot**) and heatmaps (**vrHeatmapPlot**) to further investigate the enrichment of markers across newly clustered datasets. **Note:** the vrHeatmapPlot function would require you to have the **ComplexHeatmap** package in your namespace. 

```{r eval = FALSE, class.source="watch-out", fig.align='center'}

# install patchwork package

if (!requireNamespace("ComplexHeatmap", quietly = TRUE))

  BiocManager::install("ComplexHeatmap")

library(ComplexHeatmap)

# Visualize Markers

vrHeatmapPlot(vr_merged, features = vrFeatures(vr_merged), 

              group.by = "MELC_Clusters", show_row_names = TRUE)

```

<img width="80%" height="80%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_heatmapclusters.png" class="center">

<br>

```{r eval = FALSE, class.source="watch-out", fig.align='center'}

vrViolinPlot(vr_merged, features = c("CD3", "SMA", "Pancytokeratin", "CCR2"), 

             group.by = "MELC_Clusters", ncol = 2)

```

<img width="80%" height="80%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_violinclusters.png" class="center">

<br>

## Neighborhood Analysis

We use the **vrNeighbourhoodEnrichment** function to detect cell type pairs that co occur within each others' neighborhoods. First, we establish **spatial neighborhood graphs** that determine the neighbors of each cell on tissue sections. 

[Delaunay tesselations](https://en.wikipedia.org/wiki/Delaunay_triangulation) or graphs are commonly used to determine neighbors of spatial entities. The function **getSpatialNeighbors** builds a delaunay graph of all assays of a certain type and detects neighbors of cells in a VoltRon object.

```{r eval = FALSE, class.source="watch-out", fig.align='center'}

vr_merged <- getSpatialNeighbors(vr_merged, method = "delaunay")

```

The graph **delaunay**, which we will use for spatially-aware neightborhood analysis, is now the second graph available in the VoltRon object along with **SNN**.

```{r eval = FALSE, class.source="watch-out", fig.align='center'}

vrGraphNames(vr_merged)

```

```

[1] "SNN"      "delaunay"

```

Once neighbors are founds, we can apply a **permutation test** that compares the number of cell type occurances with an expected number of these occurances under multiple permutations of labels in the tissue (fixed coordinates but cells are randomly labelled). A similar approach is used to by several spatial analysis frameworks and packages ([Schapiro et. al 2017](https://www.nature.com/articles/nmeth.4391), [Palla et. al 2022](https://www.nature.com/articles/s41592-021-01358-2)).

Here, we will use the original cell type labels annotated by [Mothes et. al 2023](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9922044/).

```{r eval = FALSE, class.source="watch-out", fig.align='center'}

neighborhood_results <- vrNeighbourhoodEnrichment(vr_merged, group.by = "Clusters")

```

The neighborhood analysis provides the results of:

* the **association** tests (whether cell types are within each other's neighborhood)

* the **segregation** tests (whether cell types are clustered separately) 

between all cell type pairs across each layers and assay. 

The number of each cell in a pair in each section is reported to assess the impact of the results of the test (i.e. low number of abundance in one cell type may indicate low impact).

```{r eval = FALSE, class.source="watch-out", fig.align='center'}

head(neighborhood_results)

```

<div><pre><code style="font-size: 10px;">          from_value     to_value   p_assoc   p_segreg p_assoc_adj p_segreg_adj n_from n_to AssayID Assay    Layer         Sample

Assay1.1 CD163+ macs  CD163+ macs 0.0000000 1.00000000      0.0000   1.00000000     41   41  Assay1  MELC Section1 control_case_3

Assay1.2 CD163+ macs CD4+ T cells 0.9380000 0.03300000      0.9980   0.09762866     41   48  Assay1  MELC Section1 control_case_3

Assay1.3 CD163+ macs  CD8+ Tcells 0.8779011 0.04339051      0.9980   0.09762866     41   11  Assay1  MELC Section1 control_case_3

Assay1.4 CD163+ macs     NK cells 0.8190000 0.08700000      0.9980   0.15660000     41   15  Assay1  MELC Section1 control_case_3

Assay1.5 CD163+ macs   endothelia 0.1230000 0.85100000      0.5535   0.95737500     41  139  Assay1  MELC Section1 control_case_3

Assay1.6 CD163+ macs    epithelia 0.9320000 0.03600000      0.9980   0.09762866     41   39  Assay1  MELC Section1 control_case_3</code></pre></div>

<br>

```{r eval = FALSE, class.source="watch-out", fig.align='center'}

vrNeighbourhoodEnrichmentPlot(neighborhood_results, assay = "Assay1", type = "assoc")

```

<img width="70%" height="70%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_neighenrichment.png" class="center">