日常吐槽。
一、写在前面
stereopy日常出bug(github issue里得有一半的问题是我提的，当然也有可能是因为我菜)，stereopy自己生成的anndata自己不能计算空间共现关系，还是靠squidpy才能计算。另外还要一些函数一开并行计算就报错，这里留一些stereopy对象的转换/输出方式，以备大家不(sui)时之(pao)需(lu)。更多空转教程可见：空间转录组学习手册合辑
二、stereopy相关数据读写
Stereopy是一个基础且全面的空间转录组分析工具，其支持空间高变基因计算、RNA Velocity、批次整合、SingleR注释、细胞通讯、基因调控网络、轨迹推断等基础/进阶分析。此外，针对多样本的数据，Stereopy还支持3D细胞通讯、3D轨迹推断、3D基因调控网络等分析。学习下面的内容前需要大家先掌握python知识：生信Python速查手册
文件读取、输出、转换：
1.1 GEM与GEF格式文件
一共有GEM、GEF两种格式的数据。GEM文件包括：切片中x、y轴空间坐标信息，MIDCount为基因表达量，ExonCount为由SAW(Version >= 5.1.3)中的spatial_RNA_visualization_v5生成的外显子表达量。
数据格式包含以下信息：
GEF分为Square Bin与Cell Bin两种，前者的前缀一般为.raw.gef/.gef/.tissue.gef，后者的前缀一般为.cellbin.gef。(Stereo-seq技术分辨率可达纳米级别bin1的纳米孔半径约为250nm。通过bin_size参数，可以将一定范围内的纳米孔数据合并为一个bin unit(相当于一个分类单元)。例如，当bin_size设置为20时，bin unit的边长大约为10/14 μm(取决于纳米孔中心的距离：500nm或715nm)。在bin_size确定后，可自动生成StereoExpData用于下游分析。由于计算性能、生物学意义等问题，在分析过程中可能需要跳回“分bin”的步骤，优化参数后继续处理。)
此外stereopy还可以读取来自Anndata、Scanpy、Seurat的h5ad格式文件。
1.2 数据读取
# 导入包
import stereo as st
import warnings
warnings.filterwarnings('ignore')1.2.1 GEM文件读取
# gem文件读取
data_path = 'Demo_MouseBrain/SS200000135TL_D1.cellbin.gem'
data = st.io.read_gem(
file_path=data_path,
sep='\t',
bin_type='cell_bins',
is_sparse=True,
)
data
# [2023-09-19 06:31:29][Stereo][3273491][MainThread][139900354834816][reader][84][INFO]: the martrix has 57133 cells, and 22406 genes.
# StereoExpData object with n_cells X n_genes = 57133 X 22406
# bin_type: cell_bins
# offset_x = 3206# offset_y = 6174
# cells: ['cell_name', 'cell_point']
# genes: ['gene_name']1.2.2 GEF文件读取
# suqre bin读取
data_path = 'Demo_MouseBrain/SS200000135TL_D1.tissue.gef'
data = st.io.read_gef(
file_path=data_path,
bin_type='bins',
bin_size=100,
is_sparse=True,
)
data
# [2023-09-19 06:33:03][Stereo][3273491][MainThread][139900354834816][reader][1001][INFO]: read_gef begin ...
# [2023-09-19 06:33:23][Stereo][3273491][MainThread][139900354834816][reader][1088][INFO]: the matrix has 9124 cells, and 24302 genes.
# [2023-09-19 06:33:25][Stereo][3273491][MainThread][139900354834816][reader][1096][INFO]: read_gef end.
# StereoExpData object with n_cells X n_genes = 9124 X 24302
# bin_type: bins
# bin_size: 100
# offset_x = 0
# offset_y = 0
# cells: ['cell_name']
# genes: ['gene_name']# cell bin读取
data_path = './Demo_MouseBrain/SS200000135TL_D1.cellbin.gef'
data = st.io.read_gef(
file_path=data_path,
is_sparse=True,
bin_type='cell_bins',
)
# [2023-09-19 06:33:54][Stereo][3273491][MainThread][139900354834816][reader][1001][INFO]: read_gef begin ...
# [2023-09-19 06:33:58][Stereo][3273491][MainThread][139900354834816][reader][1038][INFO]: the matrix has 57133 cells, and 24670 genes.
# [2023-09-19 06:33:58][Stereo][3273491][MainThread][139900354834816][reader][1161][INFO]: This is GEF file which contains cell bin infomation.
# [2023-09-19 06:33:58][Stereo][3273491][MainThread][139900354834816][reader][1162][INFO]: bin_type: cell_bins
# [2023-09-19 06:33:58][Stereo][3273491][MainThread][139900354834816][reader][1168][INFO]: Number of cells: 57133
# [2023-09-19 06:33:58][Stereo][3273491][MainThread][139900354834816][reader][1171][INFO]: Number of gene: 24670
# [2023-09-19 06:33:58][Stereo][3273491][MainThread][139900354834816][reader][1174][INFO]: Resolution: 500
# [2023-09-19 06:33:58][Stereo][3273491][MainThread][139900354834816][reader][1177][INFO]: offsetX: 0
# [2023-09-19 06:33:58][Stereo][3273491][MainThread][139900354834816][reader][1180][INFO]: offsetY: 0
# [2023-09-19 06:33:58][Stereo][3273491][MainThread][139900354834816][reader][1183][INFO]: Average number of genes: 223.460693359375
# [2023-09-19 06:33:58][Stereo][3273491][MainThread][139900354834816][reader][1186][INFO]: Maximum number of genes: 1046
# [2023-09-19 06:33:58][Stereo][3273491][MainThread][139900354834816][reader][1189][INFO]: Average expression: 399.17034912109375
# [2023-09-19 06:33:58][Stereo][3273491][MainThread][139900354834816][reader][1192][INFO]: Maximum expression: 2337
# [2023-09-19 06:33:58][Stereo][3273491][MainThread][139900354834816][reader][1096][INFO]: read_gef end.1.2.3 Stereopy来源的h5ad文件读取
data_path = './Demo_MouseBrain/StereoExpData_stereopy_mykey.h5ad'
data = st.io.read_stereo_h5ad(
file_path=data_path,
use_raw=True,
use_result=True,
)
data
# StereoExpData object with n_cells X n_genes = 35998 X 3000
# bin_type: bins
# bin_size: 50
# offset_x = None
# offset_y = None
# cells: ['cell_name', 'total_counts', 'pct_counts_mt', 'n_genes_by_counts', 'leiden', 'louvain']
# genes: ['gene_name', 'n_cells', 'n_counts']
# key_record: {'cluster': ['leiden', 'louvain'], 'gene_exp_cluster': ['gene_exp_leiden', 'gene_exp_louvain']}1.2.4 Anndata来源的h5ad文件读取
ann_h5ad = './Demo_MouseBrain/anndata_output.h5ad'
data = st.io.read_ann_h5ad(
file_path=ann_h5ad,
spatial_key=None,
)
data
# StereoExpData object with n_cells X n_genes = 35998 X 24302
# bin_type: bins
# bin_size: 50
# offset_x = None
# offset_y = None
# cells: ['cell_name']
# genes: ['gene_name']# Anndata来源的h5ad文件也可用于scanpy的分析
import scanpy as sc
# 直接利用scanpy读取：
data = sc.read_h5ad('Demo_MouseBrain/anndata_output.h5ad')
data
# AnnData object with n_obs × n_vars = 35998 × 24302
# obs: 'orig.ident', 'x', 'y'
# uns: 'bin_size', 'bin_type', 'raw_cellname', 'raw_counts', 'raw_genename', 'resolution', 'sn'
# obsm: 'spatial'1.3 文件输出
虽然输出的文件均为h5ad文件，但不同格式输出的h5ad文件具有差异，需要区别对待，输出对应软件需要的版本。
1.3.1 用于stereopy读取的h5ad输出¶
# 读取GEF文件data_path = './Demo_MouseBrain/SS200000135TL_D1.tissue.gef'data = st.io.read_gef(file_path=data_path, bin_size=50)# 预处理data.tl.cal_qc()data.tl.raw_checkpoint()data.tl.sctransform(res_key='sctransform', inplace=True)# 聚类data.tl.pca(use_highly_genes=False, n_pcs=30, res_key='pca')data.tl.neighbors(pca_res_key='pca', n_pcs=30, res_key='neighbors')data.tl.umap(pca_res_key='pca', neighbors_res_key='neighbors', res_key='umap')data.tl.leiden(neighbors_res_key='neighbors', res_key='leiden')data.tl.louvain(neighbors_res_key='neighbors', res_key='louvain')# data.tl.key_record为运行函数时自动生成的字典，收录的为res_key参数返回的内容，我们在后面正式的数据分析教程中会介绍print(data.tl.key_record)# [2023-09-19 08:09:11][Stereo][3273491][MainThread][139900354834816][st_pipeline][40][INFO]: umap end, consume time 24.3978s.# [2023-09-19 08:09:11][Stereo][3273491][MainThread][139900354834816][st_pipeline][37][INFO]: start to run leiden...# [2023-09-19 08:09:20][Stereo][3273491][MainThread][139900354834816][st_pipeline][40][INFO]: leiden end, consume time 9.2709s.# [2023-09-19 08:09:21][Stereo][3273491][MainThread][139900354834816][st_pipeline][37][INFO]: start to run louvain...# [2023-09-19 08:09:21][Stereo][3273491][MainThread][139900354834816][_louvain][100][INFO]: using the "louvain" package of Traag (2017)# [2023-09-19 08:09:35][Stereo][3273491][MainThread][139900354834816][st_pipeline][40][INFO]: louvain end, consume time 14.8895s.# {'hvg': [], 'pca': ['pca'], 'neighbors': ['neighbors'], 'umap': ['umap'], 'cluster': ['leiden', 'louvain'], 'marker_genes': [], 'sct': ['sctransform'], 'gene_exp_cluster': ['gene_exp_leiden', 'gene_exp_louvain']}
# 将StereoExpData对象写为h5ad, 如果key_record = None, 则会自动采用data.tl.key_recordst.io.write_h5ad(data,use_raw=True,use_result=True,key_record=None,output='./Demo_MouseBrain/StereoExpData_stereopy.h5ad')# 原gef大约600MB，输出后的h5ad约为3GB
# 当然你也可以建立一个字典保存指定的key_recordoutkey_record = {'cluster':['leiden','louvain'],}st.io.write_h5ad(data,use_raw=True,use_result=True,key_record=outkey_record,output='./Demo_MouseBrain/StereoExpData_stereopy_mykey.h5ad',)# 指定key_record后的输出文件大小为1.1GB，需要节省磁盘空间的同学可以尝试这种方法
1.3.2 用于Anndata读取的h5ad输出
# 读取gef文件：data_path = './Demo_MouseBrain/SS200000135TL_D1.tissue.gef'data = st.io.read_gef(file_path=data_path, bin_size=50)data.tl.raw_checkpoint()# 转换输出adata = st.io.stereo_to_anndata(data,flavor='seurat',output='Demo_MouseBrain/anndata_output.h5ad')# [2023-09-19 07:07:41][Stereo][3273491][MainThread][139900354834816][reader][1001][INFO]: read_gef begin ...# [2023-09-19 07:08:00][Stereo][3273491][MainThread][139900354834816][reader][1088][INFO]: the matrix has 35998 cells, and 24302 genes.# [2023-09-19 07:08:02][Stereo][3273491][MainThread][139900354834816][reader][1096][INFO]: read_gef end.# [2023-09-19 07:08:03][Stereo][3273491][MainThread][139900354834816][reader][756][INFO]: Adding sample in adata.obs['orig.ident'].# [2023-09-19 07:08:03][Stereo][3273491][MainThread][139900354834816][reader][759][INFO]: Adding data.position as adata.obsm['spatial'] .# [2023-09-19 07:08:03][Stereo][3273491][MainThread][139900354834816][reader][764][INFO]: Adding data.position as adata.obs['x'] and adata.obs['y'] .# [2023-09-19 07:08:03][Stereo][3273491][MainThread][139900354834816][reader][853][INFO]: Adding data.tl.raw.exp_matrix as adata.uns['raw_counts'] .# [2023-09-19 07:08:03][Stereo][3273491][MainThread][139900354834816][reader][884][INFO]: Rename QC info.# [2023-09-19 07:08:03][Stereo][3273491][MainThread][139900354834816][reader][900][INFO]: Finished conversion to anndata.# [2023-09-19 07:08:06][Stereo][3273491][MainThread][139900354834816][reader][904][INFO]: Finished output to Demo_MouseBrain/anndata_output.h5ad
1.3.3 用于Seurat读取的h5ad输出
Seurat对象可谓是单细胞分析界的硬通货，一定要学会转化。
# 读取数据data_path = './Demo_MouseBrain/SS200000135TL_D1.tissue.gef'data = st.io.read_gef(file_path=data_path, bin_size=50)# 预处理一下data.tl.cal_qc()data.tl.raw_checkpoint()# 存为Seurat友好型的h5adadata = st.io.stereo_to_anndata(data,flavor='seurat',output='Demo_MouseBrain/seurat_out.h5ad')# [2023-09-19 08:32:39][Stereo][3273491][MainThread][139900354834816][reader][1001][INFO]: read_gef begin ...# [2023-09-19 08:32:58][Stereo][3273491][MainThread][139900354834816][reader][1088][INFO]: the matrix has 35998 cells, and 24302 genes.# [2023-09-19 08:33:00][Stereo][3273491][MainThread][139900354834816][reader][1096][INFO]: read_gef end.# [2023-09-19 08:33:00][Stereo][3273491][MainThread][139900354834816][st_pipeline][37][INFO]: start to run cal_qc...# [2023-09-19 08:33:00][Stereo][3273491][MainThread][139900354834816][st_pipeline][40][INFO]: cal_qc end, consume time 0.7293s.# [2023-09-19 08:33:01][Stereo][3273491][MainThread][139900354834816][reader][756][INFO]: Adding sample in adata.obs['orig.ident'].# [2023-09-19 08:33:01][Stereo][3273491][MainThread][139900354834816][reader][759][INFO]: Adding data.position as adata.obsm['spatial'] .# [2023-09-19 08:33:01][Stereo][3273491][MainThread][139900354834816][reader][764][INFO]: Adding data.position as adata.obs['x'] and adata.obs['y'] .# [2023-09-19 08:33:01][Stereo][3273491][MainThread][139900354834816][reader][853][INFO]: Adding data.tl.raw.exp_matrix as adata.uns['raw_counts'] .# [2023-09-19 08:33:01][Stereo][3273491][MainThread][139900354834816][reader][884][INFO]: Rename QC info.# [2023-09-19 08:33:01][Stereo][3273491][MainThread][139900354834816][reader][900][INFO]: Finished conversion to anndata.# [2023-09-19 08:33:04][Stereo][3273491][MainThread][139900354834816][reader][904][INFO]: Finished output to Demo_MouseBrain/seurat_out.h5ad
后续h5ad在Seurat中的读取可以参考：单细胞对象(数据格式)转换大全|2. h5ad转Seuratobj。
当然，如果你更熟悉rds文件，也可以在shell中做如下转换：
!/usr/bin/Rscript ~/stereopy/docs/source/_static/annh5ad2rds.R --infile 'Demo_MouseBrain/seurat_out.h5ad' --outfile 'Demo_MouseBrain/seurat_out.rds'1.3.4 GEF输出
创建新的GEF用于输出：
# 数据读取data_path = './Demo_MouseBrain/SS200000135TL_D1.tissue.gef'data = st.io.read_gef(file_path=data_path, bin_size=50)# 选择对应基因取子集data.tl.filter_genes(gene_list=['H2al2a','Gm6135'], inplace=True)# 仅保存过滤后的结果st.io.write_mid_gef(data=data,output='./Demo_MouseBrain/my_filter.filtered.gef')# [2023-09-19 08:18:02][Stereo][3273491][MainThread][139900354834816][reader][1001][INFO]: read_gef begin ...# [2023-09-19 08:18:22][Stereo][3273491][MainThread][139900354834816][reader][1088][INFO]: the matrix has 35998 cells, and 24302 genes.# [2023-09-19 08:18:24][Stereo][3273491][MainThread][139900354834816][reader][1096][INFO]: read_gef end.# [2023-09-19 08:18:24][Stereo][3273491][MainThread][139900354834816][st_pipeline][37][INFO]: start to run filter_genes...# [2023-09-19 08:18:24][Stereo][3273491][MainThread][139900354834816][st_pipeline][40][INFO]: filter_genes end, consume time 0.5481s.# [2023-09-19 08:18:24][Stereo][3273491][MainThread][139900354834816][writer][287][INFO]: The output standard gef file only contains one expression matrix with mid count.Please make sure the expression matrix of StereoExpData object is mid count without normaliztion.
# 读取GEF文件data_path = './Demo_MouseBrain/SS200000135TL_D1.tissue.gef'data = st.io.read_gef(file_path=data_path, bin_size=50)# 预处理data.tl.cal_qc()data.tl.raw_checkpoint()data.tl.sctransform(res_key='sctransform', inplace=True)# 分群data.tl.pca(use_highly_genes=False, n_pcs=30, res_key='pca')data.tl.neighbors(pca_res_key='pca', n_pcs=30, res_key='neighbors')data.tl.umap(pca_res_key='pca', neighbors_res_key='neighbors', res_key='umap')data.tl.leiden(neighbors_res_key='neighbors', res_key='leiden')# 更新本地已存在的gefexist_path = './Demo_MouseBrain/my_filter.filtered.gef'st.io.update_gef(data=data,gef_file=exist_path,cluster_res_key='leiden',)
保存于已存在的GEF文件:
1.4 以AnnData格式操作
AnnData对象在Python亦是单细胞/空转分析的硬通货，Stereo中对象与其的相互转换与交互值得我们单独探究一番：
1.4.1 AnnData来源h5ad读取
import stereo as stdata = st.io.read_h5ad('Demo_MouseBrain/anndata_output.h5ad')# 此时读入的便是一个AnnData对象：data# AnnData object with n_obs × n_vars = 35998 × 24302# obs: 'orig.ident', 'x', 'y'# uns: 'bin_size', 'bin_type', 'raw_cellname', 'raw_counts', 'raw_genename', 'resolution', 'sn'# obsm: 'spatial'
# data._ann_data可以完全作为一个AnnData对象处理，而不是StereoExpData对象data._ann_data# AnnData object with n_obs × n_vars = 35998 × 24302# obs: 'orig.ident', 'x', 'y'# uns: 'bin_size', 'bin_type', 'raw_cellname', 'raw_counts', 'raw_genename', 'resolution', 'sn'# obsm: 'spatial'
1.4.2 各数据存储位置
# 表达矩阵：data._ann_data.X# <35998x24302 sparse matrix of type '<class 'numpy.float64'>'# with 41464124 stored elements in Compressed Sparse Row format>
# 空间坐标信息：data._ann_data.uns# OverloadedDict, wrapping:# {'bin_size': 50, 'bin_type': 'bins', 'raw_cellname': array(['24481313596250', '26199300516500', '26843545608050', ...,# '48533130461250', '38439957317800', '53687091208300'], dtype=object), 'raw_counts': <35998x24302 sparse matrix of type '<class 'numpy.uint32'>'# with 41464124 stored elements in Compressed Sparse Row format>, 'raw_genename': array(['0610005C13Rik', '0610006L08Rik', '0610009B22Rik', ..., 'mt-Nd4l',# 'mt-Nd5', 'mt-Nd6'], dtype=object), 'resolution': 500, 'sn': batch sn# 0 -1 SS200000135TL_D1}# With overloaded keys:# ['neighbors'].
# 细胞及相关注释信息：data._ann_data.obs
# 基因及相关注释信息：data._ann_data.var
1.4.3 数据处理
AnnData可直接参与大部分用于StereoExpData的计算
# 例如以下这些：data.tl.cal_qc()data.tl.raw_checkpoint()data.tl.normalize_total(target_sum=1e4)data.tl.log1p()data.tl.highly_variable_genes(min_mean=0.0125, max_mean=3, min_disp=0.5, res_key='highly_variable_genes', n_top_genes=None)data.tl.pca(use_highly_genes=True, hvg_res_key='highly_variable_genes', n_pcs=20, res_key='pca_test', svd_solver='arpack')data.tl.neighbors(pca_res_key='pca_test', n_pcs=30, res_key='neighbors_test', n_jobs=8)data.tl.umap(pca_res_key='pca_test', neighbors_res_key='neighbors_test', res_key='umap_test', init_pos='spectral')data.tl.leiden(neighbors_res_key='neighbors_test', res_key='leiden_test')
# 计算得到的结果也会自动保存于data._ann_data中data._ann_data# AnnData object with n_obs × n_vars = 35998 × 24302# obs: 'orig.ident', 'x', 'y', 'total_counts', 'n_genes_by_counts', 'pct_counts_mt', 'leiden_test'# var: 'n_cells', 'n_counts', 'mean_umi', 'means', 'dispersions', 'dispersions_norm', 'highly_variable'# uns: 'bin_size', 'bin_type', 'raw_cellname', 'raw_counts', 'raw_genename', 'resolution', 'sn', 'highly_variable_genes', 'pca_test', 'neighbors_test', 'umap_test', 'leiden_test', 'gene_exp_leiden_test'# obsm: 'spatial', 'X_pca_test', 'X_umap_test'# obsp: 'neighbors_test_connectivities', 'neighbors_test_distances'
# 例如这里多出的细胞信息：data._ann_data.obs
# 可视化函数一样可以利用：data.plt.umap(res_key='umap_test', cluster_key='leiden_test')
1.4.4 数据输出
用于stereopy分析的AnnData对象也可以存为h5ad，
data._ann_data.write_h5ad('./test_result/SS200000135TL_D1.stereo.h5ad')更多空间转录组分析技巧可见：空间转录组学习手册合辑