你是否为各种各样的空间转录组数据类型而眼花缭乱？你是否想要复现文章中的分析而倒在了数据读入？你是否希望学习空间转录组数据分析却又无从下手？
今天，我们就和大家一起来学习多种类型空间转录组数据的读取，迈开空间转录组数据分析的第一步。
当前常见的基于空间条形码和NGS测序的空间组转录组技术包括10x Visium、华大基因的Stereo-seq和Broad研究所的Slide-seq等，每种空间转录组技术都有特有的基因表达矩阵文件，格式各异。下面，我们就依次来介绍上述3种空间转录组技术的数据读取。
10x Visium
10x Visium是目前使用最广泛的空间组学技术，Visium数据经SpaceRanger处理，生成类似如下结构的输出文件：
在Seurat中，使用
Load10X_Spatial()函数加载10x Visium数据，该函数读入SpaceRanger的输出，并返回一个Seurat对象，该对象既包含基因表达数据，也包含组织切片的相关图像。例：
library(Seurat)
library(ggplot2)
obj <- Load10X_Spatial('/path/spaceranger_output/',
filename = "filtered_feature_bc_matrix.h5",
assay = "Spatial",
slice = "slice1",
filter.matrix = TRUE
)
obj
plot <- SpatialFeaturePlot(obj, features = "nCount_Spatial") + theme(legend.position = "right")
plot
An object of class Seurat
32285 features across 2797 samples within 1 assay
Active assay: Spatial (32285 features, 0 variable features)
1 image present: slice1
这样我们就成功导入了1个10x Visium的空间转录组数据。
如果你习惯使用Python进行数据处理，我们也可以用scanpy来读入Visium数据。scanpy使用scanpy.read_visium()函数导入Visium数据为anndata格式。
例：
import scanpy as sc
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
adata = sc.read_visium('/path/spaceranger_output/',
count_file='filtered_feature_bc_matrix.h5',
load_images=True)
adata
sc.pp.calculate_qc_metrics(adata, inplace=True)
sc.pl.spatial(adata, img_key='hires', color="n_genes_by_counts")
AnnData object with n_obs × n_vars = 2797 × 32285
obs: 'in_tissue', 'array_row', 'array_col'
var: 'gene_ids', 'feature_types', 'genome'
uns: 'spatial'
obsm: 'spatial'
在这里还遇到一个具有借鉴意义的bug。我从10x官网的公开数据中下载了数据用来测试，该数据是由SpaceRanger v2.0处理的。在使用sc.read_visium读入时，会提示“OSError: Could not find '\path\spatial\tissue_positions_list.csv'”。原因是从SpaceRanger v2.0开始，此文件已重命名为“tissue_positions.csv”并包含标题行。如果你的数据也是由SpaceRanger v2.0产生的，那么注意使用scanpy 1.9.4以上版本进行分析。
Stereo-seq
Stereo-seq是由华大基因开发的具有亚细胞分辨率的空间转录组技术，测序reads通过SAW流程分析产生表达矩阵，由同样由华大开发的Python包stereopy进行读取。Stereo-seq的基因表达矩阵格式如下：
使用stereopy包中的stereo.io.read_gem()读取矩阵为StereoExpData格式。
例：
import stereo as st
import warnings
warnings.filterwarnings('ignore')
data_path = '/path/Mouse_brain.gem'
data = st.io.read_gem(
file_path=data_path,
bin_type='bins',
bin_size=50,
is_sparse=True,
)
data.tl.cal_qc()
data.plt.spatial_scatter()
StereoExpData object with n_cells X n_genes = 38741 X 26177
bin_type: bins
bin_size: 50
offset_x = 3225
offset_y = 6175
cells: ['cell_name']
genes: ['gene_name']
通过stereopy读取后，可将StereoExpData转换为anndata从而使用scanpy进行分析，而anndata的h5ad文件又可以进一步转化为rds（https://github.com/STOmics/stereopy/blob/dev/docs/source/_static/annh5ad2rds.R）从而便于使用Seurat进行分析。
#Scanpy
adata = st.io.stereo_to_anndata(data,flavor='scanpy',output='scanpy_out.h5ad')
#Seurat
data.tl.cal_qc()
adata = st.io.stereo_to_anndata(data,flavor='seurat',output='seurat_out.h5ad')
Slide-seq
Slide-seq由Broad研究所开发，数据经由Slide-seq tools流程处理。从目前已公开的数据来看，Slide-seq的结果文件类型多样，各个公开数据库中下载到的Slide-seq数据类型都有差别，但其主要数据由1个基因表达矩阵和1个beads矩阵构成，两个文件可能为分别的两个文本文件或整合为loom格式文件等。
Seurat官方给出的读入方法是：首先使用CreateSeuratObject将表达矩阵读入，然后将坐标信息附加至SeuratObject中。
例：
library(Seurat)
library(ggplot2)
slideCount <- '/path/Slide-seq/MappedDGEForR.csv'
slideBead <- '/path/Slide-seq/BeadLocationsForR.csv'
bead <- ReadSlideSeq(slideBead, assay = 'Spatial')
count <- read.table(slideCount, header = TRUE, sep=',', row.names='Row')
obj <- CreateSeuratObject(count, assay ="Spatial")
bead <- bead[Cells(x = obj)]
obj[["slice1"]] <- bead
plot <- SpatialFeaturePlot(obj, features = "nCount_Spatial") + theme(legend.position = "right")
plot
至此，三种空间转录组技术所产生的数据就可以被读入R或Python中进行后续的数据分析，后续我们也将继续从0开始，带大家一起学习单细胞/空间转录组的分析方法。
参考资料：
https://satijalab.org/seurat/articles/spatial_vignette
https://scanpy-tutorials.readthedocs.io/en/latest/spatial/basic-analysis.html
https://stereopy.readthedocs.io/en/latest/Tutorials/IO.html
https://patrickcnmartin.github.io/Vesalius/articles/Vesalius_Analysis/Vesalius_analysis.html
https://squidpy.readthedocs.io/en/latest/notebooks/tutorials/tutorial_read_spatial.html
