Note

This page was generated from 5.1 3D reconstruction on Drosophila bin 50 data with pairwise alignment.ipynb. Interactive online version: Colab badge. Some tutorial content may look better in light mode.

5.1 3D Reconstruction with Pairwise Alignment on Drosophila Data (bin50)

In this tutorial, you will learn how to reconstruct a 3D spatial transcriptomics by aligning the consecutive spatial transcriptomics slices serially (basic).

[2]:
import os
os.environ['CUDA_VISIBLE_DEVICES'] = '1'

import torch
device = 'cuda' if torch.cuda.is_available() else 'cpu'
print("Running this notebook on: ", device)

import spateo as st
print("Last run with spateo version:", st.__version__)

# Other imports
import warnings
warnings.filterwarnings('ignore')
import matplotlib.pyplot as plt
import dynamo as dyn
import numpy as np
import anndata as ad

# Uncomment the following if running on the server
import pyvista as pv
pv.start_xvfb()

%config InlineBackend.print_figure_kwargs={'facecolor' : "w"}
%config InlineBackend.figure_format='retina'
Running this notebook on:  cuda
Last run with spateo version: 1.1.0.dev48+bbacf0a.dirty

Loading the Data

In this tutorial, we will use the Drosophila dataset at stages E14-E16, obtained using the Stereo-Seq technique. You can download the raw data from Flysta3D or from the links below. After downloading, be sure to place the data in the appropriate directory.

[5]:
adata = st.read('./data/E14-16h_a_count_normal_stereoseq.h5ad')

adata
[5]:
AnnData object with n_obs × n_vars = 15295 × 13668
    obs: 'slice_ID', 'raw_x', 'raw_y', 'new_x', 'new_y', 'new_z', 'annotation'
    obsm: 'X_umap', 'spatial'
    layers: 'raw_counts'

Preprocessing

In the first step, we perform preprocessing. Here, we use Dynamo’s preprocessing framework, which is simple and fast.

[6]:
preprocessor = dyn.preprocessing.Preprocessor()
preprocessor.preprocess_adata(adata, recipe="monocle")
|-----> Running monocle preprocessing pipeline...
|-----------> filtered out 145 outlier cells
|-----? No layers exist in adata, skipp filtering by shared counts
|-----------> filtered out 6533 outlier genes
|-----> PCA dimension reduction
|-----> <insert> X_pca to obsm in AnnData Object.
|-----> [Preprocessor-monocle] completed [21.8389s]
[7]:
dyn.tl.reduceDimension(adata, basis="pca")
|-----> retrieve data for non-linear dimension reduction...
|-----? adata already have basis umap. dimension reduction umap will be skipped!
set enforce=True to re-performing dimension reduction.
|-----> Start computing neighbor graph...
|-----------> X_data is None, fetching or recomputing...
|-----> fetching X data from layer:None, basis:pca
|-----> method arg is None, choosing methods automatically...
|-----------> method pynn selected
|-----> [UMAP] completed [29.9284s]
[8]:
# Create a list of colors for each label using rainbow colors
import matplotlib.colors as mcolors
labels = adata.obs['annotation'].cat.categories.tolist()
colors = plt.cm.rainbow(np.linspace(0, 1, len(labels)))
palette = dict(zip(labels, colors))

dyn.pl.umap(adata, color='annotation', pointsize=0.02, alpha=1, color_key=palette)
|-----------> plotting with basis key=X_umap
|-----------> skip filtering annotation by stack threshold when stacking color because it is not a numeric type
../../../_images/tutorials_notebooks_3_alignment_5.1_3D_reconstruction_on_Drosophila_bin_50_data_with_pairwise_alignment_7_1.png

In fact, they have already reconstructed the 3D Drosophila, and the reconstructed spatial coordinates are stored in new_x, new_y, and new_z. To simulate the unaligned case, we generate a random rotation for each slice.

[10]:
adata.obsm['spatial_3D'] = adata.obs[['new_x', 'new_y', 'new_z']].values.copy()
z_height = np.unique(adata.obsm['spatial_3D'][:,2])
slices = [adata[adata.obsm['spatial_3D'][:,2] == z] for z in z_height]
slices[0].obsm['spatial_2D'] = slices[0].obsm['spatial_3D'][:,:2]
for i in range(1, len(slices)):
    slices[i].obsm['spatial_2D'] = slices[i].obsm['spatial_3D'][:,:2].copy()
    # generate random rotation matrix
    st.align.rigid_transformation(slices[i], spatial_key='spatial_2D', key_added='spatial_2D')

The 2D visualization plot show the initial spatial coordinate.

[11]:
spatial_key = 'spatial_2D'
cluster_key = 'annotation'

st.pl.slices_2d(
    slices = slices,
    label_key = cluster_key,
    spatial_key = spatial_key,
    height=2,
    center_coordinate=True,
    show_legend=True,
    legend_kwargs={'loc': 'upper center', 'bbox_to_anchor': (0.5, 0) ,'ncol': 5, 'borderaxespad': -6, 'frameon': False},
    palette=palette,
)
../../../_images/tutorials_notebooks_3_alignment_5.1_3D_reconstruction_on_Drosophila_bin_50_data_with_pairwise_alignment_11_0.png

Pairwise rigid alignment

Next, we calculate the spatial transformation between each pair of adjacent slices through pairwise rigid alignment and save these transformations in transformation (with a length of len(slices)-1) using st.align.morpho_align_transformation.

[12]:
key_added = 'align_spatial'
transformation = st.align.morpho_align_transformation(
    models=slices,
    spatial_key=spatial_key,
    key_added=key_added,
    device=device,
    verbose=False,
    rep_layer='X_pca',
    rep_field='obsm',
    dissimilarity='cos',
)
|-----> [Models alignment based on morpho.] in progress: 100.0000%
|-----> [Models alignment based on morpho.] finished [29.4830s]

The output transformations include both a rotation matrix and a translation matrix

[16]:
transformation[0]
[16]:
{'Rotation': array([[-0.53751908,  0.84325158],
        [-0.84325158, -0.53751908]]),
 'Translation': array([0.71284574, 0.31073377])}

We then align these spatial transformations to unify the coordinate system of all slices, aligning them with the first slice (which remains unchanged).

[13]:
aligned_slices = st.align.morpho_align_apply_transformation(
    models=slices,
    spatial_key=spatial_key,
    key_added=key_added,
    transformation=transformation,
)
|-----> [Models alignment based on morpho, applying transformation.] in progress: 100.0000%
|-----> [Models alignment based on morpho, applying transformation.] finished [0.0098s]

Plot the aligned slices

After the alignment, we can visualize the aligned slices to check whether the pairwise alignment results are correct.

[17]:
st.pl.slices_2d(
    slices = aligned_slices,
    label_key = cluster_key,
    spatial_key = key_added,
    height=2,
    center_coordinate=False,
    show_legend=True,
    legend_kwargs={'loc': 'upper center', 'bbox_to_anchor': (0.5, 0) ,'ncol': 5, 'borderaxespad': -6, 'frameon': False},
    palette=palette,
)
../../../_images/tutorials_notebooks_3_alignment_5.1_3D_reconstruction_on_Drosophila_bin_50_data_with_pairwise_alignment_19_0.png

Since the Drosophila data is already aligned, we can compare the alignment results. Here, we select slices #5, #6, and #7 and overlay them. We can see that Spateo’s aligned results are better.

[21]:
unaligned_slices_overlap_plot = ad.concat(aligned_slices[5:8])
unaligned_slices_overlap_plot.obsm['plot_spatial'] = unaligned_slices_overlap_plot.obsm['spatial_2D']
unaligned_slices_overlap_plot.obs['title'] = 'Unaligned'

aligned_slices_overlap_plot = unaligned_slices_overlap_plot.copy()
aligned_slices_overlap_plot.obsm['plot_spatial'] = aligned_slices_overlap_plot.obsm['align_spatial']
aligned_slices_overlap_plot.obs['title'] = 'Spateo aligned'

previous_aligned_slices_overlap_plot = unaligned_slices_overlap_plot.copy()
previous_aligned_slices_overlap_plot.obsm['plot_spatial'] = previous_aligned_slices_overlap_plot.obsm['spatial_3D'][:,:2]
previous_aligned_slices_overlap_plot.obs['title'] = 'Previous aligned'

st.pl.slices_2d(
    slices = [unaligned_slices_overlap_plot, previous_aligned_slices_overlap_plot, aligned_slices_overlap_plot],
    label_key = cluster_key,
    spatial_key = 'plot_spatial',
    height=3,
    center_coordinate=False,
    show_legend=False,
    palette=palette,
    point_size=2,
    slices_key='title',
)
../../../_images/tutorials_notebooks_3_alignment_5.1_3D_reconstruction_on_Drosophila_bin_50_data_with_pairwise_alignment_21_0.png

3D Visualization

In the end, we concatenate the aligned slices and restore the z-axis data to obtain the 3D reconstruction result using pairwise rigid alignment. The 3D visualization is then performed using Spateo’s built-in plotting functions.

[22]:
aligned_adata = ad.concat(aligned_slices)
aligned_adata.obsm['aligned_spatial_3D'] = np.concatenate([aligned_adata.obsm['align_spatial'], np.array(aligned_adata.obsm['spatial_3D'][:,2])[:,None]], axis=1)
[ ]:
# save the results
aligned_adata.write_h5ad('./results/Drosophila_E14-16h_rigid_alignment.h5ad')
[23]:
aligned_pc, _ = st.tdr.construct_pc(adata=aligned_adata,spatial_key="aligned_spatial_3D",groupby=cluster_key,key_added="tissue",colormap=palette,)

aligned_pc
[23]:
HeaderData Arrays
PolyDataInformation
N Cells15150
N Points15150
N Strips0
X Bounds-1.517e+01, 1.631e+01
Y Bounds-2.867e+01, 3.068e+01
Z Bounds1.000e+00, 1.600e+01
N Arrays3
NameFieldTypeN CompMinMax
tissue_rgbaPointsfloat3240.000e+001.000e+00
tissuePoints1nannan
obs_indexPoints1nannan

Warning:

The following code block involves PyVista visualization. Before running the PyVista-based visualization function st.pl.three_d_plot, please ensure that PyVista is installed correctly, especially for server configurations. Incorrect installation may cause Jupyter Notebook to crash. There are certain specifications to consider when running PyVista on a remote server. We recommend reviewing their tutorials, particularly for remote Jupyter environments.

Tips:

You can enable interactive visualization by changing the PyVista backend via the jupyter parameter in st.pl.three_d_plot. Setting jupyter='static' will output a static image. I recommend using 'html' for lightweight interactive visualization. Another very convenient feature is that you can save the interactive 3D visualization as an HTML file by specifying plotter_filename in st.pl.three_d_plot, making it much easier to share results

st.pl.three_d_plot provides an easy way to visualize 3D plots.

[24]:
# show the aligned drosophila in 3D
st.pl.three_d_plot(model=aligned_pc,key="tissue",model_style='points',model_size=8, show_axes=True,
    jupyter="static",window_size=(1200, 1200),show_outline=True,outline_kwargs={'show_labels': False, 'outline_width': 3},
    # plotter_filename='./results/Drosophila_8_9_rigid_alignment_reconstruction.html',  # uncomment this for saving the interactive 3D visualization
)
../../../_images/tutorials_notebooks_3_alignment_5.1_3D_reconstruction_on_Drosophila_bin_50_data_with_pairwise_alignment_29_0.png

We can also compare with the spatial provided by the data itself.

[25]:
previous_pc, _ = st.tdr.construct_pc(adata=aligned_adata,spatial_key="spatial_3D",groupby=cluster_key,key_added="tissue",colormap=palette,)
[29]:
st.pl.three_d_multi_plot(model=st.tdr.collect_models([previous_pc, aligned_pc]),key="tissue",model_style='points',model_size=4,
    jupyter="static",window_size=(500, 500),show_outline=True,outline_kwargs={'show_labels': False, 'outline_width': 1},
)
../../../_images/tutorials_notebooks_3_alignment_5.1_3D_reconstruction_on_Drosophila_bin_50_data_with_pairwise_alignment_32_0.png

Let’s highlight some of the cluster to show the reconstructed 3D patterns

[31]:
highlight_tissues = ['CNS', 'midgut', 'salivary gland', ]
pc_origin_highlight = previous_pc.copy()
pc_highlight = aligned_pc.copy()
pc_highlight['tissue_rgba'][:,3] = 0.02
pc_highlight['tissue_rgba'][aligned_adata.obs[cluster_key].isin(highlight_tissues),3] = 1
pc_origin_highlight['tissue_rgba'][:,3] = 0.02
pc_origin_highlight['tissue_rgba'][aligned_adata.obs[cluster_key].isin(highlight_tissues),3] = 1

Because the coordinates are the same. For better comparison, we rotate the point cloud (this will not affect AnnData it self)

[32]:
pc_origin_highlight = pc_origin_highlight.rotate_z(205)
[34]:
# show the aligned drosophila in 3D
st.pl.three_d_multi_plot(model=st.tdr.collect_models([pc_origin_highlight, pc_origin_highlight, pc_origin_highlight, pc_highlight, pc_highlight, pc_highlight]),
                         key="tissue",model_style='points',model_size=4, cpo=['xy', 'xz', 'yz', 'xy', 'xz', 'yz', ],   jupyter="static",window_size=(500, 500),
                         text=['Previous xy', 'Previous xz', 'Previous yz', 'Spateo xy', 'Spateo xz', 'Spateo yz']
)
../../../_images/tutorials_notebooks_3_alignment_5.1_3D_reconstruction_on_Drosophila_bin_50_data_with_pairwise_alignment_37_0.png