【生信实战】从马赛克到4K，极度震撼的Visium HD空间转录组

近日10X终于发布了他们拥有亚细胞分辨率的空间转录组产品，并给出了一些列的测试数据。今天带大家进行一些基本的探索与实战

1 数据下载

这里笔者选择了10X给到的小鼠的脑组织进行测试，具体的下载路径可以见
https://www.10xgenomics.com/datasets/visium-hd-cytassist-gene-expression-libraries-of-mouse-brain-he

image.png

2 数据分析

空转转录组的分析主要使用R、python。由于Visium HD的超高的spot数量，在8bin的情况下整个表达矩阵包含了近400W个基因，因此我们后续的分析默认使用scanpy完成。

import os
import scanpy as sc
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import anndata

#使用jutypter notebook加上这两句
%matplotlib inline
%config InlineBackend.figure_format = 'retina'

#首先导入数据，我们将数据下载到/data/outs/binned_outputs中
#这里出现了一个坑，由于scanpy并没有适配Visium HD，会导致无法读取数据，这是由于坐标文件由原来的csv转为了parquet，这里首先进行转换
tissue_position_file = '/data/outs/binned_outputs/square_008um/spatial/tissue_positions.parquet'
tissue_position_csv = '/data/outs/binned_outputs/square_008um/spatial/tissue_positions_list.csv'

if  not os.path.exists(tissue_position_csv):
    tissue_position_df = pd.read_parquet(tissue_position_file)
    tissue_position_df.to_csv(tissue_position_csv, index=False, header=None)

adata = sc.read_visium(path = '/data/outs/binned_outputs/square_008um')

数据基本概况

整个数据无比庞大，我们接着往后进行细胞质控和过滤等操作

adata.var_names_make_unique() 
# qc
adata.var["mt"] = adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(
    adata, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True
)
sc.pl.violin(
    adata,
    ["n_genes_by_counts", "total_counts", "pct_counts_mt"],
    jitter=0.4,
    multi_panel=True
)
sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt')
sc.pl.scatter(adata, x="total_counts", y="n_genes_by_counts")
sc.pl.scatter(adata, "total_counts", "n_genes_by_counts", color="pct_counts_mt")

数据质控

线粒体占比

#这里当然也可以选择不进行过滤，由于数据过大，导致spot的数据相比我们常规的单细胞明显的不足，进而单个spot的转录本较少，所以比较适合较为宽松的阈值。过滤太狠的情况下出现较多的空spot
sc.pp.filter_cells(adata, min_genes=50)
sc.pp.filter_cells(adata, max_genes=8000)
adata = adata[(adata.obs.pct_counts_mt < 20)]
sc.pp.filter_genes(adata, min_cells=3)

随后就是比较常规的HVG鉴定、归一化和降维聚类，这一步需要一定的时间，请耐心等待

sc.pp.highly_variable_genes(adata, n_top_genes= 2000, flavor='seurat_v3')
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.pca(adata, n_comps=50, use_highly_variable=True)
sc.pp.neighbors(adata, n_neighbors=20, n_pcs=30)
sc.tl.umap(adata, min_dist = 0.3)
sc.tl.leiden(adata, key_added="seurat_clusters", resolution=0.5)
sc.tl.tsne(adata, n_pcs = 30)

然后我们看一下结果，保存到clustering.png文件中

fig, ax = plt.subplots(2,2,figsize=(12, 10), dpi=300)
sc.pl.tsne(adata, color = "seurat_clusters", wspace=0.4, title='clustering using tsne', ax = ax[0,0], show=False)
sc.pl.umap(adata, color = "seurat_clusters", wspace=0.4, title = 'clustering using umap', ax = ax[0,1], show=False)
point_size = 1
sc.pl.spatial(adata, img_key="hires", color="n_genes_by_counts", size=point_size, ax = ax[1,0], show=False)
sc.pl.spatial(adata, img_key="hires", color="seurat_clusters", size=point_size, ax = ax[1,1], show=False)
fig.tight_layout()
fig.savefig('clustering.png')
fig.savefig('clustering.pdf')

输出结果

小鼠大脑分区示意图

非HD的版本

整个分辨率的提升真是令人叹为观止

3 题外

如果使用Seurat进行分析时，记得使用最新的版本，安装命令如下。如果使用先前的版本也会出现不兼容的问题

# packages required for Visium HD
install.packages("hdf5r")
install.packages("arrow")
remotes::install_github(repo = "satijalab/seurat", ref = "visium-hd")