588 lines (428 with data), 27.9 kB
---
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">
Datasets
Models
