---

title: "Importing Spatial Data"

output: 

  html_document:

    toc: true

    toc_depth: 4

    toc_float:

      collapsed: false

      smooth_scroll: false

---

<style>

.title{

  display: none;

}

body {

  text-align: justify

}

.center {

  display: block;

  margin-left: auto;

  margin-right: auto;

}

table, th, td {

  border-collapse: collapse;

  align-self: center;

  padding-right: 10px;

  padding-left: 10px;

}

</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)

library(VoltRon)

library(patchwork)

library(magick)

if(!requireNamespace("ggnewscale"))

  install.packages("ggnewscale")

```

<br>

# Working with VoltRon objects

In this tutorial, we will cover some of the fundemantel built-in functions of VoltRon to manage data, images and spatial information of assays. Here is the list of example datasets included in the Package.

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

data("visium_data")

data("melc_data")

data("xenium_data")

data("merged_object")

```

<br>

## Sample Metadata {#samplemetadata}

The sample metadata is a summarized data frame which provides informations on assays, layers and sample (tissue blocks) included in the VoltRon object. You can use the **SampleMetadata** function to call this table. The row names are the unique assay IDs, following assay types (e.g. Visium), layer name and Sample (typicall tissue block with multiple layers and assays).

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

SampleMetadata(visium_data)

```

<br>

## Metadata

The **Metadata** function calls the metadata table of the main assay type (see [vrMainAssay](#assays)).

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

Metadata(visium_data)

```

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

head(Metadata(visium_data),6)

```

You can also specify the assay ID or the assay type to call metadata of a subset of spatial points.

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

Metadata(visium_data, assay = "Visium")

Metadata(visium_data, assay = "Assay1")

```

The **type** argument returns all spatial points of a given spatial point type (ROI, spot, cell, molecule or tile).

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

Metadata(visium_data, type = "spot")

```

If the type argument is specified as **all**. Then this would return the VoltRon Metadata Object (vrMetadata)

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

Metadata(visium_data, type = "all")

```

<br>

## Spatial Points {#spatialpoints}

In a VoltRon object, spatial points are spatially defined entities with coordinates, segments which are of any of five assay types (ROI, spot, cell, molecule and tile).

The **vrSpatialPoints** function return the IDs of these entities for further downstream operations, such as subsetting etc.

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

vrSpatialPoints(visium_data)

```

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

head(vrSpatialPoints(visium_data),3)

```

You can also specify the assay ID or the assay type to call metadata of a subset of spatial points.

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

vrSpatialPoints(visium_data, assay = "Visium")

vrSpatialPoints(visium_data, assay = "Assay1")

```

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

head(vrSpatialPoints(visium_data, assay = "Visium"),3)

```

<br>

## Assays {#assays}

The default (or main) assay of the VoltRon object is typically shown when printed (next to the assay name says "main" in paranthesis)

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

melc_data

```

You can also call/get the name of the default (or main) using the **vrMainAssay** function

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

vrMainAssay(melc_data)

```

You can also set the main assay yourself, but only the assay types given in the **Assay** column of **SampleMetadata(object)**.

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

SampleMetadata(melc_data)

```

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

vrMainAssay(melc_data) <- "MELC"

```

You can also get the assay IDs associated with the main assay using **vrAssayNames** function. 

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

vrAssayNames(melc_data)

```

The assay type can be provided with the **assay** arguement to get type specific assay IDs

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

vrAssayNames(melc_data, assay = "MELC")

```

<br>

## Coordinates and Segments

The **vrCoordinates** function is used to call the centroids of spots, cells and all other spatial points types with VoltRon objects. 

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

vrCoordinates(visium_data)

```

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

head(vrCoordinates(visium_data),6)

```

You can also specify the assay ID or the assay type to call coordinates of a subset of spatial points.

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

vrCoordinates(visium_data, assay = "Visium")

vrCoordinates(visium_data, assay = "Assay1")

```

<br>

Each assay a VoltRon object may incorporate indefinite number of coordinate systems. One can look for these coordinate systems using the **vrMainSpatial** function, and select one of systems to call coordinates (see [vrMainSpatial](#image))

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

vrCoordinates(visium_data, spatial_name = "main")

```

The **reg** option in the **vrCoordinates** function looks for a registered version of the main coordinate system and returns its coordinates (if there is any).

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

vrCoordinates(visium_data, spatial_name = "main", reg = TRUE)

```

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

head(vrCoordinates(visium_data, spatial_name = "main", reg = TRUE),6)

```

<br>

The arguements of the **vrSegments** functions are idential to vrCoordinates and return a list of polygon corners associated with the coordinate system of the coordinates. 

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

vrSegments(xenium_data)

vrSegments(xenium_data, assay = "Xenium")

vrSegments(xenium_data, assay = "Assay1")

vrSegments(xenium_data, spatial_name = "main")

vrSegments(xenium_data, spatial_name = "main", reg = TRUE)

```

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

head(vrSegments(xenium_data),2)

```

<br>

## Images and Coordinate Systems {#image}

Below we review some of the essential built-in functions to manipulate and manage images of a VoltRon object as well as the coordinate systems that are associated with these images.

<br>

### Spatial Coordinate Systems

In a VoltRon object, each image has a unique ID which is also associated with a coordinate (or spatial) system. The names of these spatial systems can be called using the **vrSpatialNames** function.

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

vrSpatialNames(visium_data)

```

You can also specify the assay ID or the assay type to call image names of a subset assays or any assay type (see [Sample Metadata](#samplemetadata))

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

vrSpatialNames(visium_data, assay = "Visium")

vrSpatialNames(visium_data, assay = "Assay1")

```

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

vrSpatialNames(visium_data, assay = "Visium")

```

<br>

If you wanna see the list of all spatial systems associated with all assays of the main assay type, you can call the **vrMainSpatial** function.

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

vrMainSpatial(visium_data)

```

You can also specify the assay ID or the assay type to call spatial system names of a subset assays or any assay type (see [Sample Metadata](#samplemetadata))

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

vrMainSpatial(visium_data, assay = "Visium")

vrMainSpatial(visium_data, assay = "Assay1")

```

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

vrMainSpatial(visium_data, assay = "Visium")

```

<br>

### Channel Names

Each spatial system object (vrSpatial) in VoltRon object can also include a indefinite number of channels which you can get (or request) a list of for further use.

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

vrImageChannelNames(melc_data)

vrImageChannelNames(melc_data, assay = "MELC")

vrImageChannelNames(melc_data, assay = "Assay1")

```

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

vrImageChannelNames(melc_data)

```

<br>

### Get and Set Images

In VoltRon, images can be called specifically, or return as a list. The return image is of a magick image object (see [magick](https://cran.r-project.org/web/packages/magick/vignettes/intro.html) package)

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

vrImages(visium_data)

vrImages(visium_data, assay = "Visium")

vrImages(visium_data, assay = "Assay1")

```

```{r echo = FALSE, class.source="watch-out", fig.align='center', fig.height=5, fig.width=6, out.width="60%"}

magick::image_ggplot(vrImages(visium_data))

```

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

<br>

Once you know the name of a specific channel, you can the image of a specific channel by providing the name and the associated channel.

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

vrImages(melc_data, name = "MELC", channel = "DAPI")

```

```{r echo = FALSE, class.source="watch-out", fig.align='center', fig.height=5, fig.width=6, out.width="60%"}

magick::image_ggplot(vrImages(melc_data, name = "MELC", channel = "DAPI"))

```

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

<br>

You can set up the main channel as well as the main spatial system name for later use.

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

vrMainSpatial(melc_data, assay = "Assay1") <- c("MELC", "CD45")

vrImages(melc_data)

```

```{r echo = FALSE, class.source="watch-out", fig.align='center', fig.height=5, fig.width=6, out.width="60%"}

vrMainSpatial(melc_data, assay = "Assay1") <- c("MELC", "CD45")

magick::image_ggplot(vrImages(melc_data))

```

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

<br>

You can also resize the images as they are being returned. This is usually used for visualization purposes and helps speeding up visualization for large images. This is accomplished with **scale.perc** arguement.

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

vrImages(melc_data, scale.perc = 25)

```

<br>

### Combining Image Channels

VoltRon even allows manipulation of channel images if you also provide an associate list of colors.

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

melc_data <- combineChannels(melc_data, 

                             channels = c("DAPI", "CD45"), colors = c("grey", "green"), 

                             channel_key = "combined")

```

These new images can be stored as new channels within the same image object, and called later again

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

vrImageChannelNames(melc_data)

```

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

vrImages(melc_data, channel = "combined")

```

```{r echo = FALSE, class.source="watch-out", fig.align='center', fig.height=5, fig.width=6, out.width="60%"}

magick::image_ggplot(vrImages(melc_data, channel = "combined"))

```

<br>

## Feature Matrix (Data)

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

vrData(visium_data)[3:8,3:5]

```

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

vrData(visium_data, norm = TRUE)[3:8,3:5]

```

<br>

## Embeddings {#embedding}

You can parse and even set individual embedding elements in a VoltRon object.

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

vrEmbeddingNames(xenium_data)

```

You can use these names of get the associated embedding dataset from the object. 

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

vrEmbeddings(xenium_data, type = "umap")

```

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

head(vrEmbeddings(xenium_data, type = "umap"),6)

```

You can also set and create new embedding elements in the voltron object. In this case, you have to make sure that the row names should match with the targeted spatial points.

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

new_umap_data <- vrEmbeddings(xenium_data, type = "umap")

vrEmbeddings(xenium_data, type = "new_umap") <- new_umap_data*2

```

Now we can observe changes to the new embedding data. 

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

vrEmbeddings(xenium_data, type = "new_umap")

```

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

head(vrEmbeddings(xenium_data, type = "new_umap"),6)

```

You can choose individual assay names or assay classes.

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

vrEmbeddings(xenium_data, type = "umap", assay = "Xenium")

vrEmbeddings(xenium_data, type = "pca", assay = "Assay1")

```

<br>

## Subsetting VoltRon objects

<br>

### sample/assay 

VoltRon object can be subsetted in a variety of ways using assay names, sample names, spatial point names, features (e.g. gene), image coordinates (crop boxes or bounding boxes) as well as interactively.

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

visium_data_subset <- subset(visium_data, assays = "Assay1")

```

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

visium_data_subset <- subset(visium_data, samples = "Anterior1")

visium_data_subset

```

<br>

### spatial points

You can use a list of spatial points (typically using **vrSpatialPoints** function, see [vrSpatialPoints](#spatialpoints)) 

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

selected_points <- vrSpatialPoints(visium_data)

selected_points[1:20]

```

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

selected_points <- vrSpatialPoints(visium_data)

head(selected_points,3)

```

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

visium_data_subset <- subset(visium_data, spatialpoints = selected_points[1:20])

visium_data_subset

```

<br>

### features

You can select a few number of features and subset the features given this list. However, it would only subset the main assay (see [vrMainAssay](#assays))

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

selected_features <- vrFeatures(visium_data)

selected_features[1:20]

```

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

visium_data_subset <- subset(visium_data, features = selected_features[1:20])

vrFeatures(visium_data_subset)

```

<br>

### Interactive subsetting

VoltRon allows interactively subsetting spatial data. Using the arguement **interactive = TRUE**, a mini Shiny app is triggered where users can select a bounding box to crop the spatial data.

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

visium_data_subset_info <- subset(visium_data, interactive = TRUE)

visium_data_subset <- visium_data_subset_info$subsets[[1]]

```

```

VoltRon Object 

Anterior1_subset: 

  Layers: Section1 

Assays: Visium(Main) 

```

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

vrImages(visium_data_subset)

```

<br>

## Visualization

VoltRon provides visualization utilities for both spatial and embedding level visualizations.

<br>

### Spatial Plots

**vrSpatialPlot** is the main function for visualize labels and identities of cells. This information is parsed from the Metadata of the main assay (see [vrMainAssay](#assays)). The users can also specify the assay.

```{r eval = FALSE, class.source="watch-out", fig.align='center', fig.height=5, fig.width=6, out.width="60%"}

vrSpatialPlot(xenium_data, group.by = "clusters")

```

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

<br>

You can also visualize the segments of spatial points and even get the segments transparent

```{r eval = FALSE, class.source="watch-out", fig.align='center', fig.height=5, fig.width=6, out.width="60%"}

vrSpatialPlot(xenium_data, group.by = "clusters", plot.segments = TRUE, alpha = 0.6)

```

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

<br>

The background color can be set to any color.

```{r eval = FALSE, class.source="watch-out", fig.align='center', fig.height=5, fig.width=6, out.width="60%"}

vrSpatialPlot(xenium_data, group.by = "clusters", plot.segments = TRUE, 

              background.color = "white")

```

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

<br>

In cases where there are multiple coordinate systems, you can get the background set to the name of the image (or coordinate system). You can get the name of these images from [vrSpatialNames](#image). 

If also want to select a channel from the same coordinate system, you can set background arguement as a vector of 2 where first is the name of the image (coordinate system) and the other would be channel name. You can get the name of these channels from [vrImageChannelNames](#image).

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

vrSpatialPlot(xenium_data, group.by = "clusters", plot.segments = TRUE,

              spatial = "main")

vrSpatialPlot(xenium_data, group.by = "clusters", plot.segments = TRUE,

              spatial = "main", channel = "DAPI")

```

<br>

If the visualized assay is of a type "spot", you can crop the image to encapsulate the smallest subset that include all spots.

```{r eval = FALSE, class.source="watch-out", fig.align='center', fig.height=5, fig.width=6, out.width="60%"}

vrSpatialPlot(visium_data, group.by = "Sample", crop = TRUE)

```

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

<br>

The **vrSpatialFeaturePlot** functions the same way as vrSpatialPlot but requires extra arguements such as **features** for selecting features, **norm** for normalized expression (default), and **log** for log transformed counts.

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

vrSpatialFeaturePlot(visium_data, features = "Count")

vrSpatialFeaturePlot(visium_data, features = "Stat1", norm = TRUE)

```

```{r eval = FALSE, echo = FALSE, class.source="watch-out", fig.align='center', fig.height=5, fig.width=12}

g1 <- vrSpatialFeaturePlot(visium_data, features = "Count")

g2 <- vrSpatialFeaturePlot(visium_data, features = "Stat1", norm = TRUE)

g1 | g2

```

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

<br>

For all variations of **vrSpatialPlot** and **vrSpatialFeaturePlot** functions above, you can specifiy the assay names or assay class.

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

vrSpatialPlot(visium_data, assay = "Visium", plot.segments = TRUE)

vrSpatialPlot(visium_data, assay = "Assay1", plot.segments = TRUE)

vrSpatialFeaturePlot(visium_data, assay = "Visium", features = "Count")

vrSpatialFeaturePlot(visium_data, assay = "Assay1", features = "Count")

```

<br>

### Spatial Plots (Multi-Layer)

**vrSpatialPlot** function can be paired with the **addSpatialLayer** function to add additional assays of the same spatial coordinate system as the original assay. We use the piping operator "|>" to connect to functions which overlays the results of the first image with the second. **addSpatialLayer** will check if two assays are within the same tissue block and coordinate system, and will overlay these assays if thats the case. 

This utility can be used to visualize molecules along with cell segments. Here we can also color each individual molecule (gene) type.

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

library(ggnewscale)

```

```{r eval = TRUE, class.source="watch-out", fig.height=5, fig.width=6, fig.align='center', out.width="60%"}

vrSpatialPlot(merged_object, plot.segments = TRUE, group.by = "clusters") |>

  addSpatialLayer(merged_object, assay = "Assay2", group.by = "gene", alpha = 1, 

                  colors = list(KRT15 = "blue", KRT14 = "green"))

```

<br>

You can also visualize the segments that represents annotations of the tissue sections along with the cells.

```{r eval = TRUE, class.source="watch-out", fig.height=5, fig.width=6, fig.align='center', out.width="60%"}

vrSpatialPlot(merged_object, plot.segments = FALSE, group.by = "clusters") |>

  addSpatialLayer(merged_object, assay = "Assay3", group.by = "Region", alpha = 0.3, 

                  color = list(Region1 = "blue", Region2 = "yellow"))

```

<br>

### Embedding Plots

Dimensional reduction and embedding features are also possible for the VoltRon objects.

```{r, eval = FALSE, class.source="watch-out", fig.align='center', fig.height=5, fig.width=6, out.width="60%"}

vrEmbeddingPlot(xenium_data, embedding = "umap", group.by = "clusters", label = T)

```

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

<br>

The **vrEmbeddingFeaturePlot** functions the same way as vrEmbeddingPlot but requires extra arguements such as **features** for selecting features, **norm** for normalized expression (default), and **log** for log transformed counts.

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

vrEmbeddingFeaturePlot(xenium_data, embedding = "umap", features = "Count")

vrEmbeddingFeaturePlot(xenium_data, embedding = "umap", features = "KRT5", norm = FALSE)

```

```{r echo=FALSE, eval = FALSE, class.source="watch-out", fig.height=5, fig.width=12, fig.align='center'}

g1 <- vrEmbeddingFeaturePlot(xenium_data, embedding = "umap", features = "Count")

g2 <- vrEmbeddingFeaturePlot(xenium_data, embedding = "umap", features = "KRT5", norm = FALSE)

g1 | g2

```

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

<br>

For all variations of **vrEmbeddingPlot** and **vrEmbeddingFeaturePlot** functions above, you can specifiy the assay names or assay class.

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

vrEmbeddingPlot(xenium_data, assay = "Xenium", embedding = "umap")

vrEmbeddingPlot(xenium_data, assay = "Assay1", embedding = "umap")

vrEmbeddingFeaturePlot(xenium_data, assay = "Xenium", embedding = "umap", features = "Count")

vrEmbeddingFeaturePlot(xenium_data, assay = "Assay1", embedding = "umap", features = "Count")

```