---

title: "Niche Clustering"

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>

# Spot-based Niche Clustering

Spot-based spatial transcriptomic assays capture spatially-resolved gene expression profiles that are somewhat closer to single cell resolution. However, each spot still includes a few number of cells that are likely originated from few number of cell types, hence transcriptomic profile of each spot would likely include markers from multiple cell types. Here, **RNA deconvolution** can be incorporated to estimate the percentage/abundance of cell types for each spot. We use a scRNAseq dataset as a reference to computationally estimate the relative abundance of cell types across across the spots. 

VoltRon includes wrapper commands for using popular spot-level RNA deconvolution methods such as [RCTD](https://www.nature.com/articles/s41587-021-00830-w) and return estimated abundances as additional feature sets within each layer. These estimated percentages of cell types for each spot could be incorporated to detect **niches** (i.e. small local microenvironments of cells) within the tissue. We can process cell type abundance assays and used them for clustering to detect these niches. 

<br>

## Import Visium 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/2** datasets can be also 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 brain tissue 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")

Ant_Sec2 <- importVisium("Sagittal_Anterior/Section2/", sample_name = "Anterior2")

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

Pos_Sec2 <- importVisium("Sagittal_Posterior/Section2/", sample_name = "Posterior2")

# merge datasets

MBrain_Sec_list <- list(Ant_Sec1, Ant_Sec2, Pos_Sec1, Pos_Sec2)

MBrain_Sec <- merge(MBrain_Sec_list[[1]], MBrain_Sec_list[-1], 

                    samples = c("Anterior", "Anterior", "Posterior", "Posterior"))

MBrain_Sec

```

```

VoltRon Object 

Anterior: 

  Layers: Section1 Section2 

Posterior: 

  Layers: Section1 Section2 

Assays: Visium(Main) 

Features: RNA(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/decon_first_plot.png" class="center">

<br>

## Import scRNA data

We will now import the scRNA data for reference which can be downloaded from [here](https://www.dropbox.com/s/cuowvm4vrf65pvq/allen_cortex.rds?dl=1). Specifically, we will use a scRNA data of Mouse cortical adult brain with 14,000 cells, generated with the SMART-Seq2 protocol, from the Allen Institute. This scRNA data is also used by the Spatial Data Analysis tutorial in [Seurat](https://satijalab.org/seurat/articles/spatial_vignette.html#integration-with-single-cell-data) website.  

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

# install packages if necessary

if(!requireNamespace("Seurat"))

  install.packages("Seurat")

if(!requireNamespace("dplyr"))

  install.packages("dplyr")

# import scRNA data

library(Seurat)

allen_reference <- readRDS("allen_cortex.rds")

# process and reduce dimensionality

library(dplyr)

allen_reference <- SCTransform(allen_reference, ncells = 3000, verbose = FALSE) %>%

  RunPCA(verbose = FALSE) %>%

  RunUMAP(dims = 1:30)

```

Before deconvoluting Visium spots, we correct cell types labels and drop some cell types with extremely few number of cells (e.g. "CR").  

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

# update labels and subset

allen_reference$subclass <- gsub("L2/3 IT", "L23 IT", allen_reference$subclass)

allen_reference <- allen_reference[,colnames(allen_reference)[!allen_reference@meta.data$subclass %in% "CR"]]

# visualize

Idents(allen_reference) <- "subclass"

gsubclass <- DimPlot(allen_reference, reduction = "umap", label = T) + NoLegend()

Idents(allen_reference) <- "class"

gclass <- DimPlot(allen_reference, reduction = "umap", label = T) + NoLegend()

gsubclass | gclass

```

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

<br>

## Spot Deconvolution with RCTD

In order to integrate the scRNA data and the spatial data sets within the VoltRon object and estimate relative cell type abundances for each Visium spot, we will use **RCTD** algorithm which is accessible with the [spacexr](https://github.com/dmcable/spacexr) package. 

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

if(!requireNamespace("spacexr"))

  devtools::install_github("dmcable/spacexr", build_vignettes = FALSE)

```

After running **getDeconvolution**, an additional feature set within the same Visium assay with name **Decon** will be created.

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

library(spacexr)

MBrain_Sec <- getDeconvolution(MBrain_Sec, sc.object = allen_reference, sc.cluster = "subclass", max_cores = 6)

MBrain_Sec

```

```

VoltRon Object 

Anterior: 

  Layers: Section1 Section2 

Posterior: 

  Layers: Section1 Section2 

Assays: Visium(Main) 

Features: RNA(Main) Decon 

```

We can now switch to the **Decon** feature type where features are cell types from the scRNA reference and the data values are cell types percentages in each spot. 

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

vrMainFeatureType(MBrain_Sec) <- "Decon"

vrFeatures(MBrain_Sec)

```

```

 [1] "Astro"      "Endo"       "L23 IT"     "L4"         "L5 IT"      "L5 PT"     

 [7] "L6 CT"      "L6 IT"      "L6b"        "Lamp5"      "Macrophage" "Meis2"     

[13] "NP"         "Oligo"      "Peri"       "Pvalb"      "Serpinf1"   "SMC"       

[19] "Sncg"       "Sst"        "Vip"        "VLMC"          

```

These features (i.e. cell type abundances) can be visualized like any other feature type. 

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

vrSpatialFeaturePlot(MBrain_Sec, features = c("L4", "L5 PT", "Oligo", "Vip"), 

                     crop = TRUE, ncol = 2, alpha = 1, keep.scale = "all")

```

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

<br>

## Clustering 

Relative cell type abundances that are learned by RCTD and stored within VoltRon can now be used to cluster spots. These groups or clusters of spots can often be referred to as **niches**. Here, as a definition, a niche is a region or a collection of regions within tissue that have a distinct cell type composition as opposed to the remaining parts of the tissue.   

The cell type abundances (which adds up to one for each spot) can be normalized and processed like transcriptomic and proteomic profiles prior to clustering (i.e. niche clustering). We treat cell type abundances as [compositional data](https://en.wikipedia.org/wiki/Compositional_data), hence we incorporate **centred log ratio (CLR)** transformation for normalizing them. 

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

vrMainFeatureType(MBrain_Sec) <- "Decon"

MBrain_Sec <- normalizeData(MBrain_Sec, method = "CLR")

```

The CLR normalized assay have only 25 features, each representing a cell type from the single cell reference data. Hence, we can **directly calculate UMAP reductions from this feature abundances** since we dont have much number of features which necessitates dimensionality reduction such as PCA. 

However, we may still need to reduce the dimensionality of this space with 25 features using UMAP for visualizing purposes. VoltRon is also capable of calculating the UMAP reduction from normalized data slots. Hence, we build a UMAP reduction from CLR data directly. However, UMAP will always be calculated from a PCA reduction by default (if a PCA embedding is found in the object). 

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

MBrain_Sec <- getUMAP(MBrain_Sec, data.type = "norm")

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

```

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

<br>

Using normalized cell type abundances, we can now generate k-nearest neighbor graphs and cluster the graph using leiden method.  

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

MBrain_Sec <- getProfileNeighbors(MBrain_Sec, data.type = "norm", method = "SNN")

vrGraphNames(MBrain_Sec)

```

```

[1] "SNN"

```

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

MBrain_Sec <- getClusters(MBrain_Sec, resolution = 0.6, graph = "SNN")

```

<br>

## Visualization

VoltRon incorporates distinct plotting functions for, e.g. embeddings, coordinates, heatmap and even barplots. We can now map the clusters we have generated on UMAP embeddings. 

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

# visualize 

g1 <- vrEmbeddingPlot(MBrain_Sec, embedding = "umap", group.by = "Sample")

g2 <- vrEmbeddingPlot(MBrain_Sec, embedding = "umap", group.by = "niche_clusters", label = TRUE)

g1 | g2

```

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

<br>

Mapping clusters on the spatial images and spots would show the niche structure across all four tissue sections.

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

vrSpatialPlot(MBrain_Sec, group.by = "niche_clusters", crop = TRUE, alpha = 1)

```

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

<br>

We use **vrHeatmapPlot** to investigate relative cell type abundances across these niche clusters. You will need to have **ComplexHeatmap** package in your namespace.

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

# install packages if necessary

if(!requireNamespace("ComplexHeatmap"))

  BiocManager::install("ComplexHeatmap")

# heatmap of niches

library(ComplexHeatmap)

vrHeatmapPlot(MBrain_Sec, features = vrFeatures(MBrain_Sec), group.by = "niche_clusters", 

              show_row_names = T, show_heatmap_legend = T)

```

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

<br>

# Cell-based Niche Clustering

Similar to spot-based spatial omics assays, we can build and cluster niche associated to each cell for spatial transcriptomics datasets in single cell resolution. For this, we require building niche assays for the collections of cells where a niche of cell is defined as a region of sets of regions with distinct cell type population that each of these cells belong to. 

Here, we dont require any scRNA reference dataset but we may first need to cluster and annotate cells in the RNA/transcriptome level profiles, and determine cell types. Then, we first detect the mixture of cell types within a spatial neighborhood around all cells and use that as a profile to perform clustering where these clusters will be associated with niches. 

## Import Xenium Data

For this, the data has to be already clustered (and annotated if possible). We will use the cluster labels generated at the end of the Xenium analysis workflow from [Cell/Spot Analysis](spotanalysis.html). You can download the VoltRon object with clustered and annotated Xenium cells along with the Visium assay from [here](https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/SpatialDataAlignment/Xenium_vs_Visium/VRBlock_data_clustered.rds). 

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

Xen_data <- readRDS("VRBlock_data_clustered.rds")

```

We will use all these 18 cell types used for annotating Xenium cells for detecting niches with distinct cellular type mixtures.

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

vrMainSpatial(Xen_data, assay = "Assay1") <- "main"

vrMainSpatial(Xen_data, assay = "Assay3") <- "main"

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/decon_xenium_clusters.png" class="center">

<br>

## Creating Niche Assay

For calculating niche profiles for each cell, we have to first build spatial neighborhoods around cells and capture the local cell type mixtures. Using **getSpatialNeighbors**, we build a spatial neighborhood graph to connect all cells to other cells within at most 15 distance apart. 

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

Xen_data <- getSpatialNeighbors(Xen_data, radius = 15, method = "radius")

vrGraphNames(Xen_data)

```

```

[1] "radius"

```

Now, we can build a niche assay for cells using the **getNicheAssay** function which will create an additional feature set for cells called **Niche**. Here, each cell type is a feature and the profile of a cell represents the relative abundance of cell types around each cell.

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

Xen_data <- getNicheAssay(Xen_data, label = "CellType", graph.type = "radius")

Xen_data

```

```

VoltRon Object 

10XBlock: 

  Layers: Section1 Section2 Section3 

Assays: Xenium(Main) Visium 

Features: RNA(Main) Niche

```

<br>

## Clustering

The Niche assay can be normalized similar to the spot-level niche analysis using **centred log ratio (CLR)** transformation. 

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

vrMainFeatureType(Xen_data) <- "Niche"

Xen_data <- normalizeData(Xen_data, method = "CLR")

```

Default clustering functions could be used to analyze the normalized niche profiles of cells to detect niches associated with each cell. However, we use K-means algorithm to perform the niche clustering. For this exercise, we pick an estimate of 7 clusters which will be the number of niche clusters we get.  

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

Xen_data <- getClusters(Xen_data, nclus = 7, method = "kmeans", label = "niche_clusters")

```

After the niche clustering, the metadata is updated and observed later like below.

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

head(Metadata(Xen_data))

```

<div><pre><code style="font-size: 10px;">               id Count assay_id  Assay    Layer   Sample CellType niche_clusters

1_Assay1 1_Assay1    28   Assay1 Xenium Section1 10XBlock   DCIS_1              2

2_Assay1 2_Assay1    94   Assay1 Xenium Section1 10XBlock   DCIS_2              2

3_Assay1 3_Assay1     9   Assay1 Xenium Section1 10XBlock   DCIS_1              2

4_Assay1 4_Assay1    11   Assay1 Xenium Section1 10XBlock   DCIS_1              2

5_Assay1 5_Assay1    48   Assay1 Xenium Section1 10XBlock   DCIS_2              2

6_Assay1 6_Assay1     7   Assay1 Xenium Section1 10XBlock   DCIS_1              2</code></pre></div>

<br>

## Visualization

After niche clustering, each cell in the Xenium assay will be assigned a niche which is initially a number which indicates the ID of each particular niche. It is up to the user to annotate, filter and visualize these niches moving forward. 

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

vrSpatialPlot(Xen_data, group.by = "niche_clusters", alpha = 1, legend.loc = "top")

```

We use **vrHeatmapPlot** to investigate the abundance of each cell type across the niche clusters. You will need to have **ComplexHeatmap** package in your namespace. We see that niche cluster 1 include all invasive tumor subtypes (IT 1-3). We see this for two subtypes of in situ ductal carcinoma (DCIS 1,2) subtypes as well other than a third DCIS subcluster being within proximity to myoepithelial cells. Niche cluster 6 also shows regions within the breast cancer tissue where T cells and B cells are found together abundantly. 

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

<br>

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

# install packages if necessary

if(!requireNamespace("ComplexHeatmap"))

  BiocManager::install("ComplexHeatmap")

# heatmap of niches

library(ComplexHeatmap)

vrHeatmapPlot(Xen_data, features = vrFeatures(Xen_data), group.by = "niche_clusters")

```

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

<br>