Interactive visualization of spatial omics neighborhoods

Background

The proliferation of spatial omic data has attracted attention from the modeling and visualization communities. Important themes that have emerged include the selection of spatially varying genes, derivation of spatial summary measures, and discovery of spatially consistent microenvironments. The resulting software packages allow analysts to generate overviews of spatial variation as well as focus on specific genomic features of interest.

Several studies propose feature-level models of spatial variation to select those with notable spatial expression patterns. SPARK fits a collection of random effects models with diverse sets of kernels to capture variation at several spatial scales Sun et al. (2020). Alternatively (Zhu and Sabatti, 2020), computed a measure of spatial variation based on a spatially induced graph laplacian; genes exhibiting similar patterns of spatial expression are then clustered. Alternatively (Hsu and Culhane, 2020), proposed an adaptation of Moran’s I-statistic to measure the extent of spatial clustering across cell types, highlighting the potential for the classical spatial statistics methods to support modern spatial omics analysis. Like NBFvis, these methods compute spatial statistics to summarize spatial omics datasets. However, they tend not to provide localized measures of spatial structure, focusing instead on tissue-level properties.

The Giotto package includes approaches to dimensionality reduction and interactive visualization of spatial omics data Dries et al. (2021b). Of particular interest, the package supports interactive visualization that dynamically links embeddings of expression measurements with corresponding cell locations. Note however that these embeddings are derived without reference to spatial features.

Similar to our approach, Spatial-LDA proposes a variation of the topic models that learns spatially consistent patterns of cell type mixing (Chen et al., 2020). This is achieved by tying together mixed memberships of neighboring cells in a structured prior, and the model is fitted using a custom optimization scheme. Regions with similar topic memberships can be interpreted as microenvironments. Our proposal has a similar data analytic goal; however, we aim to support more generic spatial features while preserving simplicity in implementation.

Simulation

We provided a toy simulation to clarify the differences in embeddings when neighborhood information is and is not used. We found that if only cell-level information is considered, the embeddings will be dominated by cell types and fail to reflect microenvironment structure.

Dataset construction

Assume that there is one tissue section with three cell types. For each cell, five proteins are measured. Different cell types have different protein profiles, which means the average measurements of proteins differ according to cell type. We assume that cell types are clustered spatially, but that these clusters are close enough so that some areas overlap. These overlapping areas can be considered different microenvironments since the local mixture of protein profiles is different from regions of pure cell types.

Figure 1 is the spatial plot for the simulated dataset. Two thousand “cells” are generated and divided into three different cell types according to a multinomial distribution with the probability (0.2, 0.3, 0.5) for Cell Type 1, 2, and 3,

Figure 1. Simulation of a tissue section with three different, partially overlapping cell types.

Overlapping regions can be thought of as their distinct microenvironments. The goal is to construct embeddings that reflect different mixing patterns.

For this demonstration, we imagined that protein abundances are drawn from a mixture of multivariate normals. The average of each mixture component represents the typical cell profile for each cell type. That is, for each cell, the measurements for each of the five proteins have the form,

Next, we simulated cell locations to get mixed spatial patterns. We used a different mixture of (now two-dimensional) multivariate normals. As before, component means center₁, center₂, and center₃ were drawn from a multivariate normal. Denoting the coordinates of cell i by $\left(x_i,y_i\right)$ and the spatial mean of cell type j by center_j, we drew,

After simulation, we obtained a 2000 × 5 expression matrix, each row of which corresponds to the simulated observation of one cell. We called this matrix the “single cell matrix”.

To extract neighborhood information for each cell, we first found neighborhoods with a given radius (here we used 0.2 units in length). We then calculated statistics within each neighborhood. We chose quantiles of protein content as neighborhood-based statistics, which are simple but effective. For every cell and protein, we derived 21 quantiles q₀, q_0.05, …, q₁ in the neighborhood. After calculation, an extended 2000 × 105 matrix was obtained. We called this the “neighborhood matrix”.

Comparison

We next applied dimensionality reduction methods to both the single cell and the neighborhood matrices, in order to clarify the difference in the resulting embeddings.

First, we applied uniform manifold approximation and projection (UMAP) McInnes et al. (2018) to the single cell matrix, whose low-dimensional embeddings are shown in Figure 2. Only three separate clusters are visible in the embedding plot, each corresponding to a cell type. These UMAP embeddings ignore the microenvironments of mixed cell types along cluster borders in the spatial plot. This result indicates that, when spatial information is not directly incorporated, the low-dimensional embeddings are dominated by cell types and fail to distinguish microenvironments.

Figure 2. A UMAP embedding plot on the single cell matrix.

Cell types clustered with one another, but different mixture patterns were not observed. The embeddings were dominated by the cell types, obscuring the presence of microenvironments.

In contrast, the embedding plot of the neighborhood matrix detects microenvironment structures; see Figure 3. We noticed that there were still three clusters consisting of pure cell types. However, there were additional clusters of mixed cell types. Between the three pure clusters is a region corresponding to microenvironments with mixed cell types in the spatial plot. Furthermore, we noticed that this region can be further divided into spatially consistent “subclusters”. For instance, one region with only blue and green cells was related to the blue-green spatial boundary. This can be treated as a unique microenvironment. Similar red-blue and red-green regions were also visible.

Figure 3. A UMAP embedding plot on the neighborhood matrix.

Though simple, the neighborhood quantile statistics make it possible to detect mixture microenvironments. We could further find subclusters like the red-green mixture in the central cluster.

In summary, UMAP embeddings using the single cell matrix were dominated by cell types and failed to detect microenvironments with mixed cell types. However, by simply applying UMAP to the neighborhood matrix, we were able to detect these spatially meaningful microenvironments.

Example

We next discuss a specific instantiation of this general procedure, describing the neighborhood and featurization choices used in the Results section and implemented in NBFvis. There, $N$ gives the number of cells in one tissue section, $D$ is the number of proteins measured, and $s$ is the spatial location matrix of the segmented cells. We applied a PCA to the expression matrix $X\in {\mathrm{ℝ}}^{N\times D}$ and then derived the reduced expression matrix $\widehat{X}\in {\mathrm{ℝ}}^{N\times P}$. Neighborhoods were constructed by keeping the K nearest neighbors that are also within a given radius.

We used two types of featurization functions $T_j$ – quantile features and network features. For the i^th cell’s neighborhood, Z quantiles $\left(q_1^{i,k},q_2^{i,k},\dots ,q_Z^{i,k}\right)$ were calculated for the k^th protein, where $k=1,2,\dots ,P$. It means that we derived a PZ-dimensional vector $\left(q_1^{i,1},q_2^{i,1},\dots ,q_{Z-1}^{i,P},q_Z^{i,P}\right)$ for each neighborhood. Thus, $T_{\text{quantile}}\left(X\right):{\mathrm{ℝ}}^{N\times P}\to {\mathrm{ℝ}}^{N\times \mathit{PZ}}$. After featurization, we obtained an $N\times PZ$ matrix, which we called the “quantile matrix”. Next, consider the construction of network features. Let $G_i$ denote the geometric graph associated with $m_i$, using the metric induced by $s$. Based on $G_i$, we can calculate a variety of node or edge features. The associated network featurization here is $T_{\text{network}}\left(X\right):{\mathrm{ℝ}}^{N\times P}\to {\mathrm{ℝ}}^{N\times M}$, where $M$ is the number of network statistics. For example, in the experiments below, we used the number of edges $\text{degree}\left(G_i\right)$ and a variety of centrality measures. We used an ensemble of 29 network-based statistics in our example, detailed in Table 1 and the Extended Data XTH1114 and Sankaran (2022).

The final featurization combined both quantile and network features,

$T\left(X\right)$ is a $N\times \left(\mathit{PZ},+,M\right)$ neighborhood matrix. Rescaling was applied to this neighborhood matrix so that every column was on a similar scale. This rescaled neighborhood matrix was passed to UMAP to obtain low-dimensional embeddings. These embeddings could then be clustered to identify distinct microenvironments. The whole workflow is shown in Figure 4.

Figure 4. The general workflow of the neighborhood-based featurization.

A) The spatial omics datasets are composed of two parts: spatial coordinates and expressions of each cell. We applied dimensionality reduction to the expression matrix to derive a reduced one on which features are derived. B) We treat each cell $x_i$ as a center and build a neighborhood for it. There are two general ways to construct a neighborhood: using distance between spatial coordinates and using the K-nearest neighbors. In our example, we combined these methods in the following way. First, we included cells whose distance from the central cell was less than a given length. Then we only kept the nearest K cells as its neighbors. C) To derive the quantile featurization for a gene, the quantiles of the distribution of each neighborhood’s expression for that gene are calculated. D) For the centrality featurization, we built a network within the neighborhood. Edges exist between two cells that are close enough. Then, we calculated centralities with respect to the center cell $x_i$. E) We concatenated and rescaled these derived features into what we call a neighborhood matrix. F) We applied uniform manifold approximation and projection (UMAP) to the neighborhood matrix to obtain embeddings for each cell, where we directly applied clustering algorithms. See the Example subsection below for details of the implementation.

Implementation details

Several subtle but important details are worth noting. Before we calculate a featurization matrix, a preliminary dimensionality reduction method is needed. First, applying dimensionality reduction decreases the computational burden of downstream analysis. Computing quantiles for each feature in a high-dimensional dataset further increases the dimensionality. For example, computing 10 quantiles for each of 100 variables results in 1000 columns, which significantly increases the computational burden of embedding. Second, a statistical reason for dimensionality reduction is to reduce the noise in the original high-dimensional dataset. If the original data are effectively low-rank, then dimensionality reduction method will reduce unnecessary noise while preserving most statistical information, which is beneficial for the following embedding.

Another detail is the rescaling of the neighborhood matrix. Although the neighborhood matrix could have hundreds or even thousands of columns, there is no need to apply a preliminary dimensionality reduction to it, since all values are approximately comparable. However, it is necessary to rescale the neighborhood matrix because the ranges of different statistics vary dramatically, causing one or two variables with large variance to dominate the whole UMAP embedding. For instance, the entries in the quantile matrix were between -1.5 and 1.5 in the TNBC dataset, but for the network matrix, it is common to have some network statistics larger than 10. These network statistics would dominate the UMAP embedding if no rescaling is applied.

Visualization design

We devised an interactive Shiny app (Chang et al., 2015) to analyze outputs from the neighborhood-based analysis, supporting visualization of microenvironment differences. In this subsection, we discuss the design and visual queries supported by the interface.

Figure 5 shows the first component of the Shiny app, the UMAP embedding and linked spatial plot. This is used to relate the low-dimensional embeddings of each cell’s neighborhood features to its overall spatial context. Figure 5a is the two-dimensional UMAP embedding plot derived from the neighborhood matrix. Each point corresponds to one cell. The closer these points are, the more similar their neighborhood featurizations are. To clearly visualize the distribution of cell types, the points in Figure 5a are colored according to cell types.

Figure 5. The first component: the UMAP embedding and spatial plots.

Part (a) is the two-dimensional embedding of the neighborhood matrix, and (b) is the original spatial layout of cell types.

Figure 5b is the spatial plot. Each point here represents a cell center, derived from the original cell polygon in the tissue section. As before, different cell types are distinguished by colors. Furthermore, the two panels in Figure 5 are dynamically linked. When points are selected in one plot through a mouse brush, the corresponding points will also be highlighted in the other plot. Figure 6 shows the highlighted points in these two plots after one such selection. We can click on the legend on the sidebar to deselect these cell types so that they do not appear. Figure 7 shows the embedding plot and scatter plot after deselecting the immune cells.

Figure 6. Area selection.

Figure 7. Cell types can be deselected by clicking on the legend.

The second component of the Shiny app shows the same embeddings but colored by K-means cluster rather than cell type. For example, in Figure 8, the positions of points are still the same as in Figure 5, but they are clustered into five K-means clusters. A slider is provided at the top of the second component in Figure 8, which is used for changing the value of K in the K-means clustering.

Figure 8. The second component: UMAP embedding plot and spatial plot of K-means clustering.

The third component of this Shiny app supports the comparison of expression levels across K-means clusters using a heatmap, structure plot, and histogram; see Figure 9. There are three tab panels with which we can switch between these three plots. Before making a further comparison, we can filter to cells of interest using the checkboxes at the top of Figure 9. Two groups of checkboxes are offered to select the cell types and K-means clusters to focus on. Based on the filtered cells, an expression heatmap of K-means clusters is provided in Figure 9. By default, it shows the top 10 most differentially expressed features across the selected clusters, based on the median of expression value in each cluster. A numeric input is offered above the heatmap – this controls the number of features appearing in the heatmap. The structure plot of the selected K-means clusters is provided in Figure 10, with which we can see the proportion of each cell type across every cluster. The histogram of expression is available to compare the selected feature’s expression across clusters. For example, Figure 11 is the histogram of the HLA Class 1 content in Clusters 1 and 4. Note that a selection input box is offered above the histogram to change the selected feature easily.

Figure 9. The third component: Heatmap, structure plot, and histogram.

These views help describe clusters identified by K-means.

Figure 10. Structure plot of K-means clusters.

The dominant cell types in each clusters are shown clearly.

Figure 11. Histogram of expressions among K-means clusters.

Users could select different expressions by the selection input box above the histogram.

To show the spatial distribution of a specific feature’s expression, another combination of the embedding and spatial plot is provided in Figure 12. The colors of the points in Figure 12 reflect Human Leukocyte Antigen (HLA) Class 1 content. This expression plot highlights spatial characteristics of the expression content. In this case, expression is elevated in immune cells, especially those closest to the tumor-immune boundary.

Figure 12. Expression plot of selected cells.

Instead of coloring by cell type or K-means cluster, each cell is shaded according to the selected genomic feature.

Cell-level approach

We also used the visualization tool to show the UMAP embedding and clustering results when directly applied to the original cell-level protein expression matrix. We used the same preprocessing as Keren et al. (2019). This serves as a reference point against which to compare the proposed neighborhood-based featurization.

Figure 21 gives the UMAP embedding and spatial plot using the cell-level approach. We found two clusters in the Figure 21a, one mainly made up of immune and one of tumor cells, respectively. The result was similar to the simulation, where UMAP embeddings were dominated by the differences between cell types and microenvironments were hardly distinguishable.

Figure 21. The UMAP embedding and spatial plots obtained without neighborhood features.

Cells are shaded by cell type. Compare with Figure 13.

Figure 22 shows the clustering results after K-means clustering with K = 5. The clustered microenvironments were mixed with each other; in particular, it was difficult to distinguish a tumor-immune boundary microenvironment. Figure 23 compares the clustered spatial plots based on the cell-level and neighborhood-based approaches. In Figure 23a, the cells in Region 1 were a mixture of three microenvironments derived from the cell-level approach. It was difficult to identify which microenvironment this region belonged to. Although Region 2 of Figure 23a was mainly composed of Cluster 3, there were cells from Clusters 4 and 5 distributed throughout. Though in principle it is possible to distinguish microenvironments based on particular mixture patterns across cell types, doing so requires much more effort than examining the neighborhood-based representation.

Figure 22. A version of Figure 21 where cells are shaded by K-means clusters found in the embedding on the left.

Sub-cell type variation in the embedding plot does not correspond to spatially meaningful microenvironments. Compare with Figure 15.

Figure 23. A direct comparison of the spatial plots from Figures 15 and 22.

Microenvironments with similar expression patterns (and stable cell type mixtures) are enclosed in black boxes. Microenvironments are more clearly visible when using neighborhood-based featurization.

Compared with the cell-level approach, the neighborhood-based featurization has a noticeably clearer clustering result. In Region 1 of Figure 23b, the cells in the boundary of tumor cells are spatially consistent according to their own cell types. Further, in Region 2 of Figure 23b, we observe a dominant microenvironment without needing to parse mixed patterns of cell types.

Overall, the neighborhood-based featurization provides representations with better spatial consistency, simplifying the discovery of microenvironments.

Package

We next summarize how to use NBFvis to implement the proposed workflow. We first loaded the packages and dataset we need. The dataset patient4 is a 6643 × 59 data frame of all cells in the tissue section of Patient 4 in the TNBC data (Keren et al., 2019). We added two columns named x_center and y_center, which are the coordinates of the calculated cell centers from the spatial raster data.

library (NBFvis)
library (dplyr)
data (patient4)

We selected 41 variables from dsDNA to HLA Class 1, most of which are proteins and cell type markers. The quantile_matrix function generates the quantile matrix from each cell’s neighborhood.

Quantiles_patient4 <- quantiles_matrix(
   data = patient4 %>% select (dsDNA:HLA_Class_1),
   coordinate = patient4 %>% select(x_center,y_center),
   index = patient4$index,
   NN = 40,
   distance = 60,
   min_percentile = 0.1,
   max_percentile = 0.9,
   quantile_number = 17,
   method = pca_)

The function network_matrix first builds the network inside the neighborhood and then calculates the corresponding network statistics using the argument given by fun. In this example, we used the function centralities, also exported by our package.

centrality_patient4 <- network_matrix(
   coordinate = patient4 %>% select (ends_with("_center")),
   index = patient4$index,
   radius = 60,
   NN = 40,
   edge = 30,
   fun = centralities,
   length_output = 29,
   name_output = NULL)

The scales of these two matrices are not the same, which means rescaling is needed. Here we removed Column, index and n_neighborhood in the quantile matrix so that all the columns left were quantile and network variables. Normalization and centering were applied to the centralities matrix so that they had a similar scale to the quantile matrix. We then combined the quantile matrix and the rescaled network matrix to construct an extended featurization matrix, which we called the neighborhood matrix earlier.

neighborhood_info_patient4 <- cbind(
   Quantiles_patient4 %>% select(-index, -n_neighbor),
   scale (centrality_patient4 %>% select(-index)))

The final step was to input the neighborhood matrix, the cell dataset patient4, and the names of the variable of interest in the function NBFvis. This returned an interactive Shiny app that was the source of figures in the Results section,

NBF_vis(
   matrix = neighborhood_info_patient4,
   origin_data = patient4,
   var_names = colnames (patient4)[17:57])

Underlying data

The TNBC dataset of Keren et al. (2019) can be downloaded from https://www.angelolab.com/mibi-data.

Extended data

Analysis code available from: https://github.com/XTH1114/NBFvis

Archived analysis code as at time of publication: DOI: 10.5281/zenodo.6639613

License: GNU General Public License