---
title: "Conversion"
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>

# Conversion to Other Platforms

VoltRon is capable of end-to-end spatial data analysis for all levels of spatial resolutions, including those of single cell resolution. However, VoltRon provides a ecosystem friendly infrastructure where VoltRon objects could be transformed into data structures used by popular computational platforms such as [Seurat](https://satijalab.org/seurat/), [Squidpy](https://squidpy.readthedocs.io/en/stable/) and even [Zarr](http://vitessce.io/docs/data-file-types/#anndata-zarr) for interactive spatial data visualizatiob with [Vitessce](http://vitessce.io/).

For both **Seurat (R)** and **Squidpy (Python)**, we analyse readouts of the experiments conducted on example tissue sections analysed by the [Xenium In Situ](https://www.10xgenomics.com/platforms/xenium) platform. For more information on processing and analyzing Xenium datasets, check the [Cell/Spot Analysis](spotanalysis.html) tutorial.

<br>

## Seurat

We will first see how we can transform VoltRon objects into Seurat object and use built-in functions such as **FindAllMarkers** to visualize marker genes of clusters found by VoltRon. You can find the clustered Xenium data using VoltRon [here](https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Cellanalysis/Xenium/Xenium_data_clustered.rds). 

```{r class.source="watch-out", eval = FALSE}
Xen_data <- readRDS("Xenium_data_clustered.rds")
SampleMetadata(Xen_data)
```

```
        Assay    Layer   Sample
Assay1 Xenium Section1 XeniumR1
Assay3 Xenium Section1 XeniumR2
```

<br>

We use the **as.Seurat** function to convert spatial assays of VoltRon into Seurat objects. Here, a Seurat object defines spatial components of cellular and subcellular assays as **FOV** objects, and we use the **type = "image"** argument to convert spatial coordinates of cells and molecules into individual FOV objects for each Xenium assay/layer in the VoltRon object.  

Please check the [Analysis of Image-based Spatial Data in Seurat](https://satijalab.org/seurat/articles/seurat5_spatial_vignette_2) tutorial for more information on analyzing FOV-based spatial data sets with Seurat. 

note: Use VoltRon::as.Seurat to avoid conflict with Seurat package's as.Seurat function


```{r class.source="watch-out", eval = FALSE}
library(Seurat)
Xen_data_seu <- VoltRon::as.Seurat(Xen_data, cell.assay = "Xenium", type = "image")
Xen_data_seu <- NormalizeData(Xen_data_seu)
Xen_data_seu
```

```
An object of class Seurat 
313 features across 283298 samples within 1 assay 
Active assay: Xenium (313 features, 0 variable features)
 1 layers present: counts
 2 dimensional reductions calculated: pca, umap
 2 spatial fields of view present: fov_Assay1 fov_Assay3
```

<br>

### Marker Analysis

Now that we converted VoltRon into a Seurat object, we can pick the **Clusters** metadata column indicating the clustering of cells and test for marker genes of each individual cluster. 

```{r class.source="watch-out", eval = FALSE}
Idents(Xen_data_seu) <- "Clusters"
markers <- FindAllMarkers(Xen_data_seu)
head(markers[order(markers$avg_log2FC, decreasing = TRUE),])
```

<div><pre><code style="font-size: 13px;">         p_val avg_log2FC pct.1 pct.2 p_val_adj cluster   gene
CPA3         0   7.343881 0.977 0.029         0      16   CPA3
CTSG         0   7.114698 0.878 0.011         0      16   CTSG
LILRA4.1     0   6.992717 0.939 0.015         0      19 LILRA4
ADIPOQ       0   6.860190 0.974 0.025         0       5 ADIPOQ
MS4A1        0   6.763083 0.919 0.027         0      17  MS4A1
BANK1        0   6.082192 0.889 0.037         0      17  BANK1</code></pre></div>

<br>

### Visualization

We can now pick top positive markers from each of these clusters prior to visualization.

```{r class.source="watch-out", eval = FALSE}
library(dplyr)
topmarkers <- markers %>%
  group_by(cluster) %>%
  top_n(n = 5, wt = avg_log2FC)
```

Here, VoltRon incorporates the unique markers learned by the **FindAllMarkers** function from Seurat and uses them to visualize the expression of these markers on heatmaps, and now we can also use these markers for annotating the clusters.   

```{r class.source="watch-out", eval = FALSE}
library(ComplexHeatmap)
marker_features <- unique(topmarkers$gene)
vrHeatmapPlot(Xen_data, features = marker_features, group.by = "Clusters", 
              show_row_names = TRUE, font.size = 10)
```

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

<br>

### Convert with Molecule Data

If defined, the **as.Seurat** function may also convert the molecule assay of the VoltRon object into a Seurat FOV object and allow visualizing molecules. You can find the Xenium VoltRon object with the molecule assay [here](https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Cellanalysis/Xenium/Xen_R1.rds). 

```{r class.source="watch-out", eval = FALSE}
Xen_R1 <- readRDS("Xen_R1.rds")
SampleMetadata(Xen_R1)
```

```
            Assay    Layer   Sample
Assay1     Xenium Section1 XeniumR1
Assay2 Xenium_mol Section1 XeniumR1
```

<br>

We define both the cell level assay and the molecule level assay. 

```{r class.source="watch-out", eval = FALSE}
Xen_R1_seu <- as.Seurat(Xen_R1, cell.assay = "Xenium", molecule.assay = "Xenium_mol", type = "image")
```

<br>

Now we can visualize molecules alongside with cells.

```{r class.source="watch-out", eval = FALSE}
ImageDimPlot(Xen_R1_seu, fov = "fovAssay1", molecules = "PGR", group.by = "orig.ident", cols = "lightgrey", mols.size = 1)
```

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

<br>

## SpatialExperiment

VoltRon can also convert objects in [SpatialExperiment](https://www.bioconductor.org/packages/release/bioc/html/SpatialExperiment.html) objects. We are going to use the Xenium data clustered using VoltRon [here](https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Cellanalysis/Xenium/Xenium_data_clustered.rds).

```{r class.source="watch-out", eval = FALSE}
Xen_data <- readRDS("Xenium_data_clustered.rds")
SampleMetadata(Xen_data)
```

We use the **as.SpatialExperiment** function to convert spatial assays of VoltRon into SpatialExperiment objects. Please check the [Introduction to the SpatialExperiment class](https://www.bioconductor.org/packages/release/bioc/vignettes/SpatialExperiment/inst/doc/SpatialExperiment.html) tutorial for more information. 

```{r class.source="watch-out", eval = FALSE}
library(SpatialExperiment)
spe <- as.SpatialExperiment(Xen_data, assay = "Xenium")
```

Here we can parse the image and visualize.

```{r class.source="watch-out", eval = FALSE}
img <- imgRaster(spe, 
                 sample_id = "Assay1", 
                 image_id = "main")
plot(img)
```

<br>

## Squidpy (Anndata, h5ad)

A true ecosystem friendly computational platform should support data types across multiple computing environments. By allowing users to convert VoltRon objects into annotated data matrix formats such as [anndata](https://github.com/scverse/anndata), we can use built-in spatial data analysis methods available on [squidpy](https://squidpy.readthedocs.io/en/stable/).

You can find the clustered and the annotated Xenium data using VoltRon [here](https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Cellanalysis/Xenium/Xenium_data_clustered_annotated.rds). 

Anndata objects wrapped on h5ad files are commonly used by the [scverse](https://www.nature.com/articles/s41587-023-01733-8) ecosystem for single cell analysis which bring together numeruous tools maintained and distributed by a large community effort. Both squidpy and [scanpy](https://scanpy.readthedocs.io/en/stable/) are currently available on scverse. 

```{r class.source="watch-out", eval = FALSE}
Xen_data <- readRDS("Xenium_data_clustered_annotated.rds")
as.AnnData(Xen_data, assay = "Xenium", file = "Xen_adata_annotated.h5ad")
```

<br>

### Configure Squidpy (scverse)

Here, we use the [reticulate](https://rstudio.github.io/reticulate/) package to call **scverse** module in Python through a prebuilt anaconda environment. However, any python installation with the scverse module can be incorporated by reticulate.   

```{r class.source="watch-out", eval = FALSE}
library(reticulate)
use_condaenv("scverse", required = T)
```

We import some other necessary modules such as pandas, scanpy and squidpy.

```{python class.source="watch-out", eval = FALSE}
from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import scanpy as sc
import squidpy as sq
sc.logging.print_header()
```

<br>

### Filter and Normalize 

We read the annotated Xenium object that was saved as an h5ad file using the **as.Anndata** function in VoltRon, and process before analysis. For more information using scanpy and squidpy on Xenium datasets, check the [Analyzing Xenium data](https://squidpy.readthedocs.io/en/stable/notebooks/tutorials/tutorial_xenium.html) tutorial at squidpy webpage. 

```{python class.source="watch-out", eval = FALSE}
adata = sc.read_h5ad("Xen_adata_annotated.h5ad")
adata.layers["counts"] = adata.X.copy()
sc.pp.normalize_total(adata, inplace=True)
sc.pp.log1p(adata)
```

<br>

### Visualize

We use the **squidpy.pl.spatial_scatter** functions available in squidpy to visualize the spatial localization of cell types of both Xenium replicates.

```{python class.source="watch-out", eval = FALSE}
fig, ax = plt.subplots(1, 2, figsize=(10, 7))
sq.pl.spatial_scatter(adata, library_key = "library_id", library_id = "Assay1", 
                      color=["CellType"], shape=None, size=1, img = False, ax=ax[0])
sq.pl.spatial_scatter(adata, library_key = "library_id", library_id = "Assay3", 
                      color=["CellType"], shape=None, size=1, img = False, ax=ax[1])
plt.show(ax)
```

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

<br>

### Neighborhood Enrichment

We can now use high level spatially-aware functions available in squidpy. We first establish spatial neighbors using the delaunay graphs. The spatial graph and distances will be stored under **.obsp** attribute/matrix. 

```{python class.source="watch-out", eval = FALSE}
sq.gr.spatial_neighbors(adata, coord_type="generic", delaunay=True)
adata
```

```
AnnData object with n_obs × n_vars = 283298 × 313
    obs: 'Count', 'Assay', 'Layer', 'Sample', 'Clusters', 'CellType', 'library_id'
    uns: 'log1p', 'spatial_neighbors'
    obsm: 'spatial'
    layers: 'counts'
    obsp: 'spatial_connectivities', 'spatial_distances'
```

We can now conduct the permutation test for neighborhood enrichment across cell type pairs. 

```{python class.source="watch-out", eval = FALSE}
sq.gr.nhood_enrichment(adata, cluster_key="CellType")
```

```{python class.source="watch-out", eval = FALSE}
fig, ax = plt.subplots(1, 2, figsize=(13, 7))
sq.pl.nhood_enrichment(adata, cluster_key="CellType", figsize=(8, 8), 
                       title="Neighborhood enrichment adata", ax=ax[0])
sq.pl.spatial_scatter(adata, color="CellType", library_key = "library_id", 
                      library_id = "Assay1", shape=None, size=2, ax=ax[1])
plt.show(ax)
```

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

<br>

## Vitessce (Anndata, zarr) 
 
In this section, we will transform VoltRon objects of Xenium data into zarr arrays, and use them for interactive visualization in [Vitessce](http://vitessce.io/). We should first download the vitessceR package which incorporates wrapper function to visualize zarr arrays interactively in R.

```{r class.source="watch-out", eval = FALSE}
install.packages("devtools")
devtools::install_github("vitessce/vitessceR")
```

<br>

You can find the clustered and annotated Xenium data using VoltRon [here](https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Cellanalysis/Xenium/Xenium_data_clustered_annotated.rds). 

```{r class.source="watch-out", eval = FALSE}
Xen_data <- readRDS("Xenium_data_clustered_annotated.rds")
SampleMetadata(Xen_data)
```

```
        Assay    Layer   Sample
Assay1 Xenium Section1 XeniumR1
Assay2 Xenium Section1 XeniumR2
```

<br>

### Interactive Visualization

Now we can convert the VoltRon object into a zarr array using the **as.Zarr** function which will create the array in a specified location.

```{r class.source="watch-out", eval = FALSE}
as.AnnData(Xen_data, assays = "Assay1", 
           file = "xendata_clustered_annotated.zarr", flip_coordinates = TRUE)
```

We can use the zarr file directly in the **vrSpatialPlot** function to visualize the zarr array interactively in Rstudio viewer. The **reduction** arguement allows the umap of the Xenium data to be visualized alongside with the spatial coordinates of the Xenium cells.

```{r class.source="watch-out", eval = FALSE}
vrSpatialPlot("xendata_clustered_annotated.zarr", group.by = "CellType", reduction = "umap")
```

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