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)

用中性笔做狙击枪_n在python中怎么用_用中国传统色打开春耕画卷

我们还可以把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()

n在python中怎么用_用中性笔做狙击枪_用中国传统色打开春耕画卷

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