1. 前言
在上一篇推送中,我们给大家介绍了用于整合不同条件(样本条件),技术(测序手段),发育阶段(针对发育数据集)的空间转录组数据集的算法-STAligner,今天我们在这里给出如何使用STAligner的教程。
目前STAligner已经被整合进了Omicverse这个python包,Omicverse中关于STAligner教程的链接如下:
STAligner官方教程链接如下:
2. 环境配置
conda create -n omicverse python=3.10
conda activate omicverse
conda install mamba -c conda-forge
mamba install jax jaxlib -c conda-forge
pip3 install torch torchvision torchaudio
pip3 install torch_geometric
pip3 install pyg_lib torch_scatter torch_sparse torch_cluster torch_spline_conv -f https://data.pyg.org/whl/torch-2.3.0+cu121.html
pip3 install -U omicvers
3. 载入包和数据
from scipy.sparse import csr_matrix
import omicverse as ov
import scanpy as sc
import anndata as ad
import pandas as pd
import os
ov.utils.ov_plot_set(
我们使用官方提供的数据,这是一个使用stereo-seq进行测序,和一个使用slide-seq进行测序的小鼠大脑嗅球区域的空间转录组切片,这些数据存储在谷歌云盘中,各位可以直接下载。()
对于各位自己的课题,只要确保数据格式为.h5ad,可以被读取为anndata数据结构的数据即可。
Batch_list = []
adj_list = []
section_ids = ['Slide-seqV2_MoB', 'Stereo-seq_MoB']
print(section_ids)
pathway = './scRNA-seq/scripts/STAligner'
我们首先需要对数据进行预处理。对于STAligner算法,我们需要注意两点,第一是高变基因数量的选择,第二是用于构建空间图的空间半径的选择,这两个参数会直接影响模型运行的效果。
首先我们谈高变基因数量的选择:即sc.pp.highly_variable_genes(adata, flavor=”seurat_v3″, n_top_genes=10000)函数中的n_top_genes参数,对于 STAligner,它首先计算单个样本的高变基因,然后再把不同的样本拼接起来,选择它们共有的高变基因。因此,高边基因的数量不应该选择得太低,否则,在样本量比较多的情况下,将没有足够多的特征来给staligner进行训练,进而影响下游分析的效果。
在我自己的课题中,我有14张Stereo-seq切片,每个样本选择前8000个高变基因,14个样本拼接完只剩下大约1400个基因。(当然也可以先拼接再选高变基因,但是这种情况内存占用峰值会比较高,非常容易爆内存)
然后是空间半径的选择,即ov.space.Cal_Spatial_Net(adata, rad_cutoff=50)函数中的rad_cutoff参数。空间半径是指在空间转录组切片中,以每个spots为圆心,在空间半径范围内的所有其他spots。在构建空间图时,将这些spots与圆心spots连上一条边。
在实际的使用中,代码会输出每张空间转录组学切片平均每个spots有多少个临近spots与之相连(举个官网教程的例子:The graph contains 144318 edges, 19109 cells. 7.5524 neighbors per cell on average.),在使用过程中,大家需要根据数据的不同调整rad_cutoff参数,确保每个spots平均有5-10个临近spots与之相连即可。
for section_id in section_ids:
print(section_id)
adata = sc.read_h5ad(os.path.join(pathway,section_id+".h5ad"))
# check whether the adata.X is sparse matrix
if isinstance(adata.X, pd.DataFrame):
adata.X = csr_matrix(adata.X)
else:
pass
adata.var_names_make_unique(join="++")
# make spot name unique
adata.obs_names = [x+'_'+section_id for x in adata.obs_names]
# Constructing the spatial network
ov.space.Cal_Spatial_Net(adata, rad_cutoff=50) # the spatial network are saved in adata.uns[‘adj’]
# Normalization
sc.pp.highly_variable_genes(adata, flavor="seurat_v3", n_top_genes=10000)
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
adata = adata[:, adata.var['highly_variable']]
adj_list.append(adata.uns['adj'])
Batch_list.append(adata)
接下来我们需要把每张空间转录组学切片,给拼接成一个大的anndata数据结构,用于后续的模型训练。
adata_concat = ad.concat(Batch_list, label="slice_name", keys=section_ids)
adata_concat.obs["batch_name"] = adata_concat.obs["slice_name"].astype('category')
print('adata_concat.shape: ', adata_concat.shape
然后我们就可以训练STAligner模型了,omicverse封装了训练过程的所有代码,大家直接运行下面三行代码即可运行STAligner。
%%time
# iter_comb is used to specify the order of integration. For example, (0, 1) means slice 0 will be algined with slice 1 as reference.
iter_comb = [(i, i + 1) for i in range(len(section_ids) - 1)]
# Here, to reduce GPU memory usage, each slice is considered as a subgraph for training.
STAligner_obj = ov.space.pySTAligner(adata_concat, verbose=True, knn_neigh = 100, n_epochs = 600, iter_comb = iter_comb,
batch_key = 'batch_name', key_added='STAligner', Batch_list = Batch_list)
STAligner_obj.train()
adata = STAligner_obj.predicted()
4. 下游分析
STAligner模型的输出是空间中每个spots的embedding矩阵,我们可以利用这个矩阵去进行下游聚类分析,并且确定跨切片的空间区域结构。
我们使用 UMAP 图来可视化整合了两张空间转录组学切片的效果。左图以样本类型为基准进行了着色,清晰地展示了不同切片的 spots 被很好地混合在一起,消除了分辨率的影响(slide-seq 分辨率约为 10 微米,stereo-seq 分辨率约为 25 微米,在 bin50 情况下)。右图基于STAligner生成的embedding矩阵,利用 Leiden 算法得到的聚类结果。
sc.pp.neighbors(adata, use_rep='STAligner', random_state=666)
ov.utils.cluster(adata,use_rep='STAligner',method='leiden',resolution=0.4)
sc.tl.umap(adata, random_state=666)
sc.pl.umap(adata, color=['batch_name',"leiden"],wspace=0.5)
我们还可以把leiden算法得到的聚类结果投影回空间坐标上,这样我们就可以很直观的看到STAligner识别的空间区域。
import matplotlib.pyplot as plt
spot_size = 50
title_size = 15
fig, ax = plt.subplots(1, 2, figsize=(6, 3), gridspec_kw={'wspace': 0.05, 'hspace': 0.2})
_sc_0 = sc.pl.spatial(adata[adata.obs['batch_name'] == 'Slide-seqV2_MoB'], img_key=None, color=['leiden'], title=['Slide-seqV2'],
legend_fontsize=10, show=False, ax=ax[0], frameon=False, spot_size=spot_size, legend_loc=None)
_sc_0[0].set_title('Slide-seqV2', size=title_size)
_sc_1 = sc.pl.spatial(adata[adata.obs['batch_name'] == 'Stereo-seq_MoB'], img_key=None, color=['leiden'], title=['Stereo-seq'],
legend_fontsize=10, show=False, ax=ax[1], frameon=False, spot_size=spot_size)
_sc_1[0].set_title('Stereo-seq',size=title_size)
_sc_1[0].invert_yaxis()
plt.show()
5. 总结
在上一次的推文中,我们向大家介绍了用于整合多张空间转录组学切片的 STAligner 算法。正好,Omicverse 包含了 STAligner 算法,因此在这篇笔记中,我们给大家介绍了如何在 Omicverse 中使用 STAligner 算法。Omicverse是一个整合了多种(大于30种)空间转录组,scRNA-seq和Bulk RNA-seq分析的工具包,为我们进行相关生物信息学分析提供了巨大的帮助。希望通过这篇文章,能够启发大家探索 Omicverse 的其他功能,我们也将继续分享这些笔记,与大家一同学习。
6. 参考资料
Zhou X, Dong K, Zhang S. Integrating spatial transcriptomics data across different conditions, technologies and developmental stages[J]. Nature Computational Science, 2023, 3(10): 894-906.
限时特惠:本站每日持续更新海内外内部创业教程,一年会员只需88元,全站资源免费下载点击查看详情。
站长微信:nnxmw123