ELIXIR-SCO Spatial Omics Data Analysis 2026

Spatially-resolved transcriptomics data formats

Day 1 - Practical 0

6th October 2026

Jonas Windhager (NBIS)

Adapted from the MGC/BioSB Spatial Omics Analysis course 2025 by Benedetta Manzato (LUMC) with the help of Åsa Björklund (NBIS) and Claude Fable 5.

Agenda:
  1. Introduction to AnnData
  2. Introduction to SpatialData
  3. Accessing images
In [1]:
import warnings

warnings.filterwarnings("ignore")

# Importing necessary libraries for spatial data analysis and visualization
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scanpy as sc
import seaborn as sns

# Display the version of important packages for reproducibility
sc.logging.print_header()

# Set the default aesthetics for Scanpy plots
sc.set_figure_params(facecolor='white', figsize=(8, 8))  # White background, figure size 8x8
sc.settings.verbosity = 0  # Set verbosity to show only errors (0) # errors (0), warnings (1), info (2), hints (3)

AnnData format¶

Expression matrices are stored in the anndata h5ad format and need to be downloaded to a local file system for usage.

Key Components of an AnnData Object:

- .X: The main data matrix (usually cells by genes), where the rows are cells and the columns are genes. It typically stores raw or normalized gene expression data.

- .obs: Observations (metadata about the cells), such as cell type, experimental conditions, or quality control metrics.

- .var: Variables (metadata about the genes), such as gene names or other annotations.

- .obsm: Multi-dimensional annotations of cells, such as dimensionality reduction results (e.g., PCA, UMAP coordinates).

- .uns: Unstructured annotations, typically storing experimental metadata or plot settings.

- .layers: Stores multiple versions of the data matrix (e.g., raw counts, normalized data, etc.).

AnnData format structure
From Virshup et al. (2024), JOSS
  • Question: where do you expect spatial coordinates to be stored?

Dataset: Xenium In Situ of an Alzheimer's disease mouse brain (TgCRND8, 17.9 months)¶

This dataset is a coronal section of one brain hemisphere of an adult mouse, profiled with the 10x Genomics Xenium In Situ platform (imaging-based spatial transcriptomics). The sample (FFPE) comes from a well-characterized transgenic Alzheimer's disease model (TgCRND8) at the age of 17.9 months. A panel of 347 genes (the Xenium mouse brain panel with 248 genes plus a 99-gene custom add-on) was measured across ~62,000 segmented cells and ~63 million decoded transcripts.

The Xenium output bundle (Xenium_V1_FFPE_TgCRND8_17_9_months_outs) contains, among others:

  • experiment.xenium: experiment specifications (e.g. pixel size, panel).
  • cell_feature_matrix.h5: cell-by-gene count matrix (10x HDF5 format).
  • cells.parquet: cell metadata (centroids, areas, counts).
  • cell_boundaries.parquet / nucleus_boundaries.parquet: cell and nucleus boundary polygons.
  • transcripts.parquet: decoded transcripts with their coordinates.
  • morphology_mip.ome.tif / morphology_focus.ome.tif: DAPI (nuclear) images.
  • analysis.tar.gz: onboard secondary analysis (clustering, differential expression).

The same dataset will be used in the following practicals (cell segmentation, imaging-based data analysis). The data was retrieved from 10x Genomics and is available in the shared course data directory ../data/Xenium_V1_FFPE_TgCRND8_17_9_months_outs (populated by the course's download script).

For more details: 10x Genomics dataset page and the Xenium In Situ application note.

Xenium output file formats: https://www.10xgenomics.com/support/software/xenium-onboard-analysis/latest/analysis/xoa-output-understanding-outputs

Read data¶

In [2]:
xenium_dir = "../data/Xenium_V1_FFPE_TgCRND8_17_9_months_outs"

# the cell-by-gene count matrix is stored in the 10x HDF5 format, which scanpy can read directly
adata = sc.read_10x_h5(f"{xenium_dir}/cell_feature_matrix.h5")
adata
Out[2]:
AnnData object with n_obs × n_vars = 62268 × 347
    var: 'gene_ids', 'feature_types', 'genome'
    layers: None (.X)

The count matrix alone does not tell us where the cells are. The cell metadata (including the cell centroids) is stored in a separate Parquet file, which we read with pandas and attach to .obs (the row order matches the count matrix):

In [3]:
cells = pd.read_parquet(f"{xenium_dir}/cells.parquet")
cells["cell_id"] = cells["cell_id"].str.decode("utf-8")  # cell IDs are stored as bytes
cells = cells.set_index("cell_id")

adata.obs = adata.obs.join(cells)
adata
Out[3]:
AnnData object with n_obs × n_vars = 62268 × 347
    obs: 'x_centroid', 'y_centroid', 'transcript_counts', 'control_probe_counts', 'control_codeword_counts', 'unassigned_codeword_counts', 'total_counts', 'cell_area', 'nucleus_area'
    var: 'gene_ids', 'feature_types', 'genome'
    layers: None (.X)

Access cell metadata with .obs

In [4]:
adata.obs
Out[4]:
x_centroid y_centroid transcript_counts control_probe_counts control_codeword_counts unassigned_codeword_counts total_counts cell_area nucleus_area
aaabiggh-1 831.785336 755.803772 554 0 0 0 554 509.091563 65.386250
aaacfoel-1 821.733453 768.910446 291 0 0 0 291 275.498281 36.982969
aaaeefil-1 831.932053 780.367651 220 1 0 0 221 261.635312 15.849844
aaaehidd-1 853.614108 774.764157 1029 0 0 0 1029 818.547344 129.056563
aaagcbkg-1 821.639603 799.171515 453 0 0 1 454 522.232031 96.002188
... ... ... ... ... ... ... ... ... ...
oindmjog-1 5068.380664 4446.804883 35 0 0 0 35 229.709844 15.985313
oinecaba-1 5072.841504 4459.385132 124 0 0 0 124 359.353438 18.423750
oinehfmf-1 5251.196240 4366.193091 15 0 0 0 15 686.329844 24.745625
oinemiil-1 5248.493506 4382.190356 22 0 0 0 22 333.478906 20.681563
oingfoec-1 5272.023315 4387.369629 18 0 0 0 18 750.948438 9.798906

62268 rows × 9 columns

Access feature metadata with .var

In [5]:
adata.var
Out[5]:
gene_ids feature_types genome
2010300C02Rik ENSMUSG00000026090 Gene Expression Unknown
Abca7 ENSMUSG00000035722 Gene Expression Unknown
Acsbg1 ENSMUSG00000032281 Gene Expression Unknown
Acta2 ENSMUSG00000035783 Gene Expression Unknown
Acvrl1 ENSMUSG00000000530 Gene Expression Unknown
... ... ... ...
Vwc2l ENSMUSG00000045648 Gene Expression Unknown
Wfs1 ENSMUSG00000039474 Gene Expression Unknown
Zfp366 ENSMUSG00000050919 Gene Expression Unknown
Zfp536 ENSMUSG00000043456 Gene Expression Unknown
Zfpm2 ENSMUSG00000022306 Gene Expression Unknown

347 rows × 3 columns

In [6]:
print(f"Total number of cells/spots in this section: {adata.shape[0]}")
print(f"Size of the gene panel: {adata.shape[1]}")
Total number of cells/spots in this section: 62268
Size of the gene panel: 347

In .obs we see per-cell metadata such as transcript_counts, cell_area and nucleus_area, as well as the cell centroids x_centroid and y_centroid (in µm).

The Xenium onboard analysis also provides an unsupervised (graph-based) clustering of the cells (analysis.tar.gz). We add it to .obs as a categorical annotation (cells excluded from the onboard clustering get a missing value):

In [7]:
import tarfile

with tarfile.open(f"{xenium_dir}/analysis.tar.gz") as tar:
    clusters = pd.read_csv(
        tar.extractfile("analysis/clustering/gene_expression_graphclust/clusters.csv"),
        index_col="Barcode",
    )

adata.obs["cluster"] = pd.Categorical(
    clusters["Cluster"].reindex(adata.obs_names).astype("Int64").astype(str),
    categories=[str(c) for c in sorted(clusters["Cluster"].unique())],
)
adata.obs["cluster"].value_counts(dropna=False)
Out[7]:
cluster
1      6331
2      4663
3      4548
4      4359
5      4182
6      3528
7      3384
8      3051
9      2615
10     2558
11     2264
12     2120
13     1864
14     1766
15     1557
16     1521
17     1495
18     1495
19     1450
NaN    1439
20     1030
21      997
22      911
23      899
24      621
25      555
26      538
27      527
Name: count, dtype: int64

We can plot the annotations spatially in a scatterplot (using scanpy.pl.scatter or scanpy.pl.spatial).

Scanpy (and other spatial analysis tools) require adata.obsm['X_spatial'] or adata.obsm['spatial']. It is not yet present in our object (the coordinates only exist as columns of .obs), so we add it:

In [8]:
adata.obsm["spatial"] = adata.obs[["x_centroid", "y_centroid"]].to_numpy()
adata.obsm["spatial"]
Out[8]:
array([[ 831.7853363 ,  755.80377197],
       [ 821.73345337,  768.91044617],
       [ 831.93205261,  780.36765137],
       ...,
       [5251.19624023, 4366.19309082],
       [5248.49350586, 4382.19035645],
       [5272.02331543, 4387.36962891]], shape=(62268, 2))
In [9]:
annotations = ['cluster', 'transcript_counts']

for ann in annotations:
    fig, ax = plt.subplots(figsize=(8, 8))
    sc.pl.scatter(
        adata,
        x='x_centroid',
        y='y_centroid',
        color=ann,
        ax=ax,
        size=10,
        show=False,
        title=ann.replace('_', ' ').capitalize()
    )
    ax.invert_yaxis()
    ax.axis('equal')
    ax.set_xlabel('X Coordinate')
    ax.set_ylabel('Y Coordinate')
    ax.tick_params(axis='both', which='both',
                   bottom=False, left=False,
                   labelbottom=False, labelleft=False)
    ax.grid(False)
    plt.tight_layout()
    plt.show()
No description has been provided for this image
No description has been provided for this image

Alternatively we can plot a simple scatterplot with seaborn/matplotlib directly from adata.obs.

In [10]:
# Create the subplots
fig, axes = plt.subplots(1, 2, figsize=(15, 7.5))
fig.suptitle('Scatter Plots of Different Annotations', fontsize=16)

# Plot 'cluster' (categorical)
sns.scatterplot(ax=axes[0], x='x_centroid', y=-adata.obs['y_centroid'], hue='cluster', palette='tab20', data=adata.obs, s=4, legend=False)
axes[0].set_title('Cluster')
axes[0].set_xlabel('X Coordinate')
axes[0].set_ylabel('Y Coordinate')
axes[0].set_aspect('equal')
axes[0].tick_params(axis='both', which='both', bottom=False, left=False, labelbottom=False, labelleft=False)
axes[0].grid(False)  # Remove grid

# Plot 'transcript_counts' (continuous)
sns.scatterplot(ax=axes[1], x='x_centroid', y=-adata.obs['y_centroid'], hue='transcript_counts', palette='viridis', data=adata.obs, s=4, legend=False)
axes[1].set_title('Transcript counts')
axes[1].set_xlabel('X Coordinate')
axes[1].set_ylabel('Y Coordinate')
axes[1].set_aspect('equal')
axes[1].tick_params(axis='both', which='both', bottom=False, left=False, labelbottom=False, labelleft=False)
axes[1].grid(False)  # Remove grid

plt.tight_layout(rect=[0, 0, 1, 0.96])
plt.show()
No description has been provided for this image

You can also plot the gene expression of genes of interest:

In [11]:
genes_toplot = ['Slc17a7', 'Penk', 'Gfap']

sc.pl.spatial(adata, color=genes_toplot, spot_size=20)
No description has been provided for this image

Building blocks of anndata.AnnData class¶

Metadata¶

.obs and .var: These are stored as pandas.DataFrame objects to handle metadata efficiently.

Use cases include filtering, subsetting, and quick visualizations.

#### `pandas.DataFrame`

pandas.png

Let's filter the AnnData object to keep only cells assigned to cluster 1 or 2

In [12]:
adata_filtered_obs = adata[adata.obs['cluster'].isin(['1', '2'])]
adata_filtered_obs
Out[12]:
View of AnnData object with n_obs × n_vars = 10994 × 347
    obs: 'x_centroid', 'y_centroid', 'transcript_counts', 'control_probe_counts', 'control_codeword_counts', 'unassigned_codeword_counts', 'total_counts', 'cell_area', 'nucleus_area', 'cluster'
    var: 'gene_ids', 'feature_types', 'genome'
    uns: 'cluster_colors'
    obsm: 'spatial'
    layers: None (.X)
  • Exercise: Filter the AnnData object to keep only cells assigned to cluster '3'
In [13]:
adata_filtered_obs = adata[adata.obs['cluster'] == '3']
adata_filtered_obs.obs
Out[13]:
x_centroid y_centroid transcript_counts control_probe_counts control_codeword_counts unassigned_codeword_counts total_counts cell_area nucleus_area cluster
aabblhbk-1 831.906165 844.420578 105 0 0 0 105 126.663281 16.436875 3
aacmjcce-1 822.455258 856.651364 142 0 0 0 142 197.287656 7.450781 3
aadbolko-1 590.532520 865.743988 179 0 0 0 179 622.795000 20.591250 3
aaecjbfl-1 829.358240 368.995251 124 0 0 0 124 230.071094 29.125781 3
aaeecbmn-1 826.277924 381.758247 47 0 0 0 47 282.091094 24.971406 3
... ... ... ... ... ... ... ... ... ... ...
oibmklmg-1 4416.401611 311.657242 133 0 0 0 133 260.190313 24.248906 3
oieckffg-1 4907.441797 662.475339 196 0 0 0 196 257.526094 31.473906 3
oienoode-1 4924.513184 699.942886 149 0 0 0 149 195.932969 16.527188 3
oikhpgjh-1 4206.703540 552.846558 226 0 0 0 226 381.750938 30.119219 3
oingfoec-1 5272.023315 4387.369629 18 0 0 0 18 750.948438 9.798906 3

4548 rows × 10 columns

  • Exercise: Filter the AnnData object to keep only Gad1, Cst3, Apoe genes
  • Reminder: we can think of .obs as the rows of a pandas dataframe and .var as the columns
In [14]:
adata_filtered_var = adata[...]  # replace ... with your selection
adata_filtered_var.var
Out[14]:
gene_ids feature_types genome
2010300C02Rik ENSMUSG00000026090 Gene Expression Unknown
Abca7 ENSMUSG00000035722 Gene Expression Unknown
Acsbg1 ENSMUSG00000032281 Gene Expression Unknown
Acta2 ENSMUSG00000035783 Gene Expression Unknown
Acvrl1 ENSMUSG00000000530 Gene Expression Unknown
... ... ... ...
Vwc2l ENSMUSG00000045648 Gene Expression Unknown
Wfs1 ENSMUSG00000039474 Gene Expression Unknown
Zfp366 ENSMUSG00000050919 Gene Expression Unknown
Zfp536 ENSMUSG00000043456 Gene Expression Unknown
Zfpm2 ENSMUSG00000022306 Gene Expression Unknown

347 rows × 3 columns

scipy.sparse Matrices¶

scipy_sparse.png

Large datasets (e.g., single-cell or spatial datasets) are often sparse, meaning most of the matrix entries are zero. To save memory and computation time, these matrices are stored in sparse formats.

Common Formats:

scipy.sparse matrices:

  • csr_matrix (Compressed Sparse Row): Efficient for row slicing.
  • csc_matrix (Compressed Sparse Column): Efficient for column slicing.
  • coo_matrix (Coordinate): Useful for constructing sparse matrices.
  • Is the AnnData object currently in a dense or sparse format?
In [15]:
### check adata count matrix 
### hint: where is the count matrix stored?
In [16]:
from scipy import sparse

# Convert the count matrix to a sparse CSR matrix
# (a no-op here: the 10x HDF5 format already stores the counts sparsely)
sparse_X = sparse.csr_matrix(adata.X)
adata.X = sparse_X

NumPy Arrays¶

numpy.png

Dense data operations and algorithms that don’t support sparse matrices directly often require converting the sparse .X matrix to a NumPy array.

Convert sparse to dense matrix:

In [17]:
import scipy.sparse as sp

dense_array = adata.X.toarray() if sp.issparse(adata.X) else adata.X
dense_array
Out[17]:
array([[4., 2., 1., ..., 0., 0., 1.],
       [0., 0., 1., ..., 0., 0., 1.],
       [0., 0., 5., ..., 0., 0., 2.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]], shape=(62268, 347), dtype=float32)

Convert AnnData count matrix to pandas dataframe:

In [18]:
adata_df = adata.to_df()
adata_df
Out[18]:
2010300C02Rik Abca7 Acsbg1 Acta2 Acvrl1 Adamts2 Adamtsl1 Adgrl4 Aldh1a2 Aldh1l1 ... Unc13c Vat1l Vcan Vim Vip Vwc2l Wfs1 Zfp366 Zfp536 Zfpm2
aaabiggh-1 4.0 2.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 4.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 1.0
aaacfoel-1 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 2.0 ... 0.0 0.0 0.0 1.0 0.0 2.0 0.0 0.0 0.0 1.0
aaaeefil-1 0.0 0.0 5.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 2.0
aaaehidd-1 6.0 2.0 3.0 0.0 0.0 1.0 0.0 1.0 0.0 0.0 ... 0.0 8.0 0.0 1.0 0.0 4.0 4.0 0.0 0.0 2.0
aaagcbkg-1 2.0 3.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0 3.0 ... 0.0 2.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 2.0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
oindmjog-1 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 2.0 0.0 0.0 0.0 0.0 0.0 0.0
oinecaba-1 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 ... 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0
oinehfmf-1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
oinemiil-1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
oingfoec-1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0

62268 rows × 347 columns

Access 10X Genomics Visium data with sc.datasets.visium_sge¶

In [19]:
sc.settings.datasetdir = "../data"  # download into the shared course data directory

adata = sc.datasets.visium_sge(sample_id="V1_Mouse_Brain_Sagittal_Posterior")
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)
  0%|          | 0.00/9.26M [00:00<?, ?B/s]
  0%|          | 0.00/20.1M [00:00<?, ?B/s]
In [20]:
adata
Out[20]:
AnnData object with n_obs × n_vars = 3355 × 32285
    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'
    layers: None (.X)

Here images are stored in: adata.uns['spatial']['V1_Mouse_Brain_Sagittal_Posterior']['images']

In [21]:
adata.uns['spatial']['V1_Mouse_Brain_Sagittal_Posterior']['images'].keys()
Out[21]:
dict_keys(['hires', 'lowres'])
In [22]:
adata.uns['spatial']['V1_Mouse_Brain_Sagittal_Posterior']['images']['hires'][0].shape
Out[22]:
(2000, 3)

Given that the image is available, we can overlay spot‑level metadata on top of it to visualize spatial patterns.

In [23]:
plt.rcParams["figure.figsize"] = (8, 8)
sc.pl.spatial(adata, img_key="hires", color=["total_counts", "n_genes_by_counts"])
No description has been provided for this image

For detailed instructions on how to add tissue images to an AnnData object using Squidpy, see the tutorial here: Import spatial data in AnnData & Squidpy

Introduction to SpatialData objects¶

spatialdata.webp

Reading common raw data formats¶

The spatialdata-io library allows for reading common raw data formats. Here we show an example where 3 very different technologies from 10x Genomics are read into a general data structure, given by the SpatialData class. These are the elements that can be contained inside of a SpatialData object:

Screenshot 2025-01-07 at 3.58.39 PM.png

SpatialData Data Types¶

Data Type Recommended Object Type Description
Images xarray.DataArray or xarray.DataTree For single-scale or multi-scale images
Labels xarray.DataArray or xarray.DataTree Contain integer codes for different labels
Points dask.DataFrame Lazy version of pandas.DataFrame, used for point coordinates
Shapes geopandas.GeoDataFrame Store geometric objects such as polygons
Tables anndata.AnnData Used for tabular data with annotations

Before any discussion on file structure, let's see how we can download any dataset of interest from 10x (this can be used as a recipe for loading any similar dataset:

In [24]:
import matplotlib.pyplot as plt
import spatialdata as sd
import spatialdata_plot  # noqa: F401 (registers the `.pl` plotting accessor)
from spatialdata_io import visium_hd, xenium
from spatialdata.models import PointsModel
from spatialdata.transformations import get_transformation

Raw data can be downloaded directly from the vendor (or from a data repository), e.g. with curl/wget in a terminal. For this course, the datasets used below were already downloaded and extracted into the shared course data directory ../data by the course's download script, so the following commands are only shown for reference:

In [25]:
# Xenium (10x Genomics, ~6 GB):
# curl -fL -o Xenium_V1_FFPE_TgCRND8_17_9_months_outs.zip https://cf.10xgenomics.com/samples/xenium/1.4.0/Xenium_V1_FFPE_TgCRND8_17_9_months/Xenium_V1_FFPE_TgCRND8_17_9_months_outs.zip
# unzip Xenium_V1_FFPE_TgCRND8_17_9_months_outs.zip -d ../data/Xenium_V1_FFPE_TgCRND8_17_9_months_outs
In [26]:
# Visium HD (Zenodo, ~0.8 GB; 8 µm bins and segmented cells only):
# curl -fL -o data.zip "https://zenodo.org/records/22081261/files/data.zip?download=1"
# unzip data.zip -d ../data/data

Xenium¶

Let's start with the Xenium dataset from the beginning of this practical. The xenium reader of spatialdata-io takes the Xenium output directory and reads all its components into one SpatialData object: the cell-by-gene table, the morphology images, the cell/nucleus segmentation masks (labels) and boundaries (shapes), and the individual transcripts (points). Large components (images, transcripts) are loaded lazily.

In [27]:
sdata_xenium = xenium(xenium_dir)
sdata_xenium
Out[27]:
SpatialData object
├── Images
│     ├── 'morphology_focus': DataTree[cyx] (1, 24864, 31348), (1, 12432, 15674), (1, 6216, 7837), (1, 3108, 3918), (1, 1554, 1959)
│     └── 'morphology_mip': DataTree[cyx] (1, 24864, 31348), (1, 12432, 15674), (1, 6216, 7837), (1, 3108, 3918), (1, 1554, 1959)
├── Labels
│     ├── 'cell_labels': DataTree[yx] (24864, 31348), (12432, 15674), (6216, 7837), (3108, 3918), (1554, 1959)
│     └── 'nucleus_labels': DataTree[yx] (24864, 31348), (12432, 15674), (6216, 7837), (3108, 3918), (1554, 1959)
├── Points
│     └── 'transcripts': DataFrame with shape: (62751928, 10) (3D points)
├── Shapes
│     ├── 'cell_boundaries': GeoDataFrame shape: (62268, 1) (2D shapes)
│     └── 'nucleus_boundaries': GeoDataFrame shape: (62268, 1) (2D shapes)
└── Tables
      └── 'table': AnnData (62268, 347)
with coordinate systems:
    ▸ 'global', with elements:
        morphology_focus (Images), morphology_mip (Images), cell_labels (Labels), nucleus_labels (Labels), transcripts (Points), cell_boundaries (Shapes), nucleus_boundaries (Shapes)

The SpatialData Tables slot holds the same cell-by-gene matrix we assembled manually above (as an AnnData object), but here it is linked to the other elements, e.g. the cell labels via adata.obs['cell_labels']. We save it to disk in the h5ad format for later practicals:

In [28]:
adata = sdata_xenium["table"]
adata.write_h5ad("../data/Xenium_V1_FFPE_TgCRND8_17_9_months.h5ad")
adata
Out[28]:
AnnData object with n_obs × n_vars = 62268 × 347
    obs: 'cell_id', 'transcript_counts', 'control_probe_counts', 'control_codeword_counts', 'unassigned_codeword_counts', 'total_counts', 'cell_area', 'nucleus_area', 'region', 'z_level', 'cell_labels'
    var: 'gene_ids', 'feature_types', 'genome'
    uns: 'spatialdata_attrs'
    obsm: 'spatial'
    layers: None (.X)

Images (and labels) can have multiple scales, which are stored in a pyramid of downsampled representations. This allows plotting the whole section efficiently, even though the full-resolution image has ~780 million pixels:

In [29]:
sdata_xenium["morphology_focus"]
Out[29]:
<xarray.DataTree>
Group: /
├── Group: /scale0
│       Dimensions:  (c: 1, y: 24864, x: 31348)
│       Coordinates:
│         * c        (c) int64 8B 0
│         * y        (y) float64 199kB 0.5 1.5 2.5 3.5 ... 2.486e+04 2.486e+04 2.486e+04
│         * x        (x) float64 251kB 0.5 1.5 2.5 3.5 ... 3.135e+04 3.135e+04 3.135e+04
│       Data variables:
│           image    (c, y, x) uint16 2GB dask.array<chunksize=(1, 4096, 4096), meta=np.ndarray>
├── Group: /scale1
│       Dimensions:  (c: 1, y: 12432, x: 15674)
│       Coordinates:
│         * c        (c) int64 8B 0
│         * y        (y) float64 99kB 1.0 3.0 5.0 7.0 ... 2.486e+04 2.486e+04 2.486e+04
│         * x        (x) float64 125kB 1.0 3.0 5.0 7.0 ... 3.134e+04 3.134e+04 3.135e+04
│       Data variables:
│           image    (c, y, x) uint16 390MB dask.array<chunksize=(1, 4096, 4096), meta=np.ndarray>
├── Group: /scale2
│       Dimensions:  (c: 1, y: 6216, x: 7837)
│       Coordinates:
│         * c        (c) int64 8B 0
│         * y        (y) float64 50kB 2.0 6.0 10.0 ... 2.485e+04 2.486e+04 2.486e+04
│         * x        (x) float64 63kB 2.0 6.0 10.0 ... 3.134e+04 3.134e+04 3.135e+04
│       Data variables:
│           image    (c, y, x) uint16 97MB dask.array<chunksize=(1, 4096, 4096), meta=np.ndarray>
├── Group: /scale3
│       Dimensions:  (c: 1, y: 3108, x: 3918)
│       Coordinates:
│         * c        (c) int64 8B 0
│         * y        (y) float64 25kB 4.0 12.0 20.0 ... 2.484e+04 2.485e+04 2.486e+04
│         * x        (x) float64 31kB 4.001 12.0 20.0 ... 3.133e+04 3.134e+04 3.134e+04
│       Data variables:
│           image    (c, y, x) uint16 24MB dask.array<chunksize=(1, 3108, 3918), meta=np.ndarray>
└── Group: /scale4
        Dimensions:  (c: 1, y: 1554, x: 1959)
        Coordinates:
          * c        (c) int64 8B 0
          * y        (y) float64 12kB 8.0 24.0 40.0 ... 2.482e+04 2.484e+04 2.486e+04
          * x        (x) float64 16kB 8.001 24.0 40.01 ... 3.131e+04 3.132e+04 3.134e+04
        Data variables:
            image    (c, y, x) uint16 6MB dask.array<chunksize=(1, 1554, 1959), meta=np.ndarray>
xarray.DataTree
      • c: 1
      • y: 24864
      • x: 31348
      • c
        (c)
        int64
        0
        long_name :
        c
        units :
        array([0])
      • y
        (y)
        float64
        0.5 1.5 2.5 ... 2.486e+04 2.486e+04
        array([5.00000e-01, 1.50000e+00, 2.50000e+00, ..., 2.48615e+04, 2.48625e+04,2.48635e+04], shape=(24864,))
      • x
        (x)
        float64
        0.5 1.5 2.5 ... 3.135e+04 3.135e+04
        array([5.00000e-01, 1.50000e+00, 2.50000e+00, ..., 3.13455e+04, 3.13465e+04,3.13475e+04], shape=(31348,))
      • image
        (c, y, x)
        uint16
        dask.array<chunksize=(1, 4096, 4096), meta=np.ndarray>
        transform :
        {'global': Identity }
        Array Chunk
        Bytes 1.45 GiB 32.00 MiB
        Shape (1, 24864, 31348) (1, 4096, 4096)
        Dask graph 56 chunks in 3 graph layers
        Data type uint16 numpy.ndarray
        31348248641
      • c: 1
      • y: 12432
      • x: 15674
      • c
        (c)
        int64
        0
        long_name :
        c
        units :
        array([0])
      • y
        (y)
        float64
        1.0 3.0 5.0 ... 2.486e+04 2.486e+04
        array([1.0000e+00, 3.0000e+00, 5.0000e+00, ..., 2.4859e+04, 2.4861e+04,2.4863e+04], shape=(12432,))
      • x
        (x)
        float64
        1.0 3.0 5.0 ... 3.134e+04 3.135e+04
        array([1.0000e+00, 3.0000e+00, 5.0000e+00, ..., 3.1343e+04, 3.1345e+04,3.1347e+04], shape=(15674,))
      • image
        (c, y, x)
        uint16
        dask.array<chunksize=(1, 4096, 4096), meta=np.ndarray>
        transform :
        {'global': SequenceScale (y, x)[2. 2.]Identity }
        Array Chunk
        Bytes 371.66 MiB 32.00 MiB
        Shape (1, 12432, 15674) (1, 4096, 4096)
        Dask graph 16 chunks in 8 graph layers
        Data type uint16 numpy.ndarray
        15674124321
      • c: 1
      • y: 6216
      • x: 7837
      • c
        (c)
        int64
        0
        long_name :
        c
        units :
        array([0])
      • y
        (y)
        float64
        2.0 6.0 ... 2.486e+04 2.486e+04
        array([2.0000e+00, 6.0000e+00, 1.0000e+01, ..., 2.4854e+04, 2.4858e+04,2.4862e+04], shape=(6216,))
      • x
        (x)
        float64
        2.0 6.0 ... 3.134e+04 3.135e+04
        array([2.0000e+00, 6.0000e+00, 1.0000e+01, ..., 3.1338e+04, 3.1342e+04,3.1346e+04], shape=(7837,))
      • image
        (c, y, x)
        uint16
        dask.array<chunksize=(1, 4096, 4096), meta=np.ndarray>
        transform :
        {'global': SequenceScale (y, x)[4. 4.]Identity }
        Array Chunk
        Bytes 92.92 MiB 32.00 MiB
        Shape (1, 6216, 7837) (1, 4096, 4096)
        Dask graph 4 chunks in 13 graph layers
        Data type uint16 numpy.ndarray
        783762161
      • c: 1
      • y: 3108
      • x: 3918
      • c
        (c)
        int64
        0
        long_name :
        c
        units :
        array([0])
      • y
        (y)
        float64
        4.0 12.0 ... 2.485e+04 2.486e+04
        array([4.0000e+00, 1.2000e+01, 2.0000e+01, ..., 2.4844e+04, 2.4852e+04,2.4860e+04], shape=(3108,))
      • x
        (x)
        float64
        4.001 12.0 ... 3.134e+04 3.134e+04
        array([4.000510e+00, 1.200153e+01, 2.000255e+01, ..., 3.132800e+04,3.133600e+04, 3.134400e+04], shape=(3918,))
      • image
        (c, y, x)
        uint16
        dask.array<chunksize=(1, 3108, 3918), meta=np.ndarray>
        transform :
        {'global': SequenceScale (y, x)[8. 8.00102093]Identity }
        Array Chunk
        Bytes 23.23 MiB 23.23 MiB
        Shape (1, 3108, 3918) (1, 3108, 3918)
        Dask graph 1 chunks in 20 graph layers
        Data type uint16 numpy.ndarray
        391831081
      • c: 1
      • y: 1554
      • x: 1959
      • c
        (c)
        int64
        0
        long_name :
        c
        units :
        array([0])
      • y
        (y)
        float64
        8.0 24.0 ... 2.484e+04 2.486e+04
        array([8.0000e+00, 2.4000e+01, 4.0000e+01, ..., 2.4824e+04, 2.4840e+04,2.4856e+04], shape=(1554,))
      • x
        (x)
        float64
        8.001 24.0 ... 3.132e+04 3.134e+04
        array([8.001021e+00, 2.400306e+01, 4.000510e+01, ..., 3.130799e+04,3.132400e+04, 3.134000e+04], shape=(1959,))
      • image
        (c, y, x)
        uint16
        dask.array<chunksize=(1, 1554, 1959), meta=np.ndarray>
        transform :
        {'global': SequenceScale (y, x)[16. 16.00204186]Identity }
        Array Chunk
        Bytes 5.81 MiB 5.81 MiB
        Shape (1, 1554, 1959) (1, 1554, 1959)
        Dask graph 1 chunks in 24 graph layers
        Data type uint16 numpy.ndarray
        195915541
In [30]:
sdata_xenium.pl.render_images("morphology_focus").pl.show(title="Full morphology image")
No description has been provided for this image

Subsetting the section¶

For some tasks (e.g. cell segmentation in the next practical) the full section is too large for the compute resources available during the course. We therefore define a region of interest (ROI) in the upper-left part of the section (covering the cortex and hippocampus) that covers about 5% of the total area, i.e. each side of the ROI is sqrt(0.05) ≈ 22% of the full width/height. The ROI is defined deterministically (no interactive selection), so everyone works on the same crop. Note that the global coordinate system of the Xenium reader is in pixels of the morphology image (1 pixel = 0.2125 µm).

In [31]:
crop_fraction = 0.05  # fraction of the total area to keep
roi_center = (0.30, 0.25)  # center of the ROI as fraction of the full (width, height)

# extent of the full section (in the 'global' coordinate system)
extent = sd.get_extent(sdata_xenium["morphology_focus"], coordinate_system="global")
x_min, x_max = extent["x"]
y_min, y_max = extent["y"]

# size and position of an ROI covering `crop_fraction` of the area, centered at `roi_center`
roi_width = int((x_max - x_min) * np.sqrt(crop_fraction))
roi_height = int((y_max - y_min) * np.sqrt(crop_fraction))
x_start = int(x_min + (x_max - x_min) * roi_center[0] - roi_width / 2)
y_start = int(y_min + (y_max - y_min) * roi_center[1] - roi_height / 2)
x_end = x_start + roi_width
y_end = y_start + roi_height

print(f"Full section: {int(x_max - x_min)} x {int(y_max - y_min)} pixels")
print(f"ROI: x = {x_start} - {x_end}, y = {y_start} - {y_end} ({roi_width} x {roi_height} pixels)")
Full section: 31348 x 24864 pixels
ROI: x = 5899 - 12908, y = 3436 - 8995 (7009 x 5559 pixels)
In [32]:
from matplotlib import patches

fig, ax = plt.subplots(figsize=(8, 8))
sdata_xenium.pl.render_images("morphology_focus").pl.show(ax=ax, title="Full morphology image with ROI")
ax.add_patch(patches.Rectangle((x_start, y_start), roi_width, roi_height, linewidth=2, edgecolor="red", facecolor="none"))
plt.show()
No description has been provided for this image

A bounding box query crops all elements of the SpatialData object to the ROI at once: images and labels are cropped, shapes and points are filtered, and the table is subset to the remaining cells.

In [33]:
sdata_xenium_crop = sdata_xenium.query.bounding_box(
    axes=["x", "y"],
    min_coordinate=[x_start, y_start],
    max_coordinate=[x_end, y_end],
    target_coordinate_system="global",
)
sdata_xenium_crop
Out[33]:
SpatialData object
├── Images
│     ├── 'morphology_focus': DataTree[cyx] (1, 5559, 7009), (1, 2780, 3505), (1, 1390, 1752), (1, 695, 876), (1, 347, 438)
│     └── 'morphology_mip': DataTree[cyx] (1, 5559, 7009), (1, 2780, 3505), (1, 1390, 1752), (1, 695, 876), (1, 347, 438)
├── Labels
│     ├── 'cell_labels': DataTree[yx] (5559, 7009), (2780, 3505), (1390, 1752), (695, 876), (347, 438)
│     └── 'nucleus_labels': DataTree[yx] (5559, 7009), (2780, 3505), (1390, 1752), (695, 876), (347, 438)
├── Points
│     └── 'transcripts': DataFrame with shape: (<Delayed>, 10) (3D points)
├── Shapes
│     ├── 'cell_boundaries': GeoDataFrame shape: (5569, 1) (2D shapes)
│     └── 'nucleus_boundaries': GeoDataFrame shape: (5422, 1) (2D shapes)
└── Tables
      └── 'table': AnnData (5568, 347)
with coordinate systems:
    ▸ 'global', with elements:
        morphology_focus (Images), morphology_mip (Images), cell_labels (Labels), nucleus_labels (Labels), transcripts (Points), cell_boundaries (Shapes), nucleus_boundaries (Shapes)

Before saving, we make two changes so that the cropped data also opens in TissUUmaps:

  • The transcripts are a dask DataFrame split into several partitions, each written to its own Parquet file (part.0.parquet, part.1.parquet, ...). TissUUmaps reads only the first file, so we merge the partitions into one.
  • The count matrix X is stored row by row (CSR, one row per cell). TissUUmaps reads the expression of one gene at a time, which is much faster when the matrix is stored column by column (CSC, one column per gene).
In [34]:
# merge the transcript partitions, so they are written to a single Parquet file
transcripts = sdata_xenium_crop.points["transcripts"]
sdata_xenium_crop.points["transcripts"] = PointsModel.parse(
    transcripts.compute(),
    feature_key=transcripts.attrs["spatialdata_attrs"]["feature_key"],
    transformations=get_transformation(transcripts, get_all=True),
)

# store the count matrix by column (one gene at a time) instead of by row
table = sdata_xenium_crop.tables["table"]
table.X = sp.csc_matrix(table.X)

We save the cropped object to disk in the SpatialData Zarr format for the following practicals (this can take a few minutes, mostly for filtering the ~63 million transcripts):

In [ ]:
sdata_xenium_crop.write("../data/Xenium_V1_FFPE_TgCRND8_17_9_months_cropped.zarr", overwrite=True)

Let's have a look at the cropped morphology image, the cell segmentation and the transcripts of a single gene:

In [ ]:
fig, axes = plt.subplots(1, 2, figsize=(16, 8))
sdata_xenium_crop.pl.render_images("morphology_focus").pl.show(ax=axes[0], title="Cropped morphology image", colorbar=False)
sdata_xenium_crop.pl.render_images("morphology_focus").pl.render_shapes(
    "cell_boundaries", fill_alpha=0, outline_alpha=1, outline_color="red", outline_width=0.5
).pl.show(ax=axes[1], title="Cell boundaries", colorbar=False)
In [ ]:
fig, axes = plt.subplots(1, 2, figsize=(16, 8))
sdata_xenium_crop.pl.render_labels("cell_labels").pl.show(ax=axes[0], title="Cell labels")
sdata_xenium_crop.pl.render_labels("nucleus_labels").pl.render_points(
    "transcripts", color="feature_name", groups="Pvalb", palette="orange"
).pl.show(ax=axes[1], title="Nucleus labels and Pvalb transcripts")

Visium HD¶

Next, we read a sequencing-based dataset: a Visium HD 3' mouse brain section (10x Genomics). To save space, only the 8 µm bins and the segmented cells (as provided by Space Ranger) are included in the course data (Zenodo record 22081261). The visium_hd reader works on the Space Ranger output directory and detects the available bin sizes, images and segmentations.

In [38]:
sdata_visium_hd = visium_hd("../data/data", bin_size=8)
In [39]:
sdata_visium_hd
Out[39]:
SpatialData object
├── Images
│     ├── 'Visium_HD_3prime_Mouse_Brain_hires_image': DataArray[cyx] (3, 5492, 6000)
│     └── 'Visium_HD_3prime_Mouse_Brain_lowres_image': DataArray[cyx] (3, 549, 600)
├── Shapes
│     ├── 'Visium_HD_3prime_Mouse_Brain_cell_segmentations': GeoDataFrame shape: (84031, 2) (2D shapes)
│     └── 'Visium_HD_3prime_Mouse_Brain_square_008um': GeoDataFrame shape: (376419, 1) (2D shapes)
└── Tables
      ├── 'cell_segmentations': AnnData (84031, 33696)
      └── 'square_008um': AnnData (376419, 33696)
with coordinate systems:
    ▸ 'Visium_HD_3prime_Mouse_Brain', with elements:
        Visium_HD_3prime_Mouse_Brain_hires_image (Images), Visium_HD_3prime_Mouse_Brain_lowres_image (Images), Visium_HD_3prime_Mouse_Brain_cell_segmentations (Shapes), Visium_HD_3prime_Mouse_Brain_square_008um (Shapes)
    ▸ 'Visium_HD_3prime_Mouse_Brain_downscaled_hires', with elements:
        Visium_HD_3prime_Mouse_Brain_hires_image (Images), Visium_HD_3prime_Mouse_Brain_cell_segmentations (Shapes), Visium_HD_3prime_Mouse_Brain_square_008um (Shapes)
    ▸ 'Visium_HD_3prime_Mouse_Brain_downscaled_lowres', with elements:
        Visium_HD_3prime_Mouse_Brain_lowres_image (Images), Visium_HD_3prime_Mouse_Brain_cell_segmentations (Shapes), Visium_HD_3prime_Mouse_Brain_square_008um (Shapes)

Similarly to the regular Visium data, we can read and inspect the individual modalities contained in the data. We can see that the Visium HD data contains a lot more elements in the Shapes slot (376,419 8 µm bins vs. 3,355 spots) due to the increased resolution. In addition to the bins, the reader picked up the segmented cells (cell_segmentations), which come with their own table.

We are making unique all tables for all resolutions

In [40]:
for table_name, adata in sdata_visium_hd.tables.items():
    print(f"Processing table: {table_name}")
    adata.var_names_make_unique()
    
    
for table_name, adata in sdata_visium_hd.tables.items():
    duplicates = adata.var.index[adata.var.index.duplicated()]
    if duplicates.empty:
        print(f"All variable names are unique in {table_name}.")
    else:
        print(f"Duplicate variable names still exist in {table_name}: {duplicates}")
Processing table: square_008um
Processing table: cell_segmentations
All variable names are unique in square_008um.
All variable names are unique in cell_segmentations.

SpatialData Tables is simply an AnnData object

In [41]:
adata = sdata_visium_hd["square_008um"]
adata
Out[41]:
AnnData object with n_obs × n_vars = 376419 × 33696
    obs: 'in_tissue', 'array_row', 'array_col', 'location_id', 'region'
    var: 'gene_ids', 'feature_types', 'genome'
    uns: 'spatialdata_attrs'
    obsm: 'spatial'
    layers: None (.X)

SpatialData Images

Images can have multiple scales, which are stored in a pyramid for downsampled representations useful for more efficient visualization and analysis.

In [42]:
image_h = sdata_visium_hd["Visium_HD_3prime_Mouse_Brain_hires_image"]
image_h
Out[42]:
<xarray.DataArray 'image' (c: 3, y: 5492, x: 6000)> Size: 99MB
dask.array<transpose, shape=(3, 5492, 6000), dtype=uint8, chunksize=(3, 5492, 6000), chunktype=numpy.ndarray>
Coordinates:
  * c        (c) <U1 12B 'r' 'g' 'b'
  * y        (y) float64 44kB 0.5 1.5 2.5 3.5 ... 5.49e+03 5.49e+03 5.492e+03
  * x        (x) float64 48kB 0.5 1.5 2.5 3.5 ... 5.998e+03 5.998e+03 6e+03
Attributes:
    transform:  {'Visium_HD_3prime_Mouse_Brain_downscaled_hires': Identity , ...
xarray.DataArray
'image'
  • c: 3
  • y: 5492
  • x: 6000
  • dask.array<chunksize=(3, 5492, 6000), meta=np.ndarray>
    Array Chunk
    Bytes 94.28 MiB 94.28 MiB
    Shape (3, 5492, 6000) (3, 5492, 6000)
    Dask graph 1 chunks in 4 graph layers
    Data type uint8 numpy.ndarray
    6000 5492 3
    • c
      (c)
      <U1
      'r' 'g' 'b'
      long_name :
      c
      units :
      array(['r', 'g', 'b'], dtype='<U1')
    • y
      (y)
      float64
      0.5 1.5 2.5 ... 5.49e+03 5.492e+03
      array([5.0000e-01, 1.5000e+00, 2.5000e+00, ..., 5.4895e+03, 5.4905e+03,
             5.4915e+03], shape=(5492,))
    • x
      (x)
      float64
      0.5 1.5 2.5 ... 5.998e+03 6e+03
      array([5.0000e-01, 1.5000e+00, 2.5000e+00, ..., 5.9975e+03, 5.9985e+03,
             5.9995e+03], shape=(6000,))
  • transform :
    {'Visium_HD_3prime_Mouse_Brain_downscaled_hires': Identity , 'Visium_HD_3prime_Mouse_Brain': Scale (x, y) [3.55299992 3.55299992]}
In [43]:
image_l = sdata_visium_hd["Visium_HD_3prime_Mouse_Brain_lowres_image"]
image_l
Out[43]:
<xarray.DataArray 'image' (c: 3, y: 549, x: 600)> Size: 988kB
dask.array<transpose, shape=(3, 549, 600), dtype=uint8, chunksize=(3, 549, 600), chunktype=numpy.ndarray>
Coordinates:
  * c        (c) <U1 12B 'r' 'g' 'b'
  * y        (y) float64 4kB 0.5 1.5 2.5 3.5 4.5 ... 545.5 546.5 547.5 548.5
  * x        (x) float64 5kB 0.5 1.5 2.5 3.5 4.5 ... 596.5 597.5 598.5 599.5
Attributes:
    transform:  {'Visium_HD_3prime_Mouse_Brain_downscaled_lowres': Identity ,...
xarray.DataArray
'image'
  • c: 3
  • y: 549
  • x: 600
  • dask.array<chunksize=(3, 549, 600), meta=np.ndarray>
    Array Chunk
    Bytes 0.94 MiB 0.94 MiB
    Shape (3, 549, 600) (3, 549, 600)
    Dask graph 1 chunks in 4 graph layers
    Data type uint8 numpy.ndarray
    600 549 3
    • c
      (c)
      <U1
      'r' 'g' 'b'
      long_name :
      c
      units :
      array(['r', 'g', 'b'], dtype='<U1')
    • y
      (y)
      float64
      0.5 1.5 2.5 ... 546.5 547.5 548.5
      array([5.000e-01, 1.500e+00, 2.500e+00, ..., 5.465e+02, 5.475e+02, 5.485e+02],
            shape=(549,))
    • x
      (x)
      float64
      0.5 1.5 2.5 ... 597.5 598.5 599.5
      array([5.000e-01, 1.500e+00, 2.500e+00, ..., 5.975e+02, 5.985e+02, 5.995e+02],
            shape=(600,))
  • transform :
    {'Visium_HD_3prime_Mouse_Brain_downscaled_lowres': Identity , 'Visium_HD_3prime_Mouse_Brain': Scale (x, y) [35.52999922 35.52999922]}
In [44]:
sdata_visium_hd.pl.render_images().pl.show("Visium_HD_3prime_Mouse_Brain_downscaled_hires")
No description has been provided for this image
In [45]:
print(sdata_visium_hd.coordinate_systems)
['Visium_HD_3prime_Mouse_Brain_downscaled_hires', 'Visium_HD_3prime_Mouse_Brain_downscaled_lowres', 'Visium_HD_3prime_Mouse_Brain']

Geopandas¶

GeoPandas spatial data frame
From Introduction to GeoPandas

A GeoDataFrame is an object that inherits from a pandas.DataFrame, and as such it supports the APIs that pandas provides, but also new powerful APIs for geometric manipulation.

The core data structure in GeoPandas is the geopandas.GeoDataFrame, that can store geometry columns and perform spatial operations. The geopandas.GeoSeries handles the geometries. Therefore, your GeoDataFrame is a combination of pandas.Series, with traditional data (numerical, boolean, text etc.), and geopandas.GeoSeries, with geometries (points, polygons etc.).

The geometry column is what contains the spatial information.

In SpatialData, only dataframes where all the geometries are polygons/multipolygons are supported, either where all the geometries are shapely.Point and a column called radius is present. One can add new regular columns to the dataframe.

In [46]:
sdata_visium_hd['Visium_HD_3prime_Mouse_Brain_square_008um']
Out[46]:
geometry
location_id
0 POLYGON ((6185.027 7974.939, 6185.027 7945.728...
1 POLYGON ((3773.779 16118.174, 3773.779 16088.9...
2 POLYGON ((14157.852 7641.397, 14157.852 7612.1...
3 POLYGON ((4094.873 12784.773, 4094.873 12755.5...
4 POLYGON ((19208.134 3886.285, 19208.134 3857.0...
... ...
376414 POLYGON ((4216.623 16640.133, 4216.623 16610.9...
376415 POLYGON ((18752.337 18410.528, 18752.337 18381...
376416 POLYGON ((8438.564 5121.085, 8438.564 5091.873...
376417 POLYGON ((13278.098 7269.386, 13278.098 7240.1...
376418 POLYGON ((2570.9 12243.171, 2570.9 12213.959, ...

376419 rows × 1 columns

In [47]:
sdata_visium_hd.pl.render_shapes("Visium_HD_3prime_Mouse_Brain_square_008um").pl.show()
No description has been provided for this image

Read more: https://spatialdata.scverse.org/en/latest/tutorials/notebooks/notebooks/examples/models1.html and https://geopandas.org/en/stable/getting_started/introduction.html

In [ ]: