Mouse Brain (Visium, 3D alignment of multiple slices)
[2]:
import numpy as np
import pandas as pd
import pathlib
import matplotlib.pyplot as plt
import matplotlib as mpl
import seaborn as sns
import scanpy as sc
import sys
sys.path.append('../src')
from utils import *
import warnings
warnings.filterwarnings('ignore')
Preprocessing scRNA-ref
Ideally, the scRNAseq reference should from the exact same tissue that contains all the celltypes you expect to find in ST spots.
We will use a reference scRNA-seq dataset of ~14,000 adult mouse cortical cell taxonomy from the Allen Institute, generated with the SMART-Seq2 protocol.
The raw dataset in h5ad format (adata_processed_sc.h5ad) provided by scanpy is available in here.
Conveniently, we provided the pre-pocessed scRNA-ref (Ckpts_scRefs/VISp/Ref_scRNA_VISp_qc2_2Kgenes.h5ad) as well as other relevent materials involved in the following example in here, so you can skip this part
[20]:
ad_sc = sc.read('adata_processed_sc.h5ad')
ad_sc = ad_sc[ad_sc.obs['dissected_region']=='VISp']
ad_sc = ad_sc[ad_sc.obs['cell_subclass']!='nan']
ad_sc = ad_sc[ad_sc.obs['cell_subclass']!='No Class']
ad_sc.obs['cell_subclass'] = ad_sc.obs['cell_subclass'].astype(str)
ad_sc.obs.loc[ad_sc.obs.cell_subclass=='NP','cell_subclass'] = 'L5 NP'
cell_type_column = 'cell_subclass'
ad_sc.X.max(),ad_sc.shape
[20]:
(13.006981, (13360, 36577))
[21]:
# transform to count scale
if ad_sc.X.max()<20:
ad_sc.X = np.exp(ad_sc.X)-1
sc.pp.normalize_total(ad_sc, target_sum=1e5)
sc.pp.filter_genes(ad_sc, min_cells=10)
ad_sc = ad_sc[:,~np.array([_.startswith('MT-') for _ in ad_sc.var.index])]
ad_sc = ad_sc[:,~np.array([_.startswith('mt-') for _ in ad_sc.var.index])]
ad_sc.shape
[21]:
(13360, 34235)
[22]:
ad_sc.X.max(),ad_sc.X.min()
[22]:
(ArrayView(25776.516, dtype=float32), ArrayView(0., dtype=float32))
[23]:
sc.pp.highly_variable_genes(ad_sc, flavor='seurat_v3',n_top_genes=1000)
ad_sc.uns['log1p'] = {}
ad_sc.uns['log1p']['base'] = None
sc.tl.rank_genes_groups(ad_sc, groupby=cell_type_column, method='wilcoxon')
markers_df = pd.DataFrame(ad_sc.uns["rank_genes_groups"]["names"]).iloc[0:30, :]
markers = list(np.unique(markers_df.melt().value.values))
markers = list(set(ad_sc.var.loc[ad_sc.var['highly_variable']==1].index)|set(markers)) # highly variable genes + cell type marker genes
len(markers)
[23]:
1389
[24]:
ligand_recept = list(set(pd.read_csv('../extdata/ligand_receptors.txt',sep='\t').melt()['value'].values))
# if scRNA-seq reference is from human tissue, run following code to make gene name consistent
# ligand_recept = [_.upper() for _ in ligand_recept]
markers = markers+ligand_recept
ad_sc.var.loc[ad_sc.var.index.isin(markers),'Marker'] = True
ad_sc.var['Marker'] = ad_sc.var['Marker'].fillna(False)
ad_sc.var['highly_variable'] = ad_sc.var['Marker']
[8]:
sc.pp.log1p(ad_sc)
sc.pp.pca(ad_sc,svd_solver='arpack', n_comps=30, use_highly_variable=True)
ad_sc.X.max()
WARNING: adata.X seems to be already log-transformed.
[8]:
10.157258
[9]:
sc.pp.neighbors(ad_sc, metric='cosine', n_neighbors=30, n_pcs = 30)
sc.tl.umap(ad_sc, min_dist = 0.3, spread = 1, maxiter=100)
[10]:
fig, axs = plt.subplots(1, 1, figsize=(10, 10))
sc.pl.umap(
ad_sc, color=cell_type_column, size=10, frameon=False, show=False, ax=axs,legend_loc='on data'
)
plt.tight_layout()
[12]:
#ad_sc.write('../Ckpts_scRefs/VISp/Ref_scRNA_VISp_qc2_2Kgenes.h5ad')
[ ]:
Multiple ST slices 3D alignment using PASTE
[16]:
import paste as pst
from scipy.spatial.distance import pdist
import matplotlib.patches as mpatches
Load data
[351]:
slice1_file = '../demo_data/Visium_MouseBrain_Cortex_section1.h5ad'
slice2_file = '../demo_data/Visium_MouseBrain_Cortex_section2.h5ad'
# These two .h5ad file had been performed Step 1 (Nuclei segmentation), see more detial in Human-Heart.ipynb
slice1 = sc.read(slice1_file)
slice2 = sc.read(slice2_file)
/home/jxiaoae/.conda/envs/SpatialScope/lib/python3.9/site-packages/anndata/_core/anndata.py:1830: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.
utils.warn_names_duplicates("var")
/home/jxiaoae/.conda/envs/SpatialScope/lib/python3.9/site-packages/anndata/_core/anndata.py:1830: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.
utils.warn_names_duplicates("var")
[263]:
slice1.uns['cell_locations'].head()
[263]:
| x | y | spot_index | cell_index | cell_nums | |
|---|---|---|---|---|---|
| 0 | 7910.500000 | 3181.500000 | AAACAGAGCGACTCCT-1 | AAACAGAGCGACTCCT-1_0 | 8 |
| 1 | 7915.866197 | 3149.855634 | AAACAGAGCGACTCCT-1 | AAACAGAGCGACTCCT-1_1 | 8 |
| 2 | 7942.919540 | 3123.827586 | AAACAGAGCGACTCCT-1 | AAACAGAGCGACTCCT-1_2 | 8 |
| 3 | 7948.000000 | 3205.769231 | AAACAGAGCGACTCCT-1 | AAACAGAGCGACTCCT-1_3 | 8 |
| 4 | 7939.994505 | 3180.730769 | AAACAGAGCGACTCCT-1 | AAACAGAGCGACTCCT-1_4 | 8 |
[264]:
fig, ax = plt.subplots(1,1,figsize=(12, 6),dpi=100)
PlotVisiumCells(slice1,"cell_nums",0.23,0.35,0.3,ax=ax)
/import/home/jxiaoae/spatial/SpatialScope-Beta/demos/../src/utils.py:77: 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.
test = sc.AnnData(np.zeros(merged_df.shape), obs=merged_df)
/home/jxiaoae/.conda/envs/SpatialScope/lib/python3.9/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.
warnings.warn("Transforming to str index.", ImplicitModificationWarning)
[265]:
fig, ax = plt.subplots(1,1,figsize=(12, 6),dpi=100)
PlotVisiumCells(slice2,"cell_nums",0.23,0.35,0.3,ax=ax)
/import/home/jxiaoae/spatial/SpatialScope-Beta/demos/../src/utils.py:77: 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.
test = sc.AnnData(np.zeros(merged_df.shape), obs=merged_df)
/home/jxiaoae/.conda/envs/SpatialScope/lib/python3.9/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.
warnings.warn("Transforming to str index.", ImplicitModificationWarning)
[ ]:
Align two slices with PASTE
[266]:
slice1.obs.index = ['slice1-'+_ for _ in slice1.obs.index]
slice2.obs.index = ['slice2-'+_ for _ in slice2.obs.index]
slice1.uns['cell_locations']['slice_index'] = 'slice1'
slice2.uns['cell_locations']['slice_index'] = 'slice2'
slice1.uns['cell_locations']['spot_index'] = slice1.uns['cell_locations']['spot_index'].map(lambda x:'slice1-'+x)
slice2.uns['cell_locations']['spot_index'] = slice2.uns['cell_locations']['spot_index'].map(lambda x:'slice2-'+x)
slice1.uns['cell_locations']['cell_index'] = slice1.uns['cell_locations']['cell_index'].map(lambda x:'slice1-'+x)
slice2.uns['cell_locations']['cell_index'] = slice2.uns['cell_locations']['cell_index'].map(lambda x:'slice2-'+x)
slice1.obs_names_make_unique()
slice2.obs_names_make_unique()
slice1.var_names_make_unique()
slice2.var_names_make_unique()
[267]:
slice1.uns['cell_locations'].head()
[267]:
| x | y | spot_index | cell_index | cell_nums | slice_index | |
|---|---|---|---|---|---|---|
| 0 | 7910.500000 | 3181.500000 | slice1-AAACAGAGCGACTCCT-1 | slice1-AAACAGAGCGACTCCT-1_0 | 8 | slice1 |
| 1 | 7915.866197 | 3149.855634 | slice1-AAACAGAGCGACTCCT-1 | slice1-AAACAGAGCGACTCCT-1_1 | 8 | slice1 |
| 2 | 7942.919540 | 3123.827586 | slice1-AAACAGAGCGACTCCT-1 | slice1-AAACAGAGCGACTCCT-1_2 | 8 | slice1 |
| 3 | 7948.000000 | 3205.769231 | slice1-AAACAGAGCGACTCCT-1 | slice1-AAACAGAGCGACTCCT-1_3 | 8 | slice1 |
| 4 | 7939.994505 | 3180.730769 | slice1-AAACAGAGCGACTCCT-1 | slice1-AAACAGAGCGACTCCT-1_4 | 8 | slice1 |
[268]:
slice2.uns['cell_locations'].head()
[268]:
| x | y | spot_index | cell_index | cell_nums | slice_index | |
|---|---|---|---|---|---|---|
| 0 | 8007.241667 | 3293.700000 | slice2-AAACAGAGCGACTCCT-1 | slice2-AAACAGAGCGACTCCT-1_0 | 9 | slice2 |
| 1 | 8023.668033 | 3280.393443 | slice2-AAACAGAGCGACTCCT-1 | slice2-AAACAGAGCGACTCCT-1_1 | 9 | slice2 |
| 2 | 8070.021978 | 3271.000000 | slice2-AAACAGAGCGACTCCT-1 | slice2-AAACAGAGCGACTCCT-1_2 | 9 | slice2 |
| 3 | 8051.697674 | 3271.736434 | slice2-AAACAGAGCGACTCCT-1 | slice2-AAACAGAGCGACTCCT-1_3 | 9 | slice2 |
| 4 | 8058.884444 | 3330.088889 | slice2-AAACAGAGCGACTCCT-1 | slice2-AAACAGAGCGACTCCT-1_4 | 9 | slice2 |
[269]:
# Pairwise align the slices
pi12 = pst.pairwise_align(slice1, slice2)
Using selected backend cpu. If you want to use gpu, set use_gpu = True.
[270]:
slices, pis = [slice1, slice2], [pi12]
new_slices = pst.stack_slices_pairwise(slices, pis)
[271]:
slice_colors = ['#e41a1c','#377eb8']
slice_colors = ['teal', 'orange']
slice_colors = ['#377eb8', 'orange']
plt.figure(figsize=(11,7))
for i in range(len(new_slices)):
pst.plot_slice(new_slices[i],slice_colors[i],s=350)
plt.legend(handles=[mpatches.Patch(color=slice_colors[0], label='Cortex slice 1'),mpatches.Patch(color=slice_colors[1], label='Cortex slice 2')],
fontsize=29,bbox_to_anchor=(.96,.72))
plt.gca().invert_yaxis()
plt.axis('off')
plt.show()
[272]:
new_slices[0].obsm['spatial'] = np.hstack((new_slices[0].obsm['spatial'],
np.zeros((new_slices[0].obsm['spatial'].shape[0],1))))
new_slices[1].obsm['spatial'] = np.hstack((new_slices[1].obsm['spatial'],
(pdist(new_slices[0].obsm['spatial']).min()**2/2)**.5/2 * np.ones((new_slices[1].obsm['spatial'].shape[0],1))))
[273]:
slice1.obsm['spatial_align'] = new_slices[0].obsm['spatial']
slice2.obsm['spatial_align'] = new_slices[1].obsm['spatial']
[274]:
slice1.obsm['spatial']
[274]:
array([[7950, 3163],
[9257, 4241],
[5747, 4361],
...,
[5610, 3642],
[5885, 4361],
[6780, 5199]])
[275]:
slice1.obsm['spatial_align']
[275]:
array([[ 1592.37822014, -1400.27868852, 0. ],
[ 2899.37822014, -322.27868852, 0. ],
[ -610.62177986, -202.27868852, 0. ],
...,
[ -747.62177986, -921.27868852, 0. ],
[ -472.62177986, -202.27868852, 0. ],
[ 422.37822014, 635.72131148, 0. ]])
[276]:
def align_cell_locations(slice1):
spot_indexs = slice1.obs.index.values
cell_locations_align = slice1.uns['cell_locations'].copy()
for row in cell_locations_align.iterrows():
spot_id = np.where(spot_indexs==row[1]['spot_index'])[0][0]
cell_locations_align.loc[row[0],'x'] = row[1]['x'] - slice1.obsm['spatial'][spot_id,0] + slice1.obsm['spatial_align'][spot_id,0]
cell_locations_align.loc[row[0],'y'] = row[1]['y'] - slice1.obsm['spatial'][spot_id,1] + slice1.obsm['spatial_align'][spot_id,1]
cell_locations_align['z'] = slice1.obsm['spatial_align'][spot_id,2]
slice1.uns['cell_locations'] = cell_locations_align
align_cell_locations(slice1)
align_cell_locations(slice2)
[277]:
slice1.uns['cell_locations'].head()
[277]:
| x | y | spot_index | cell_index | cell_nums | slice_index | z | |
|---|---|---|---|---|---|---|---|
| 0 | 1552.878220 | -1381.778689 | slice1-AAACAGAGCGACTCCT-1 | slice1-AAACAGAGCGACTCCT-1_0 | 8 | slice1 | 0.0 |
| 1 | 1558.244417 | -1413.423055 | slice1-AAACAGAGCGACTCCT-1 | slice1-AAACAGAGCGACTCCT-1_1 | 8 | slice1 | 0.0 |
| 2 | 1585.297760 | -1439.451102 | slice1-AAACAGAGCGACTCCT-1 | slice1-AAACAGAGCGACTCCT-1_2 | 8 | slice1 | 0.0 |
| 3 | 1590.378220 | -1357.509458 | slice1-AAACAGAGCGACTCCT-1 | slice1-AAACAGAGCGACTCCT-1_3 | 8 | slice1 | 0.0 |
| 4 | 1582.372726 | -1382.547919 | slice1-AAACAGAGCGACTCCT-1 | slice1-AAACAGAGCGACTCCT-1_4 | 8 | slice1 | 0.0 |
[278]:
spatial_align = np.vstack((slice1.obsm['spatial_align'],slice2.obsm['spatial_align']))
cell_locations_align = slice1.uns['cell_locations'].append(slice2.uns['cell_locations']).reset_index(drop=True)
/home/share/xiaojs/software/netMHCpan-4.1/Linux_x86_64/tmp/ipykernel_4891/1740432794.py:2: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
cell_locations_align = slice1.uns['cell_locations'].append(slice2.uns['cell_locations']).reset_index(drop=True)
Combine two slices into a unified 3D ST data with aligned spatial locations
[279]:
del slice1.raw
del slice2.raw
slice_3D = slice1.concatenate(
slice2,
batch_key=None,
uns_merge="unique",
index_unique=None
)
/home/jxiaoae/.conda/envs/SpatialScope/lib/python3.9/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],
.obsm[‘spatial’] and .uns[‘cell_locations’] are most important
[280]:
slice_3D.obsm['spatial'] = spatial_align
slice_3D.uns['cell_locations'] = cell_locations_align
[281]:
slice_3D.obsm['spatial']
[281]:
array([[ 1592.37822014, -1400.27868852, 0. ],
[ 2899.37822014, -322.27868852, 0. ],
[ -610.62177986, -202.27868852, 0. ],
...,
[-1084.26238373, -397.01811884, 48.43681451],
[-2336.26848367, 44.07793161, 48.43681451],
[ -830.24461369, 329.05004085, 48.43681451]])
[282]:
slice_3D.uns['cell_locations'].head()
[282]:
| x | y | spot_index | cell_index | cell_nums | slice_index | z | |
|---|---|---|---|---|---|---|---|
| 0 | 1552.878220 | -1381.778689 | slice1-AAACAGAGCGACTCCT-1 | slice1-AAACAGAGCGACTCCT-1_0 | 8 | slice1 | 0.0 |
| 1 | 1558.244417 | -1413.423055 | slice1-AAACAGAGCGACTCCT-1 | slice1-AAACAGAGCGACTCCT-1_1 | 8 | slice1 | 0.0 |
| 2 | 1585.297760 | -1439.451102 | slice1-AAACAGAGCGACTCCT-1 | slice1-AAACAGAGCGACTCCT-1_2 | 8 | slice1 | 0.0 |
| 3 | 1590.378220 | -1357.509458 | slice1-AAACAGAGCGACTCCT-1 | slice1-AAACAGAGCGACTCCT-1_3 | 8 | slice1 | 0.0 |
| 4 | 1582.372726 | -1382.547919 | slice1-AAACAGAGCGACTCCT-1 | slice1-AAACAGAGCGACTCCT-1_4 | 8 | slice1 | 0.0 |
[294]:
slice_3D.write('../output/cortex/sp_adata_ns.h5ad')
[353]:
slice_3D.X.max()
[353]:
7.3130717
Run SpatialScope (Step2: Cell type identification)
python ./src/Cell_Type_Identification.py --tissue cortex --out_dir ./output --ST_Data ./output/cortex/sp_adata_ns.h5ad --SC_Data ./Ckpts_scRefs/VISp/Ref_scRNA_VISp_qc2_2Kgenes.h5ad --cell_class_column cell_subclass --nu 15
[354]:
celltype_res = pd.read_csv('../output/cortex/CellTypeLabel_nu15.csv',index_col=0)
[355]:
celltype_res.head()
[355]:
| x | y | spot_index | cell_index | cell_nums | slice_index | z | discrete_label | discrete_label_ct | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 1552.878220 | -1381.778689 | slice1-AAACAGAGCGACTCCT-1 | slice1-AAACAGAGCGACTCCT-1_0 | 8 | slice1 | 0.0 | 3 | L2/3 IT |
| 1 | 1558.244417 | -1413.423055 | slice1-AAACAGAGCGACTCCT-1 | slice1-AAACAGAGCGACTCCT-1_1 | 8 | slice1 | 0.0 | 3 | L2/3 IT |
| 2 | 1585.297760 | -1439.451102 | slice1-AAACAGAGCGACTCCT-1 | slice1-AAACAGAGCGACTCCT-1_2 | 8 | slice1 | 0.0 | 3 | L2/3 IT |
| 3 | 1590.378220 | -1357.509458 | slice1-AAACAGAGCGACTCCT-1 | slice1-AAACAGAGCGACTCCT-1_3 | 8 | slice1 | 0.0 | 0 | Astro |
| 4 | 1582.372726 | -1382.547919 | slice1-AAACAGAGCGACTCCT-1 | slice1-AAACAGAGCGACTCCT-1_4 | 8 | slice1 | 0.0 | 3 | L2/3 IT |
[356]:
slice1_file = '../demo_data/Visium_MouseBrain_Cortex_section1.h5ad'
slice2_file = '../demo_data/Visium_MouseBrain_Cortex_section2.h5ad'
slice1 = sc.read(slice1_file)
slice2 = sc.read(slice2_file)
/home/jxiaoae/.conda/envs/SpatialScope/lib/python3.9/site-packages/anndata/_core/anndata.py:1830: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.
utils.warn_names_duplicates("var")
/home/jxiaoae/.conda/envs/SpatialScope/lib/python3.9/site-packages/anndata/_core/anndata.py:1830: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.
utils.warn_names_duplicates("var")
[357]:
slice1.uns['cell_locations']['discrete_label_ct'] = celltype_res.loc[celltype_res['slice_index']=='slice1','discrete_label_ct'].values
slice2.uns['cell_locations']['discrete_label_ct'] = celltype_res.loc[celltype_res['slice_index']=='slice2','discrete_label_ct'].values
[358]:
slice1.uns['cell_locations'].head()
[358]:
| x | y | spot_index | cell_index | cell_nums | discrete_label_ct | |
|---|---|---|---|---|---|---|
| 0 | 7910.500000 | 3181.500000 | AAACAGAGCGACTCCT-1 | AAACAGAGCGACTCCT-1_0 | 8 | L2/3 IT |
| 1 | 7915.866197 | 3149.855634 | AAACAGAGCGACTCCT-1 | AAACAGAGCGACTCCT-1_1 | 8 | L2/3 IT |
| 2 | 7942.919540 | 3123.827586 | AAACAGAGCGACTCCT-1 | AAACAGAGCGACTCCT-1_2 | 8 | L2/3 IT |
| 3 | 7948.000000 | 3205.769231 | AAACAGAGCGACTCCT-1 | AAACAGAGCGACTCCT-1_3 | 8 | Astro |
| 4 | 7939.994505 | 3180.730769 | AAACAGAGCGACTCCT-1 | AAACAGAGCGACTCCT-1_4 | 8 | L2/3 IT |
[359]:
color_dict = {'Astro': '#1f77b4',
'CR': '#aec7e8',
'Endo': '#ffbb78',
'L2/3 IT': '#ff7f0e',
'L4': '#2ca02c',
'L5 IT': '#bcbd22',
'L5 NP': '#d62728',
'L5 PT': '#ff9896',
'L6 CT': '#17becf',
'L6 IT': '#c5b0d5',
'L6b': '#8c564b',
'Lamp5': '#c49c94',
'Macrophage': '#e377c2',
'Meis2': '#f7b6d2',
'Oligo': '#7f7f7f',
'Peri': '#c7c7c7',
'Pvalb': '#98df8a',
'SMC': '#dbdb8d',
'Serpinf1': '#9467bd',
'Sncg': '#9edae5',
'Sst': '#82A8CE',
'VLMC': '#B87BCE',
'Vip': '#3FC241'}
[360]:
fig, ax = plt.subplots(1,1,figsize=(10, 8),dpi=130)
PlotVisiumCells(slice1,"discrete_label_ct",size=0.4,alpha_img=0.1,lw=0.7,palette=color_dict,ax=ax)
/import/home/jxiaoae/spatial/SpatialScope-Beta/demos/../src/utils.py:77: 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.
test = sc.AnnData(np.zeros(merged_df.shape), obs=merged_df)
/home/jxiaoae/.conda/envs/SpatialScope/lib/python3.9/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.
warnings.warn("Transforming to str index.", ImplicitModificationWarning)
[361]:
fig, ax = plt.subplots(1,1,figsize=(10, 8),dpi=130)
PlotVisiumCells(slice2,"discrete_label_ct",size=0.4,alpha_img=0.1,lw=0.7,palette=color_dict,ax=ax)
/import/home/jxiaoae/spatial/SpatialScope-Beta/demos/../src/utils.py:77: 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.
test = sc.AnnData(np.zeros(merged_df.shape), obs=merged_df)
/home/jxiaoae/.conda/envs/SpatialScope/lib/python3.9/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.
warnings.warn("Transforming to str index.", ImplicitModificationWarning)
[362]:
xyz1 = celltype_res.loc[celltype_res['slice_index']=='slice1',['x','y','z']].values
xyz2 = celltype_res.loc[celltype_res['slice_index']=='slice2',['x','y','z']].values
x1,y1,z1 = xyz1[:,0],xyz1[:,1],xyz1[:,2]
x2,y2,z2 = xyz2[:,0],xyz2[:,1],xyz2[:,2]
[363]:
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(10,5),dpi=200)
ax = fig.gca(projection='3d')# fig.add_subplot(projection='3d')
ax.get_proj = lambda: np.dot(Axes3D.get_proj(ax), np.diag([0.85, 1.6, 0.8, 1]))
ax.scatter(x1, y1, z1, c=slice1.uns['cell_locations']['discrete_label_ct'].map(color_dict), s=1.5, alpha=1)
ax.scatter(x2, y2, z2, c=slice2.uns['cell_locations']['discrete_label_ct'].map(color_dict), s=1.5, alpha=1)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_zlabel('$z$')
#ax.set_box_aspect([np.ptp(i) for i in data]) # equal aspect ratio
# add a legend
#handles = [Line2D([0], [0], marker='o', color='w', markerfacecolor=v, label=k, markersize=8) for k, v in color_dict.items()]
#ax.legend(title='color', handles=handles, bbox_to_anchor=(1.05, 1), loc='upper left')
ax.view_init(40,0)
#fig.colorbar(p, ax=ax)
plt.show()
/home/share/xiaojs/software/netMHCpan-4.1/Linux_x86_64/tmp/ipykernel_4891/1359269678.py:6: MatplotlibDeprecationWarning: Calling gca() with keyword arguments was deprecated in Matplotlib 3.4. Starting two minor releases later, gca() will take no keyword arguments. The gca() function should only be used to get the current axes, or if no axes exist, create new axes with default keyword arguments. To create a new axes with non-default arguments, use plt.axes() or plt.subplot().
ax = fig.gca(projection='3d')# fig.add_subplot(projection='3d')
[364]:
fig, ax = plt.subplots(1,1,figsize=(10, 8),dpi=130)
l_class = ['L4']#, 'L6 CT', 'L6b']
PlotVisiumCells(slice2,"discrete_label_ct",size=0.7,alpha_img=0.1,lw=0.7,subset=l_class,palette=color_dict,ax=ax)
/import/home/jxiaoae/spatial/SpatialScope-Beta/demos/../src/utils.py:77: 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.
test = sc.AnnData(np.zeros(merged_df.shape), obs=merged_df)
/home/jxiaoae/.conda/envs/SpatialScope/lib/python3.9/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.
warnings.warn("Transforming to str index.", ImplicitModificationWarning)
[365]:
fig, ax = plt.subplots(1,1,figsize=(10, 8),dpi=130)
l_class = ['L6 CT']#, 'L6 CT', 'L6b']
PlotVisiumCells(slice2,"discrete_label_ct",size=0.7,alpha_img=0.1,lw=0.7,subset=l_class,palette=color_dict,ax=ax)
/import/home/jxiaoae/spatial/SpatialScope-Beta/demos/../src/utils.py:77: 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.
test = sc.AnnData(np.zeros(merged_df.shape), obs=merged_df)
/home/jxiaoae/.conda/envs/SpatialScope/lib/python3.9/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.
warnings.warn("Transforming to str index.", ImplicitModificationWarning)
[ ]:
[ ]:
Learning the gene expression distribution of scRNA-seq reference using score-based model
We use four GPUs to train scRNA-seq reference in parallel.
python ./src/Train_scRef.py \
--ckpt_path ./Ckpts_scRefs/VISp \
--scRef ./Ckpts_scRefs/VISp/Ref_scRNA_VISp_qc2_2Kgenes.h5ad \
--cell_class_column cell_subclass \
--gpus 0,1,2,3
The checkpoints and sampled psuedo-cells will be saved in ./Ckpts_scRefs/VISp, e.g, model_5000.pt, model_5000.h5ad
Conveniently, we provided the pre-trained checkpoint (Ckpts_scRefs/VISp/model_5000.pt) in here, so you can skip this part
[ ]:
Run SpatialScope (Step3: Gene expression decomposition)
python ./src/Decomposition.py --tissue cortex --out_dir ./output --SC_Data ./Ckpts_scRefs/VISp/Ref_scRNA_VISp_qc2_2Kgenes.h5ad --cell_class_column cell_subclass --ckpt_path ./Ckpts_scRefs/VISp/model_5000.pt --spot_range 0,500 --gpu 2,3,4,8
python ./src/Decomposition.py --tissue cortex --out_dir ./output --SC_Data ./Ckpts_scRefs/VISp/Ref_scRNA_VISp_qc2_2Kgenes.h5ad --cell_class_column cell_subclass --ckpt_path ./Ckpts_scRefs/VISp/model_5000.pt --spot_range 500,1000 --gpu 2,3,4,8
python ./src/Decomposition.py --tissue cortex --out_dir ./output --SC_Data ./Ckpts_scRefs/VISp/Ref_scRNA_VISp_qc2_2Kgenes.h5ad --cell_class_column cell_subclass --ckpt_path ./Ckpts_scRefs/VISp/model_5000.pt --spot_range 1000,1715 --gpu 1,2,4,8
scRef and sampled pseudo-cells
[345]:
ad_sc = sc.read('../Ckpts_scRefs/VISp/Ref_scRNA_VISp_qc2_2Kgenes.h5ad')
sampled_cells = sc.read('../Ckpts_scRefs/VISp/model_5000.h5ad') # 2K cells sampled from the learned gene expression distribution of scRef
cell_type_key='cell_subclass'
sns.set_context('paper',font_scale=1.5)
adata_all = PlotSampledData(sampled_cells,ad_sc,cell_type_key,palette=color_dict)
/home/jxiaoae/.conda/envs/SpatialScope/lib/python3.9/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],
decomposed single-cell resolution ST data
[383]:
ad_scsts_list = []
ad_scsts_list.append(sc.read('../output/cortex/generated_cells_spot0_500.h5ad'))
ad_scsts_list.append(sc.read('../output/cortex/generated_cells_spot500_1000.h5ad'))
ad_scsts_list.append(sc.read('../output/cortex/generated_cells_spot1000_1715.h5ad'))
ad_scsts = ad_scsts_list[0].concatenate(
ad_scsts_list[1:],
batch_key="_",
uns_merge="unique",
index_unique=None
)
ad_scst_slice2 = ad_scsts[ad_scsts.obs.slice_index=='slice2']
/home/jxiaoae/.conda/envs/SpatialScope/lib/python3.9/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],
[401]:
slice2_file = '../demo_data/Visium_MouceBrain_Cortex_section2.h5ad'
slice2 = sc.read(slice2_file)
slice2.X = np.exp(slice2.X.A)-1
/home/jxiaoae/.conda/envs/SpatialScope/lib/python3.9/site-packages/anndata/_core/anndata.py:1830: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.
utils.warn_names_duplicates("var")
[402]:
slice2.X.max(),ad_scst_slice2.X.max()
[402]:
(3490.1306, ArrayView(14820.539, dtype=float32))
[405]:
gene = 'Rgs4' #[L2/3: 'Cux2','Otof', 'Lamp5', 'Nov' | L4: 'Rgs4', 'Rorb', 'Otof' | L5 IT: Ptn, Sulf2, Fezf2 | L6: Pcp4, Osr1 | Oligo: Sox10, Aqp4 | VLMC: Ogn, Lum, and Dcn]
sns.set_context('paper',font_scale=2)
fig, ax = plt.subplots(2,1,figsize=(12, 12),dpi=100)
plt.rcParams["legend.markerscale"] = 2
PlotVisiumGene(slice2,gene,size=1.7,alpha_img=0.05,perc=0,ax=ax[0], palette = sns.cubehelix_palette(light = 1, as_cmap = True))
PlotVisiumGene(ad_scst_slice2,gene,size=2.0,alpha_img=0.05,perc=0,ax=ax[1], palette = sns.cubehelix_palette(light = 1, as_cmap = True))
plt.tight_layout()
/home/jxiaoae/.conda/envs/SpatialScope/lib/python3.9/site-packages/anndata/_core/anndata.py:1830: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.
utils.warn_names_duplicates("var")
[407]:
gene = 'Pcp4' #[L2/3: 'Cux2','Otof', 'Lamp5', 'Nov' | L4: 'Rgs4', 'Rorb', 'Otof' | L5 IT: Ptn, Sulf2, Fezf2 | L6: Pcp4, Osr1 | Oligo: Sox10, Aqp4 | VLMC: Ogn, Lum, and Dcn]
sns.set_context('paper',font_scale=2)
fig, ax = plt.subplots(2,1,figsize=(12, 12),dpi=100)
plt.rcParams["legend.markerscale"] = 2
PlotVisiumGene(slice2,gene,size=1.7,alpha_img=0.05,perc=0,ax=ax[0], palette = sns.cubehelix_palette(light = 1, as_cmap = True))
PlotVisiumGene(ad_scst_slice2,gene,size=2.0,alpha_img=0.05,perc=0,ax=ax[1], palette = sns.cubehelix_palette(light = 1, as_cmap = True))
plt.tight_layout()
/home/jxiaoae/.conda/envs/SpatialScope/lib/python3.9/site-packages/anndata/_core/anndata.py:1830: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.
utils.warn_names_duplicates("var")
[408]:
gene = 'Sox10' #[L2/3: 'Cux2','Otof', 'Lamp5', 'Nov' | L4: 'Rgs4', 'Rorb', 'Otof' | L5 IT: Ptn, Sulf2, Fezf2 | L6: Pcp4, Osr1 | Oligo: Sox10, Aqp4 | VLMC: Ogn, Lum, and Dcn]
sns.set_context('paper',font_scale=2)
fig, ax = plt.subplots(2,1,figsize=(12, 12),dpi=100)
plt.rcParams["legend.markerscale"] = 2
PlotVisiumGene(slice2,gene,size=1.7,alpha_img=0.05,perc=0,ax=ax[0], palette = sns.cubehelix_palette(light = 1, as_cmap = True))
PlotVisiumGene(ad_scst_slice2,gene,size=2.0,alpha_img=0.05,perc=0,ax=ax[1], palette = sns.cubehelix_palette(light = 1, as_cmap = True))
plt.tight_layout()
/home/jxiaoae/.conda/envs/SpatialScope/lib/python3.9/site-packages/anndata/_core/anndata.py:1830: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.
utils.warn_names_duplicates("var")
[ ]: