Alignment_Paste_DLPFC

This tutorial demonstrates alignment on DLPFC data using SODB and Paste.

A reference paper can be found at https://www.nature.com/articles/s41592-022-01459-6.

This tutorial refers to the following tutorial at https://github.com/raphael-group/paste_reproducibility/blob/main/notebooks/DLPFC_pairwise.ipynb. At the same time, the way of loadding data is modified by using SODB.

# imports various libraries and packages for data analysis and visualization.
# math: provides mathematical functions such as logarithms, trigonometric functions, etc.
# pandas: used for data manipulation and analysis.
# numpy: used for numerical computing, including mathematical operations on arrays and matrices.
# scipy: used for scientific computing, including functions for optimization, linear algebra, statistics, and signal processing.
# seaborn: used for statistical data visualization, providing high-level interfaces for creating informative and attractive visualizations.
# matplotlib: a comprehensive library for creating static, animated, and interactive visualizations in Python.
# matplotlib.patches: provides classes for creating graphical objects such as rectangles, circles, and polygons.
# style: a module within matplotlib that allows users to customize the style of plots.
# time: provides time-related functions, such as measuring execution time and converting between time formats.
# scanpy: a Python package for single-cell gene expression analysis, including preprocessing, clustering, and differential expression analysis.
# sklearn: a machine learning library with various tools for classification, regression, clustering, and dimensionality reduction.
# networkx: a Python package for creating, manipulating, and studying complex networks.
# ot: a Python package for optimal transport (OT) computations, including OT-based algorithms for data analysis and visualization.
import math
import pandas as pd
import numpy as np
import scipy
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib import style
import matplotlib
import time
import scanpy as sc
import sklearn
import networkx as nx
import ot

/tmp/ipykernel_92630/1186083132.py:19: MatplotlibDeprecationWarning: The seaborn styles shipped by Matplotlib are deprecated since 3.6, as they no longer correspond to the styles shipped by seaborn. However, they will remain available as 'seaborn-v0_8-<style>'. Alternatively, directly use the seaborn API instead.
  style.use('seaborn-white')

# import the paste package
import paste as pst

"""
#%load_ext autoreload
#%autoreload 2
#style.use('seaborn-dark')
style.use('seaborn-white')

#报错ModuleNotFoundError: No module named 'ot'
#需要pip install POT
"""

# create a list containing 12 samples
sample_list = ["151507", "151508", "151509","151510", "151669", "151670","151671", "151672", "151673","151674", "151675", "151676"]

# Import the pysodb library
# pysodb is a Python package that provides a set of tools for working with SODB (Store On-Demand Block) databases.
# SODB is a format used to store data in memory-mapped files for efficient access and querying.
# This library allows users to interact with SODB files using Python.
import pysodb

# Initialization
sodb = pysodb.SODB()

# Get the list of datasets with specific category.
# categories ["Spatial Transcriptomics", "Spatial Proteomics", "Spatial Metabolomics", "Spatial Genomics", "Spatial MultiOmics"]
sodb.list_dataset_by_category('Spatial Transcriptomics')

['moffitt2018molecular',
 'chen2022spatiotemporal',
 'Vickovic2019high',
 'eng2019transcriptome',
 'carlberg2019exploring',
 'chen2022spatiotemporal_compre_20',
 'Navarro2020Spatial',
 'gracia2021genome',
 'zhang2021spatially',
 'xia2019spatial',
 'mantri2021spatiotemporal',
 'Merfish_Visp',
 'he2020integrating',
 'Shah2016InSitu',
 'maynard2021trans',
 'DARTFISH',
 'fawkner2021spatiotemporal',
 'ortiz2020molecular',
 'Ratz2022Clonal',
 '10x',
 'Juntaro2022MEK',
 'wei2022single',
 'Lebrigand2022The',
 'Tower2021Spatial',
 'liu2022spatiotemporal',
 'scispace',
 'bergenstrahle2021super',
 'Marshall2022High_human',
 'Gouin2021An',
 'wang2022high',
 'backdahl2021spatial',
 'Wang2018Three_1k',
 'Dixon2022Spatially',
 'hildebrandt2021spatial',
 'parigi2022the',
 'Alon2021Expansion',
 'Pascual2021Dietary',
 'Sun2021Integrating',
 'Buzzi2022Spatial',
 'asp2017spatial',
 'Sanchez2021A',
 'Fang2022Conservation',
 'rodriques2019slide',
 'seqFISH_VISp',
 'moncada2020integrating',
 'Biermann2022Dissecting',
 'Kleshchevnikov2022Cell2location',
 'Garcia2021Mapping',
 'chen2021decoding',
 'kvastad2021the',
 'stickels2020highly',
 'guilliams2022spatial',
 'Joglekar2021A',
 'wang2021easi',
 'xia2022the',
 'Misra2021Characterizing',
 'Marshall2022High_mouse',
 'maniatis2019spatiotemporal',
 'Barkley2022Cancer',
 'thrane2018spatially',
 'Sun2022Excitatory',
 'lohoff2021integration',
 'Dhainaut2022Spatial',
 'Borm2022Scalable',
 'Wang2018three',
 'Vickovic2019high_update',
 'chen2021dissecting',
 'Konieczny2022Interleukin',
 'asp2019a',
 'Kadur2022Human',
 'ji2020multimodal',
 'stahl2016visualization',
 'berglund2018spatial',
 'Visium_Allen',
 'hunter2021spatially',
 'chen2020spatial',
 'Wu2022spatial',
 'Melo2021Integrating',
 'codeluppi2018spatial']

# Get the list of datasets
adata_list = sodb.load_dataset('maynard2021trans')

load experiment[151508] in dataset[maynard2021trans]
load experiment[151671] in dataset[maynard2021trans]
load experiment[151507] in dataset[maynard2021trans]
load experiment[151674] in dataset[maynard2021trans]
load experiment[151670] in dataset[maynard2021trans]
load experiment[151669] in dataset[maynard2021trans]
load experiment[151676] in dataset[maynard2021trans]
load experiment[151675] in dataset[maynard2021trans]
load experiment[151509] in dataset[maynard2021trans]
load experiment[151673] in dataset[maynard2021trans]
load experiment[151672] in dataset[maynard2021trans]
load experiment[151510] in dataset[maynard2021trans]

# creates a new dictionary called "adatas" by removing observations with missing values in the "Region" column from each dataset in the original dictionary "adata_list".
adatas = {}
for key in adata_list.keys():
    a = adata_list[key]
    a = a[np.logical_not(a.obs['Region'].isna())]
    adatas[key] = a

# defines groups of samples based on IDs
sample_groups = [["151507", "151508", "151509","151510"],[ "151669", "151670","151671", "151672"],[ "151673","151674", "151675", "151676"]]
# creates a list of lists called layer_groups
# where each sub-list contains the AnnData objects for each sample in the corresponding group from sample_groups.
layer_groups = [[adatas[sample_groups[j][i]] for i in range(len(sample_groups[j]))] for j in range(len(sample_groups))]
# create a dictionary that maps the layer number to a color from the default Seaborn color palette.
# The layer number is represented as a string in the format "Layer{layer_number}"
layer_to_color_map = {'Layer{0}'.format(i+1):sns.color_palette()[i] for i in range(6)}
# adds an additional key "WM" to the layer_to_color_map dictionary
layer_to_color_map['WM'] = sns.color_palette()[6]

# Visualize the different slices of the DLPFC brain region mapped by layer_groups
slice_map = {0:'A',1:'B',2:'C',3:'D'}
fig, axs = plt.subplots(3, 4,figsize=(15,11.5))
for j in range(len(layer_groups)):
    axs[j,0].text(-0.1, 0.5, 'Sample '+slice_map[j],fontsize=12,rotation='vertical',transform = axs[j,0].transAxes,verticalalignment='center')
    for i in range(len(layer_groups[j])):
        adata = adatas[sample_list[j*4+i]]
        colors = list(adata.obs['Region'].astype('str').map(layer_to_color_map))
        colors = [(r, g, b) for r, g, b in colors]
        axs[j,i].scatter(layer_groups[j][i].obsm['spatial'][:,0],layer_groups[j][i].obsm['spatial'][:,1],linewidth=0,s=20, marker=".",
                         color=colors
                        )
        axs[j,i].set_title('Slice '+ slice_map[i],size=12)
        axs[j,i].invert_yaxis()
        axs[j,i].axis('off')
        if i<3:
            s = '300$\mu$m' if i==1 else "10$\mu$m"
            delta = 0.05 if i==1 else 0
            axs[j,i].annotate('',xy=(1-delta, 0.5), xytext=(1.2+delta, 0.5),xycoords=axs[j,i].transAxes,textcoords=axs[j,i].transAxes,arrowprops=dict(arrowstyle='<->',lw=1))
            axs[j,0].text(1.1, 0.55, s,fontsize=9,transform = axs[j,i].transAxes,horizontalalignment='center')
    axs[j,3].legend(handles=[mpatches.Patch(color=layer_to_color_map[adata.obs['Region'].cat.categories[i]], label=adata.obs['Region'].cat.categories[i]) for i in range(len(adata.obs['Region'].cat.categories))],fontsize=10,title='Cortex layer',title_fontsize=12,bbox_to_anchor=(1, 1))

plt.savefig('figures/DLPFC_before.pdf',bbox_inches='tight',transparent=True)

# calculates the maximum accuracy of a binary classifier that distinguishes between two sets of labels
def max_accuracy(labels1,labels2):
    w = min(1/len(labels1),1/len(labels2))
    cats = set(pd.unique(labels1)).union(set(pd.unique(labels1)))
    return sum([w * min(sum(labels1==c),sum(labels2==c)) for c in cats])
# calculates the mapping accuracy between two sets of labels using the optimal transport (OT) algorithm
def mapping_accuracy(labels1,labels2,pi):
    mapping_dict = {'Layer1':1, 'Layer2':2, 'Layer3':3, 'Layer4':4, 'Layer5':5, 'Layer6':6, 'WM':7}
    return np.sum(pi*(scipy.spatial.distance_matrix(np.matrix(labels1.map(mapping_dict) ).T,np.matrix(labels2.map(mapping_dict)).T)==0))
#
import itertools
# calculates the optimal transport plan pi between two sets of labels using the OT algorithm
def max_accuracy_mapping(labels1,labels2):
    n1,n2=len(labels1),len(labels2)
    mapping_dict = {'Layer1':1, 'Layer2':2, 'Layer3':3, 'Layer4':4, 'Layer5':5, 'Layer6':6, 'WM':7}
    dist = np.array(scipy.spatial.distance_matrix(np.matrix(labels1.map(mapping_dict)).T,np.matrix(labels2.map(mapping_dict)).T)!=0,dtype=float)
    pi = ot.emd(np.ones(n1)/n1, np.ones(n2)/n2, dist)
    return pi

alpha = 0.1
res_df = pd.DataFrame(columns=['Sample','Pair','Kind','Time','Accuracy'])
pis = [[None for i in range(len(layer_groups[j])-1)] for j in range(len(layer_groups))]
for j in range(len(layer_groups)):
    for i in range(len(layer_groups[j])-1):
        pi0 = pst.match_spots_using_spatial_heuristic(layer_groups[j][i].obsm['spatial'],layer_groups[j][i+1].obsm['spatial'],use_ot=True)
        start = time.time()
        pis[j][i] = pst.pairwise_align(layer_groups[j][i], layer_groups[j][i+1],alpha=alpha,G_init=pi0,norm=True,verbose=False)
        tt = time.time()-start
        acc = mapping_accuracy(layer_groups[j][i].obs['Region'],layer_groups[j][i+1].obs['Region'],pis[j][i])
        print(j,i,'Accuracy',acc,'time',tt)
        # np.savetxt('../data/DLPFC/saved_results/init_{0}_{1}_{2}.gz'.format(j,i,'ot'), pis[j][i], delimiter=',')
        res_df.loc[len(res_df)] = [j,i,'PASTE',tt,acc]

Using selected backend cpu. If you want to use gpu, set use_gpu = True.
0 0 Accuracy 0.815543482357772 time 333.4515166282654
Using selected backend cpu. If you want to use gpu, set use_gpu = True.

RESULT MIGHT BE INACURATE
Max number of iteration reached, currently 100000. Sometimes iterations go on in cycle even though the solution has been reached, to check if it's the case here have a look at the minimal reduced cost. If it is very close to machine precision, you might actually have the correct solution, if not try setting the maximum number of iterations a bit higher
/home/yzy/anaconda3/envs/paste/lib/python3.8/site-packages/ot/lp/__init__.py:343: UserWarning: numItermax reached before optimality. Try to increase numItermax.
  result_code_string = check_result(result_code)
RESULT MIGHT BE INACURATE
Max number of iteration reached, currently 100000. Sometimes iterations go on in cycle even though the solution has been reached, to check if it's the case here have a look at the minimal reduced cost. If it is very close to machine precision, you might actually have the correct solution, if not try setting the maximum number of iterations a bit higher
RESULT MIGHT BE INACURATE
Max number of iteration reached, currently 100000. Sometimes iterations go on in cycle even though the solution has been reached, to check if it's the case here have a look at the minimal reduced cost. If it is very close to machine precision, you might actually have the correct solution, if not try setting the maximum number of iterations a bit higher
RESULT MIGHT BE INACURATE
Max number of iteration reached, currently 100000. Sometimes iterations go on in cycle even though the solution has been reached, to check if it's the case here have a look at the minimal reduced cost. If it is very close to machine precision, you might actually have the correct solution, if not try setting the maximum number of iterations a bit higher

0 1 Accuracy 0.22047228891676612 time 52.269108057022095
Using selected backend cpu. If you want to use gpu, set use_gpu = True.

RESULT MIGHT BE INACURATE
Max number of iteration reached, currently 100000. Sometimes iterations go on in cycle even though the solution has been reached, to check if it's the case here have a look at the minimal reduced cost. If it is very close to machine precision, you might actually have the correct solution, if not try setting the maximum number of iterations a bit higher
/home/yzy/anaconda3/envs/paste/lib/python3.8/site-packages/ot/lp/__init__.py:343: UserWarning: numItermax reached before optimality. Try to increase numItermax.
  result_code_string = check_result(result_code)
RESULT MIGHT BE INACURATE
Max number of iteration reached, currently 100000. Sometimes iterations go on in cycle even though the solution has been reached, to check if it's the case here have a look at the minimal reduced cost. If it is very close to machine precision, you might actually have the correct solution, if not try setting the maximum number of iterations a bit higher

0 2 Accuracy 0.8713495745166362 time 32.520468950271606
Using selected backend cpu. If you want to use gpu, set use_gpu = True.
1 0 Accuracy 0.9151285966713762 time 598.9820036888123
Using selected backend cpu. If you want to use gpu, set use_gpu = True.

RESULT MIGHT BE INACURATE
Max number of iteration reached, currently 100000. Sometimes iterations go on in cycle even though the solution has been reached, to check if it's the case here have a look at the minimal reduced cost. If it is very close to machine precision, you might actually have the correct solution, if not try setting the maximum number of iterations a bit higher
/home/yzy/anaconda3/envs/paste/lib/python3.8/site-packages/ot/lp/__init__.py:343: UserWarning: numItermax reached before optimality. Try to increase numItermax.
  result_code_string = check_result(result_code)
RESULT MIGHT BE INACURATE
Max number of iteration reached, currently 100000. Sometimes iterations go on in cycle even though the solution has been reached, to check if it's the case here have a look at the minimal reduced cost. If it is very close to machine precision, you might actually have the correct solution, if not try setting the maximum number of iterations a bit higher
RESULT MIGHT BE INACURATE
Max number of iteration reached, currently 100000. Sometimes iterations go on in cycle even though the solution has been reached, to check if it's the case here have a look at the minimal reduced cost. If it is very close to machine precision, you might actually have the correct solution, if not try setting the maximum number of iterations a bit higher

1 1 Accuracy 0.5920829519638594 time 30.46793842315674
Using selected backend cpu. If you want to use gpu, set use_gpu = True.
1 2 Accuracy 0.861938266075075 time 152.26987481117249
Using selected backend cpu. If you want to use gpu, set use_gpu = True.
2 0 Accuracy 0.8589111598101083 time 178.80665159225464
Using selected backend cpu. If you want to use gpu, set use_gpu = True.
2 1 Accuracy 0.827013263737229 time 208.154137134552
Using selected backend cpu. If you want to use gpu, set use_gpu = True.
2 2 Accuracy 0.8283188989963737 time 259.1206352710724

res_df

	Sample	Pair	Kind	Time	Accuracy
0	0	0	PASTE	333.451517	0.815543
1	0	1	PASTE	52.269108	0.220472
2	0	2	PASTE	32.520469	0.871350
3	1	0	PASTE	598.982004	0.915129
4	1	1	PASTE	30.467938	0.592083
5	1	2	PASTE	152.269875	0.861938
6	2	0	PASTE	178.806652	0.858911
7	2	1	PASTE	208.154137	0.827013
8	2	2	PASTE	259.120635	0.828319

paste_layer_groups = [pst.stack_slices_pairwise(layer_groups[j], pis[j]) for j in range(len(layer_groups)) ]

def plot_slices_overlap(groups, adatas, sample_list, layer_to_color_map,save=None,):
    marker_list = ['.','*','x','+']

    for j in range(len(groups)):
        plt.figure(figsize=(10,10))
        for i in range(len(groups[j])):
            adata = adatas[sample_list[j*4+i]]
            colors = list(adata.obs['Region'].astype('str').map(layer_to_color_map))
            plt.scatter(groups[j][i].obsm['spatial'][:,0],groups[j][i].obsm['spatial'][:,1],linewidth=1,s=80, marker=marker_list[i],color=colors,alpha=0.7)
        plt.legend(handles=[mpatches.Patch(color=layer_to_color_map[adata.obs['Region'].cat.categories[i]], label=adata.obs['Region'].cat.categories[i]) for i in range(len(adata.obs['Region'].cat.categories))],fontsize=10,title='Cortex layer',title_fontsize=15,bbox_to_anchor=(1, 1))
        plt.gca().invert_yaxis()
        plt.axis('off')
        if save is None:
            plt.show()
        else:
            plt.savefig(f'{save}_{j}.pdf',bbox_inches='tight',transparent=True)

# Plot Stacking of Four slices without alignment
plot_slices_overlap(layer_groups, adatas, sample_list, layer_to_color_map,save='figures/DLPFC_before')

# Plot Stacking of Four slices with PASTE alignment
plot_slices_overlap(paste_layer_groups, adatas, sample_list, layer_to_color_map,save='figures/DLPFC_after')

dataset_name = 'dlpfc'
for i in range(len(layer_groups)):
    rsta_list = []

    for j in range(len(layer_groups[i])):
        a= layer_groups[i][j]
        spatial = a.obsm['spatial']
        spatial_z = np.ones(shape=(spatial.shape[0],1))*j*100
        spatial_3d = np.hstack([spatial,spatial_z])
        a.obsm['spatial_3d'] = spatial_3d
        rsta_list.append(a)
    a_concat = rsta_list[0].concatenate(rsta_list[1:])
    a_concat.write_h5ad(f'{dataset_name}_sample{i}_3d.h5ad')

/home/yzy/anaconda3/envs/paste/lib/python3.8/site-packages/anndata/_core/anndata.py:1785: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.
  [AnnData(sparse.csr_matrix(a.shape), obs=a.obs) for a in all_adatas],
/home/yzy/anaconda3/envs/paste/lib/python3.8/site-packages/anndata/_core/anndata.py:1785: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.
  [AnnData(sparse.csr_matrix(a.shape), obs=a.obs) for a in all_adatas],
/home/yzy/anaconda3/envs/paste/lib/python3.8/site-packages/anndata/_core/anndata.py:1785: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.
  [AnnData(sparse.csr_matrix(a.shape), obs=a.obs) for a in all_adatas],
/home/yzy/anaconda3/envs/paste/lib/python3.8/site-packages/anndata/_core/anndata.py:1785: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.
  [AnnData(sparse.csr_matrix(a.shape), obs=a.obs) for a in all_adatas],
/home/yzy/anaconda3/envs/paste/lib/python3.8/site-packages/anndata/_core/anndata.py:1785: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.
  [AnnData(sparse.csr_matrix(a.shape), obs=a.obs) for a in all_adatas],
/home/yzy/anaconda3/envs/paste/lib/python3.8/site-packages/anndata/_core/anndata.py:1785: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.
  [AnnData(sparse.csr_matrix(a.shape), obs=a.obs) for a in all_adatas],
/home/yzy/anaconda3/envs/paste/lib/python3.8/site-packages/anndata/_core/anndata.py:1785: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.
  [AnnData(sparse.csr_matrix(a.shape), obs=a.obs) for a in all_adatas],
/home/yzy/anaconda3/envs/paste/lib/python3.8/site-packages/anndata/_core/anndata.py:1785: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.
  [AnnData(sparse.csr_matrix(a.shape), obs=a.obs) for a in all_adatas],
/home/yzy/anaconda3/envs/paste/lib/python3.8/site-packages/anndata/_core/anndata.py:1785: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.
  [AnnData(sparse.csr_matrix(a.shape), obs=a.obs) for a in all_adatas],
/home/yzy/anaconda3/envs/paste/lib/python3.8/site-packages/anndata/_core/anndata.py:1785: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.
  [AnnData(sparse.csr_matrix(a.shape), obs=a.obs) for a in all_adatas],
/home/yzy/anaconda3/envs/paste/lib/python3.8/site-packages/anndata/_core/anndata.py:1785: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.
  [AnnData(sparse.csr_matrix(a.shape), obs=a.obs) for a in all_adatas],
/home/yzy/anaconda3/envs/paste/lib/python3.8/site-packages/anndata/_core/anndata.py:1785: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.
  [AnnData(sparse.csr_matrix(a.shape), obs=a.obs) for a in all_adatas],

dataset_name = 'dlpfc'
for i in range(len(paste_layer_groups)):
    rsta_list = []

    for j in range(len(paste_layer_groups[i])):
        a= paste_layer_groups[i][j]
        spatial = a.obsm['spatial']
        spatial_z = np.ones(shape=(spatial.shape[0],1))*j*100
        spatial_3d = np.hstack([spatial,spatial_z])
        a.obsm['spatial_3d'] = spatial_3d
        rsta_list.append(a)
    a_concat = rsta_list[0].concatenate(rsta_list[1:])
    a_concat.write_h5ad(f'{dataset_name}_sample{i}_3d.h5ad')

8517.0