MISTy is designed to analyze spatial omics datasets within and between distinct spatial contexts referred to as views. This analysis can focus solely on structural information. Spatial transcriptomic methods such as Visium capture information from areas containing multiple cells. Then, deconvolution is applied to relate the measured data of the spots back to individual cells. In this vignette we will use the R package DOT
for deconvolution.
This vignette presents a workflow for the analysis of structural data, guiding users through the application of mistyR
to the results of DOT
deconvolution.
The package DOT
can be installed from Github remotes::install_github("saezlab/DOT")
.
Load the necessary packages:
# MISTy library(mistyR) library(future) # DOT library(DOT) # Loading experiment data library(Seurat) library(SeuratObject) # Data manipulation library(tidyverse) # Distances library(distances)
For this showcase, we use a 10X Visium spatial slide from Kuppe et al., 2022, where they created a spatial multi-omic map of human myocardial infarction. The tissue example data comes from the human heart of patient 14 which is in a later state after myocardial infarction. The Seurat object contains, among other things, the spot coordinates on the slides which we will need for decomposition. First, we have to download and extract the file:
# Download the data download.file("https://zenodo.org/records/6580069/files/10X_Visium_ACH005.tar.gz?download=1", destfile = "10X_Visium_ACH005.tar.gz", method = "curl") untar("10X_Visium_ACH005.tar.gz")
The next step is to load the data and extract the location of the spots. The rows are shifted, which means that the real distances between two spots are not always the same. It is therefore advantageous to use the pixel coordinates instead of row and column numbers, as the distances between these are represented accurately.
spatial_data <- readRDS("ACH005/ACH005.rds") geometry <- GetTissueCoordinates(spatial_data, cols = c("imagerow", "imagecol"), scale = NULL)
For deconvolution, we additionally need a reference single-cell data set containing a gene x cell count matrix and a vector containing the corresponding cell annotations. Kuppe et al., 2022, isolated nuclei from each sample's remaining tissue for snRNA-seq. The data corresponding to the same patient as the spatial data will be used as reference data in DOT
. First download the file:
download.file("https://www.dropbox.com/scl/fi/sq24xaavxplkc98iimvpz/hca_p14.rds?rlkey=h8cyxzhypavkydbv0z3pqadus&dl=1", destfile = "hca_p14.rds", mode = "wb")
Now load the data. From this, we retrieve a gene x cell count matrix and the respective cell annotations.
ref_data <- readRDS("hca_p14.rds") ref_counts_P14 <- ref_data$counts ref_ct <- ref_data$celltypes
Next, we need to set up the DOT object. The inputs we need are the count matrix and pixel coordinates of the spatial data and the count matrix and cell annotations of the single-cell reference data.
dot.srt <-setup.srt(srt_data = spatial_data@assays$Spatial@counts, srt_coords = geometry) dot.ref <- setup.ref(ref_data = ref_counts_P14, ref_annotations = ref_ct, 10) dot <- create.DOT(dot.srt, dot.ref)
Now we can carry out deconvolution:
# Run DOT dot <- run.DOT.lowresolution(dot)
The results can be found under dot@weights
. To obtain the calculated cell-type proportion per spot, we normalize the result to a row sum of 1.
# Normalize DOT results DOT_weights <- sweep(dot@weights, 1, rowSums(dot@weights), "/")
Now we can visually explore the slide itself and the abundance of cell types at each spot.
# Tissue Slide SpatialPlot(spatial_data, keep.scale = NULL, alpha = 0) # Results DOT draw_maps(geometry, DOT_weights, background = "white", normalize = FALSE, ncol = 3, viridis_option = "viridis")
Based on the plots, we can observe that some cell types are found more frequently than others. Additionally, we can identify patterns in the distribution of cells, with some being widespread across the entire slide while others are concentrated in specific areas. Furthermore, there are cell types that share a similar distribution.
First, we need to define an intraview that captures the cell type proportions within a spot. To capture the distribution of cell type proportions in the surrounding tissue, we add a paraview. For this vignette, the radius we choose is the mean of the distance to the nearest neighbor plus the standard deviation. We calculate the weights of each spot with family = gaussian. Then we run MISTy and collect the results.
# Calculating the radius geom_dist <- as.matrix(distances(geometry)) dist_nn <- apply(geom_dist, 1, function(x) (sort(x)[2])) paraview_radius <- ceiling(mean(dist_nn+ sd(dist_nn))) # Create views heart_views <- create_initial_view(as.data.frame(DOT_weights)) %>% add_paraview(geometry, l= paraview_radius, family = "gaussian") # Run misty and collect results run_misty(heart_views, "result/vignette_structural_pipeline") misty_results <- collect_results("result/vignette_structural_pipeline")
With the collected results, we can now answer the following questions:
Here we can look at two different statistics: multi.R2
shows the total variance explained by the multiview model. gain.R2
shows the increase in explainable variance from the paraview.
misty_results %>% plot_improvement_stats("multi.R2") %>% plot_improvement_stats("gain.R2")
The paraview particularly increases the explained variance for adipocytes. In general, the significant gain in R^2^ can be interpreted as the following:
"We can better explain the expression of marker X when we consider additional views other than the intrinsic view."
To explain the contributions, we can visualize the importance of each cell type in predicting the cell type distribution for each view separately. With trim
, we display only targets with a value above 50% for multi.R2
. To set an importance threshold we would apply cutoff
.
First, for the intrinsic view:
misty_results %>% plot_interaction_heatmap(view = "intra", clean = TRUE, trim.measure = "multi.R2", trim = 50)
We can observe that cardiomyocytes are a significant predictor for some cell types when in the same spot. To identify the target with the best prediction by cardiomyocytes, we can view the importance values as follows:
misty_results$importances.aggregated %>% filter(view == "intra", Predictor == "CM") %>% arrange(-Importance)
Let's take a look at the spatial distribution of Cardiomyocytes and their most important target, fibroblasts, in the tissue slide:
draw_maps(geometry, DOT_weights[, c("Fib", "CM")], background = "white", size = 1.25, normalize = FALSE, ncol = 1, viridis_option = "viridis")
We can observe that areas with high proportions of cardiomyocytes have low proportions of fibroblasts and vice versa.
Now we repeat this analysis with the paraview:
misty_results %>% plot_interaction_heatmap(view = "para.126", clean = TRUE, trim = 0.1, trim.measure = "gain.R2")
Here, we select the target adipocytes, as we know from previous analysis that adipocytes have the highest gain.R2
. The best predictor for adipocytes are Myeloid cells. To better identify the localization of the two cell types, we set the color scaling to a smaller range, as there are a few spots with a high proportion, which makes the distribution of spots with a low proportion difficult to recognize.
draw_maps(geometry, DOT_weights[, c("Myeloid","Adipo")], background = "white", size = 1.25, normalize = FALSE, ncol = 1, viridis_option = "viridis") + scale_colour_viridis_c(limits = c(0,0.33))
The plots show us that, in some places, the localization of the two cell types overlap.
browseVignettes("mistyR")
Here is the output of sessionInfo()
at the point when this document was compiled.
sessionInfo()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.