Niche Clustering

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 include a few number of cells that are likely from a combination of cell types within the tissue of origin. RNA deconvolution is then incorporated to estimate the percentage/abundance of these cell types within each spot/ROI given a reference scRNAseq dataset.

VoltRon includes wrapper commands for using popular RNA deconvolution methods such as RCTD (spot) and return estimated abundances as additional assays within each layer. These estimated percentages of cell types 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.

Import ST Data

For this tutorial we will analyze spot-based transcriptomic assays from Mouse Brain generated by the Visium instrument. The Mouse Brain Serial Section 1/2 datasets can be downloaded from here (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.

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) 

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

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

Import scRNA data

We will now import the scRNA data for reference which can be downloaded from here. Specifically, we will use a scRNA reference of Mouse cortical adult brain with 14,000 cells from the Allen Institute, generated with the SMART-Seq2 protocol. This scRNA data is also used by the Spatial Data Analysis tutorial in Seurat website.

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

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)

# 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

Spot Deconvolution with RCTD

In order to integrate the scRNA data and the Visium data sets within the VoltRon objects, we will use RCTD algorithm which is accessible with the spacexr package. For all layers with a Visium data, an additional assay within the same layer with **_decon** postfix will be created.

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

We can now switch to the Visium_decon assay where features are cell types from the scRNA reference and the data values are cell types percentages in each spot.

vrMainAssay(MBrain_Sec) <- "Visium_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 all other feature types.

vrSpatialFeaturePlot(MBrain_Sec, features = c("L4", "L5 PT", "Oligo", "Vip"), 
                     crop = TRUE, ncol = 2, alpha = 1, keep.scale = "all")

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 also be realized as niches, regions within tissue sections that have a distinct cell type composition compared to the remaining parts of the tissue.

We will now normalize and process relative cell type abundances to learn/detect these niche clusters. We treat cell type abundances as compositional data, hence we incorporate centred log ratio (CLR) transformation for normalizing them.

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 set without relying on a prior dimensionality reduction such as PCA.

VoltRon is also capable of calculating the UMAP reduction from normalized data slots. Hence, we build a UMAP reduction from CLR data. However, UMAP will always be calculated from a PCA reduction by default (if a PCA embedding is found in the object).

MBrain_Sec <- getUMAP(MBrain_Sec, assay = "Visium_decon", data.type = "norm")
vrEmbeddingPlot(MBrain_Sec, embedding = "umap", group.by = "Sample")

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

MBrain_Sec <- getProfileNeighbors(MBrain_Sec, assay = "Visium_decon", data.type = "norm", method = "SNN")
MBrain_Sec <- getClusters(MBrain_Sec, resolution = 0.6, graph = "SNN")

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.

# visualize 
g1 <- vrEmbeddingPlot(MBrain_Sec, embedding = "umap", group.by = "Sample")
g2 <- vrEmbeddingPlot(MBrain_Sec, embedding = "umap", group.by = "clusters", label = TRUE)
g1 | g2

Mapping clusters on the spatial images and spots would show the niche structure across the tissues and layers.

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

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

library(ComplexHeatmap)
vrHeatmapPlot(MBrain_Sec, features = vrFeatures(MBrain_Sec), group.by = "clusters", 
              show_row_names = T, show_heatmap_legend = T)