Spatially variable genes

name: spatial
channels:
  - defaults
  - conda-forge
dependencies:
  - conda-forge::python=3.12.12
  - conda-forge::scanpy=1.12
  - conda-forge::leidenalg
  - conda-forge::pytorch-cpu
  - conda-forge::pip
  - pip:
      - squidpy==1.8.1
      - SpaGCN==1.2.7
      - SpatialDE==1.1.3
      - tangram-sc==1.0.4
      - cell2location==0.1.5

Motivation¶

One main analysis step for single-cell data is to identify highly-variable genes (HVGs) and perform feature selection to reduce the dimensionality of the dataset. HVGs are genes which show significantly different expression profiles between cells or distinct groups. Methods designed for this task, however, neglect the spatial context of cells and can therefore not identify spatial variation. A gene might for example be highly variable, but not show a distinct spatial pattern and is therefore not spatially variable.

Spatial variation can be caused by differences in cell-type composition, overall functional dependencies or cell-cell communication events and help to understand the underlying tissue biology. Methods designed to identify spatially variable genes (SVGs) are designed to quantify whether a gene shows a significant spatial pattern by typically decomposing spatial and non-spatial variation in the dataset Walker et al., 2022.

Several methods have been proposed for this task with varying complexity and different assumptions. Currently there is no consensus on which method works best and how to define spatial variability in general. SpatialDE Svensson et al., 2018, SpatialDE2 Kats et al., 2021 and SPARKZhu et al., 2021 Sun et al., 2020 use spatial correlation testing, SepalAndersson & Lundeberg, 2021 leverages a Gaussian diffusion on spatial expression, scGCOZhang et al., 2022 utilizes a graph cut method and SpaGCN Hu et al., 2021 identifies SVGs based on spatial domains identified through a graph convolutional neural network.

In this notebook we provide a pedagogical example using Squidpy Palla et al., 2022 and its implementation of Moran’s I to find SVGs and subsequently an example workflow for SpatialDE.

Environment setup and data¶

We first load the respective packages needed in this tutorial and the dataset.

The dataset used in this tutorial consists of 1 tissue slides from 1 mouse and is provided by 10x Genomics Space Ranger 1.1.0. The dataset was pre-processed in Squidpy, which provides a loading function for this dataset.

Moran’s I score in Squidpy¶

One approach for the identification of spatially variable genes is the Moran’s I score, a measure of spatial autocorrelation (correlation of signal, such as gene expression, in observations close in space).

It is defined as: $I = \frac{n}{W}\frac{{\mathop {\sum }\nolimits_{i = 1}^n \mathop {\sum }\nolimits_{j = 1}^n w_{i,j}z_iz_j}}{{\mathop {\sum }\nolimits_{i = 1}^n z_i^2}}$ where

• $z_{i}$ is the deviation of the feature from the mean $\left( {x_i - \bar X} \right)$

• $w_{i,j}$ is the spatial weight between observations

• $n$ is the number of spatial units

• $W$ is the sum of all $w_{i,j}$

It can be computed with Squidpy with 1 line. For the purpose of the example, we will compute it only for a few genes.

The method adds a dataframe to adata.uns under the key moranI. We can inspect the result now:

The Squidpy implementation of Moran’s I computed for every gene

• I so the Moran’s I,

• pval_norm a p-value under normality assumption.

• var_norm the variance of the Moran’s I under normality assumption.

• {p_val}_{corr_method} the corrected p-values.

Let us look at two of the identified and significant genes, for example Nrgn and Ttr with corrected p-values of 0.0.

We can see that the expression of both of these genes seems to show a distinct localization in the tissue. It should be noted that they might (or might not) be also marker genes for specific cell clusters. One interpretation of spatially variable genes identification is that it is an orthogonal way to perform feature selection, by selecting genes that show a variability in space (instead of, for instance, across observations, as it is usually done).

SpatialDE¶

SpatialDE identifies spatially variable genes through Gaussian process regression. The spatial method decomposes each gene’s expression variability into a spatial and nonspatial component. It then computes the ratio between the spatial and nonspatial variance term to quantify the overall spatial variance present in the dataset. In order to identify significant spatially variable genes, SpatialDE compares the full model that has access to the spatial component to a model without this term.

We are using the same dataset as we did for computing Moran’s I. As SpatialDE requires the counts table saved as a DataFrame with unique variable names, we first make all variable names unique with the help of the respective Scanpy function.

Next, we collect the raw count table with all barcode names and variable names saved as dataframe indexes and columns. Scanpy provides an efficient function for this, get.obs_df which collects the respective keys stored in adata.

SpatialDE additionally requires the total counts and the spatial coordinates in the form of a DataFrame. We can use the same Scanpy function of collecting these items as before.

SpatialDE assumes a normally distributed noise. Since we just extracted raw counts, which empirically follow a negative binomial distribution, the counts data first needs to be transformed to a normal distributed noise. For this purpose SpatialDE uses a technique based on Anscombe:

The transformed data still might include a varying library size across the spatial samples which can introduce a bias in the gene expression. SpatialDE recommends to account for this before actually performing the spatial test and regressing it out with the provided function:

We can now run the actual spatial test by passing the spatial coordinates and the normalized counts to SpatialDE. On the dataset, we are using here SpatialDE takes roughly 15 minutes when running it on all genes.

We can now inspect the result:

The resulting DataFrame contains the following important columns:

• g, the gene name

• l, a parameter indicating the genes distance scale a gene changes expression over

• pval, the p-value for spatial differential expression

• qval, the corrected p-value after correcting for multiple testing

We can now sort the result based on the corrected p-values (qval) and to better read the DataFrame, subset the table to only show g, l and qval. We will additionally save the object as top10 to conveniently use it for downstream plotting.

We can now plot the top-three significant genes and inspect their spatial pattern. We additionally plot the clusters of the dataset to analyze whether the detected genes might be linked to specific clusters.

As we can observe, all three genes show spatial patterns. Esrra does not seem to be associated with a specific cluster, but shows a spatial pattern in primarily the cortex layers, the thalamus and hypothalamus. Fbxo31 and Jph3 are primarily expressed in the pyramidal layer.

Contributors¶

Authors¶

• Giovanni Palla

• Anna Schaar

Reviewers¶

• Lukas Heumos