Analysis and visualization of spatial transcriptomics data#

Author: Giovanni Palla

For up-to-date analysis tutorials, kindly check out SquidPy tutorials.

  • This tutorial demonstrates how to work with spatial transcriptomics data within Scanpy.

  • We focus on 10x Genomics Visium data, and provide an example for MERFISH.

[1]:

import scanpy as sc
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

[2]:

sc.logging.print_versions()
sc.set_figure_params(facecolor="white", figsize=(8, 8))
sc.settings.verbosity = 3

scanpy==1.5.0 anndata==0.7.1 umap==0.4.2 numpy==1.18.1 scipy==1.4.1 pandas==1.0.3 scikit-learn==0.22.1 statsmodels==0.11.0

Reading the data#

We will use a Visium spatial transcriptomics dataset of the human lymphnode, which is publicly available from the 10x genomics website: link.

The function datasets.visium_sge() downloads the dataset from 10x Genomics and returns an AnnData object that contains counts, images and spatial coordinates. We will calculate standards QC metrics with pp.calculate_qc_metrics and percentage of mitochondrial read counts per sample.

When using your own Visium data, use sc.read_visium() function to import it.

[3]:

adata = sc.datasets.visium_sge(sample_id="V1_Human_Lymph_Node")
adata.var_names_make_unique()
adata.var["mt"] = adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], inplace=True)

reading /Users/giovanni.palla/Projects/scanpy-tutorials/spatial/data/V1_Human_Lymph_Node/filtered_feature_bc_matrix.h5
 (0:00:01)

This is how the adata structure looks like for Visium data

[4]:

adata

[4]:

AnnData object with n_obs × n_vars = 4039 × 33538
    obs: 'in_tissue', 'array_row', 'array_col', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes', 'pct_counts_in_top_200_genes', 'pct_counts_in_top_500_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt'
    var: 'gene_ids', 'feature_types', 'genome', 'mt', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts'
    uns: 'spatial'
    obsm: 'spatial'

QC and preprocessing#

We perform some basic filtering of spots based on total counts and expressed genes

[5]:

fig, axs = plt.subplots(1, 4, figsize=(15, 4))
sns.distplot(adata.obs["total_counts"], kde=False, ax=axs[0])
sns.distplot(adata.obs["total_counts"][adata.obs["total_counts"] < 10000], kde=False, bins=40, ax=axs[1])
sns.distplot(adata.obs["n_genes_by_counts"], kde=False, bins=60, ax=axs[2])
sns.distplot(adata.obs["n_genes_by_counts"][adata.obs["n_genes_by_counts"] < 4000], kde=False, bins=60, ax=axs[3])

[5]:

<matplotlib.axes._subplots.AxesSubplot at 0x13045d0d0>
../_images/spatial_basic-analysis_12_1.png 
[6]:

sc.pp.filter_cells(adata, min_counts=5000)
sc.pp.filter_cells(adata, max_counts=35000)
adata = adata[adata.obs["pct_counts_mt"] < 20]
print(f"#cells after MT filter: {adata.n_obs}")
sc.pp.filter_genes(adata, min_cells=10)

filtered out 51 cells that have less than 5000 counts
filtered out 26 cells that have more than 35000 counts
#cells after MT filter: 3962
filtered out 15071 genes that are detected in less than 10 cells

We proceed to normalize Visium counts data with the built-in normalize_total method from Scanpy, and detect highly-variable genes (for later). Note that there are alternatives for normalization (see discussion in [Luecken19], and more recent alternatives such as SCTransform or GLM-PCA).

[7]:

sc.pp.normalize_total(adata, inplace=True)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, flavor="seurat", n_top_genes=2000)

normalizing counts per cell
    finished (0:00:00)
extracting highly variable genes
    finished (0:00:01)
--> added
    'highly_variable', boolean vector (adata.var)
    'means', float vector (adata.var)
    'dispersions', float vector (adata.var)
    'dispersions_norm', float vector (adata.var)

Manifold embedding and clustering based on transcriptional similarity#

To embed and cluster the manifold encoded by transcriptional similarity, we proceed as in the standard clustering tutorial.

[8]:

sc.pp.pca(adata)
sc.pp.neighbors(adata)
sc.tl.umap(adata)
sc.tl.leiden(adata, key_added="clusters")

computing PCA
    on highly variable genes
    with n_comps=50
    finished (0:00:01)
computing neighbors
    using 'X_pca' with n_pcs = 50
    finished: added to `.uns['neighbors']`
    `.obsp['distances']`, distances for each pair of neighbors
    `.obsp['connectivities']`, weighted adjacency matrix (0:00:02)
computing UMAP
    finished: added
    'X_umap', UMAP coordinates (adata.obsm) (0:00:06)
running Leiden clustering
    finished: found 10 clusters and added
    'clusters', the cluster labels (adata.obs, categorical) (0:00:01)

We plot some covariates to check if there is any particular structure in the UMAP associated with total counts and detected genes.

[9]:

plt.rcParams["figure.figsize"] = (4, 4)
sc.pl.umap(adata, color=["total_counts", "n_genes_by_counts", "clusters"], wspace=0.4)
../_images/spatial_basic-analysis_20_0.png

Visualization in spatial coordinates#

Let us now take a look at how total_counts and n_genes_by_counts behave in spatial coordinates. We will overlay the circular spots on top of the Hematoxylin and eosin stain (H&E) image provided, using the function sc.pl.spatial.

[10]:

plt.rcParams["figure.figsize"] = (8, 8)
sc.pl.spatial(adata, img_key="hires", color=["total_counts", "n_genes_by_counts"])
../_images/spatial_basic-analysis_23_0.png

The function sc.pl.spatial accepts 4 additional parameters:

  • img_key: key where the img is stored in the adata.uns element

  • crop_coord: coordinates to use for cropping (left, right, top, bottom)

  • alpha_img: alpha value for the transcparency of the image

  • bw: flag to convert the image into gray scale

Furthermore, in sc.pl.spatial, the size parameter changes its behaviour: it becomes a scaling factor for the spot sizes.

Before, we performed clustering in gene expression space, and visualized the results with UMAP. By visualizing clustered samples in spatial dimensions, we can gain insights into tissue organization and, potentially, into inter-cellular communication.

[11]:

sc.pl.spatial(adata, img_key="hires", color="clusters", size=1.5)
../_images/spatial_basic-analysis_25_0.png

Spots belonging to the same cluster in gene expression space often co-occur in spatial dimensions. For instance, spots belonging to cluster 5 are often surrounded by spots belonging to cluster 0.

We can zoom in specific regions of interests to gain qualitative insights. Furthermore, by changing the alpha values of the spots, we can visualize better the underlying tissue morphology from the H&E image.

[12]:

sc.pl.spatial(adata, img_key="hires", color="clusters", groups=["0", "5"], crop_coord=[1200, 1700, 1900, 1000], alpha=0.5, size=1.3)
../_images/spatial_basic-analysis_28_0.png

Cluster marker genes#

Let us further inspect cluster 5, which occurs in small groups of spots across the image.

Compute marker genes and plot a heatmap with expression levels of its top 10 marker genes across clusters.

[13]:

sc.tl.rank_genes_groups(adata, "clusters", method="t-test")
sc.pl.rank_genes_groups_heatmap(adata, groups="5", n_genes=10, groupby="clusters")

ranking genes
    finished: added to `.uns['rank_genes_groups']`
    'names', sorted np.recarray to be indexed by group ids
    'scores', sorted np.recarray to be indexed by group ids
    'logfoldchanges', sorted np.recarray to be indexed by group ids
    'pvals', sorted np.recarray to be indexed by group ids
    'pvals_adj', sorted np.recarray to be indexed by group ids (0:00:01)
WARNING: dendrogram data not found (using key=dendrogram_clusters). Running `sc.tl.dendrogram` with default parameters. For fine tuning it is recommended to run `sc.tl.dendrogram` independently.
    using 'X_pca' with n_pcs = 50
Storing dendrogram info using `.uns['dendrogram_clusters']`
WARNING: Groups are not reordered because the `groupby` categories and the `var_group_labels` are different.
categories: 0, 1, 2, etc.
var_group_labels: 5
../_images/spatial_basic-analysis_31_1.png

We see that CR2 recapitulates the spatial structure.

[14]:

sc.pl.spatial(adata, img_key="hires", color=["clusters", "CR2"])
../_images/spatial_basic-analysis_33_0.png

Spatially variable genes#

Spatial transcriptomics allows researchers to investigate how gene expression trends varies in space, thus identifying spatial patterns of gene expression. For this purpose, we use SpatialDE Svensson18 (code), a Gaussian process-based statistical framework that aims to identify spatially variable genes:

pip install spatialde

Recently, other tools have been proposed for the identification of spatially variable genes, such as:

First, we convert normalized counts and coordinates to pandas dataframe, needed for inputs to spatialDE.

Running SpatialDE takes considerable time.

[15]:

import SpatialDE

[16]:

%%time
counts = pd.DataFrame(adata.X.todense(), columns=adata.var_names, index=adata.obs_names)
coord = pd.DataFrame(adata.obsm['spatial'], columns=['x_coord', 'y_coord'], index=adata.obs_names)
results = SpatialDE.run(coord, counts)


CPU times: user 1h 22min 58s, sys: 1min 8s, total: 1h 24min 6s
Wall time: 21min 54s

We concatenate the results with the DataFrame of annotations of variables: adata.var.

[17]:

results.index = results["g"]
adata.var = pd.concat([adata.var, results.loc[adata.var.index.values, :]], axis=1)

Then we can inspect significant genes that varies in space and visualize them with sc.pl.spatial function.

[18]:

results.sort_values("qval").head(10)

[18]:
FSV M g l max_delta max_ll max_mu_hat max_s2_t_hat model n s2_FSV s2_logdelta time BIC max_ll_null LLR pval qval
g
EPHA1 0.092282 4 EPHA1 474.169746 9.649444 -912.478913 0.137565 0.009403 SE 3962 0.000010 0.001217 0.003795 1858.095843 -1002.530494 90.051581 0.0 0.0
CPVL 0.081288 4 CPVL 474.169746 11.087121 -2683.019176 0.398468 0.021507 SE 3962 0.000010 0.001476 0.007926 5399.176370 -2765.035269 82.016093 0.0 0.0
PPIA 0.298773 4 PPIA 474.169746 2.302420 -292.511593 3.092615 0.162267 SE 3962 0.000013 0.000332 0.003308 618.161202 -694.652548 402.140955 0.0 0.0
SAMD9L 0.125077 4 SAMD9L 474.169746 6.862115 -3003.385012 0.585090 0.041092 SE 3962 0.000009 0.000678 0.003048 6039.908042 -3202.809391 199.424379 0.0 0.0
COL1A2 0.090461 4 COL1A2 474.169746 9.863383 -3227.520356 0.740816 0.035660 SE 3962 0.000013 0.001612 0.002506 6488.178728 -3309.702658 82.182302 0.0 0.0
ARPC1B 0.120060 4 ARPC1B 474.169746 7.189890 -2285.547284 1.971545 0.074483 SE 3962 0.000013 0.001034 0.002559 4604.232584 -2412.439643 126.892359 0.0 0.0
ATP5MF 0.069258 4 ATP5MF 474.169746 13.183290 -2610.782278 1.671498 0.048891 SE 3962 0.000009 0.001708 0.009008 5254.702574 -2681.246440 70.464162 0.0 0.0
SYPL1 0.086861 4 SYPL1 474.169746 10.312907 -3329.098692 0.934305 0.039826 SE 3962 0.000012 0.001649 0.007649 6691.335402 -3401.604080 72.505388 0.0 0.0
PARP12 0.072946 4 PARP12 474.169746 12.467190 -3120.296512 0.678568 0.027263 SE 3962 0.000008 0.001420 0.008977 6273.731040 -3202.528124 82.231612 0.0 0.0
RARRES2 0.088121 4 RARRES2 474.169746 10.151330 -3289.541156 1.156461 0.045433 SE 3962 0.000010 0.001292 0.004011 6612.220329 -3381.437743 91.896587 0.0 0.0 
[19]:

sc.pl.spatial(adata, img_key="hires", color=["COL1A2", "SYPL1"], alpha=0.7)
../_images/spatial_basic-analysis_43_0.png

MERFISH example#

In case you have spatial data generated with FISH-based techniques, just read the cordinate table and assign it to the adata.obsm element.

Let’s take a look at the example from Xia et al. 2019.

First, we need to download the coordinate and counts data from the original publication.

[20]:

import urllib.request

[21]:

url_coord = "https://www.pnas.org/highwire/filestream/887973/field_highwire_adjunct_files/15/pnas.1912459116.sd15.xlsx"
filename_coord = "pnas.1912459116.sd15.xlsx"
urllib.request.urlretrieve(url_coord, filename_coord)

[21]:

('pnas.1912459116.sd15.xlsx', <http.client.HTTPMessage at 0x12e73a150>)

[22]:

url_counts = "https://www.pnas.org/highwire/filestream/887973/field_highwire_adjunct_files/12/pnas.1912459116.sd12.csv"
filename_counts = "pnas.1912459116.sd12.csv"
urllib.request.urlretrieve(url_counts, filename_counts)

[22]:

('pnas.1912459116.sd12.csv', <http.client.HTTPMessage at 0x12e74a4d0>)

And read the data in a AnnData object.

[23]:

coordinates = pd.read_excel("./pnas.1912459116.sd15.xlsx", index_col=0)
counts = sc.read_csv("./pnas.1912459116.sd12.csv").transpose()

[24]:

adata_merfish = counts[coordinates.index, :]
adata_merfish.obsm["spatial"] = coordinates.to_numpy()

We will perform standard preprocessing and dimensionality reduction.

[25]:

sc.pp.normalize_per_cell(adata_merfish, counts_per_cell_after=1e6)
sc.pp.log1p(adata_merfish)
sc.pp.pca(adata_merfish, n_comps=15)
sc.pp.neighbors(adata_merfish)
sc.tl.umap(adata_merfish)
sc.tl.leiden(adata_merfish, key_added="clusters", resolution=0.5)

normalizing by total count per cell
    finished (0:00:00): normalized adata.X and added    'n_counts', counts per cell before normalization (adata.obs)
computing PCA
    with n_comps=15
    finished (0:00:00)
computing neighbors
    using 'X_pca' with n_pcs = 15
    finished: added to `.uns['neighbors']`
    `.obsp['distances']`, distances for each pair of neighbors
    `.obsp['connectivities']`, weighted adjacency matrix (0:00:01)
computing UMAP
    finished: added
    'X_umap', UMAP coordinates (adata.obsm) (0:00:01)
running Leiden clustering
    finished: found 6 clusters and added
    'clusters', the cluster labels (adata.obs, categorical) (0:00:00)

The experiment consisted in measuring gene expression counts from a single cell type (cultured U2-OS cells). Clusters consist of cell states at different stages of the cell cycle. We don’t expect to see specific structure in spatial dimensions given the experimental setup.

We can visualize the clusters obtained from running Leiden in UMAP space and spatial coordinates like this.

[26]:

adata_merfish

[26]:

AnnData object with n_obs × n_vars = 645 × 12903
    obs: 'n_counts', 'clusters'
    uns: 'log1p', 'pca', 'neighbors', 'umap', 'leiden'
    obsm: 'spatial', 'X_pca', 'X_umap'
    varm: 'PCs'
    obsp: 'distances', 'connectivities'

[27]:

sc.pl.umap(adata_merfish, color="clusters")
sc.pl.embedding(adata_merfish, basis="spatial", color="clusters")
../_images/spatial_basic-analysis_56_0.png
../_images/spatial_basic-analysis_56_1.png

We hope you found the tutorial useful!

Report back to us which features/external tools you would like to see in Scanpy.

We are extending Scanpy and AnnData to support other spatial data types, such as Imaging Mass Cytometry and extend data structure to support spatial graphs and additional features. Stay tuned!