目前空间转录组技术基本上还未到单细胞水平，一般一个空间位置spots上有多个细胞，10X Genomics Visum一般为1-10个细胞，空间转录组学能够提供空间位置，但是分辨率达不到单细胞水平，单细胞转录组学分辨率能够达到单细胞分辨率，但是没有空间位置信息，因此将空间转录组与单细胞转录组结合的话，那就既有空间位置信息，又达到了单细胞分辨率，这为科学界提供了有力的工具，大大的提高了生物体的组织、器官等研究的分辨率和准确性。

我们导入scanpy工具：

import scanpy as sc
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as snssc.logging.print_versions()
sc.set_figure_params(facecolor="white", figsize=(8, 8))
sc.settings.verbosity = 3"""
anndata     0.8.0
scanpy      1.9.1
"""

Reading data

我们将使用人类淋巴结（human lymphnode）的Visium空间转录组学数据集，该数据集发布在10x genomics website：link；

函数 datasets.visium_sge() 从10x Genomics下载数据集，并返回一个包含计数矩阵counts和空间坐标 spatital coordinates 的AnnData对象。

adata = sc.datasets.visium_sge(sample_id="V1_Human_Lymph_Node")
adata.var_names_make_unique()adata
"""
AnnData object with n_obs × n_vars = 4035 × 36601obs: 'in_tissue', 'array_row', 'array_col'var: 'gene_ids', 'feature_types', 'genome'uns: 'spatial'obsm: 'spatial'
"""

我们将使用 pp.calculate_qc_metrics 并基于线粒体 read counts 计算QC指标。

# 将线粒体基因组保存为注释 var.mt
adata.var["mt"] = adata.var_names.str.startswith("MT-")
adata.var["mt"]
"""
MIR1302-2HG    False...
AC007325.2     False
Name: mt, Length: 36601, dtype: bool
"""# 计算指标, qc的var选择 var.mt
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], inplace=True)adata
"""
AnnData object with n_obs × n_vars = 4035 × 36601obs: 'in_tissue', 'array_row', 'array_col', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes', 'pct_counts_in_top_200_genes', 'pct_counts_in_top_500_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt'var: 'gene_ids', 'feature_types', 'genome', 'mt', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts'uns: 'spatial'obsm: 'spatial'
"""

回顾Scanpy(三)处理3k PBMCs并聚类中的内容，我们可以了解到某些质量的含义：

n_genes_by_counts：每个细胞中，有表达的基因的个数；
total_counts：每个细胞的基因总计数（总表达量）；
pct_counts_mt：每个细胞中，线粒体基因表达量占该细胞所有基因表达量的百分比；

QC and preprocessing

我们根据 total_counts 和 n_genes_by_counts 进行一些基本过滤。我们先将这些质量指标可视化为直方图：

fig, axs = plt.subplots(1, 4, figsize=(15, 4))
sns.distplot(adata.obs["total_counts"], kde=False, ax=axs[0])
sns.distplot(adata.obs["total_counts"][adata.obs["total_counts"] < 10000], kde=False, bins=40, ax=axs[1])
sns.distplot(adata.obs["n_genes_by_counts"], kde=False, bins=60, ax=axs[2])
sns.distplot(adata.obs["n_genes_by_counts"][adata.obs["n_genes_by_counts"] < 4000], kde=False, bins=60, ax=axs[3])

第二个图是第一个图在0到10000区间的可视化，total_counts 代表一个细胞的基因总表达量。第四个图是第三个图在0到4000的可视化，n_genes_by_counts 代表一个细胞中，有表达的基因的个数。注意，计数矩阵为4035×36601。

通过上面的可视化结果，我们保留 total_counts 在5000到35000的细胞，下一步，我们保留线粒体基因 pct_counts_mt 占比小于20%的细胞。最后，做基因上的过滤：

sc.pp.filter_cells(adata, min_counts=5000)
sc.pp.filter_cells(adata, max_counts=35000)
adata = adata[adata.obs["pct_counts_mt"] < 20]
print(f"#cells after MT filter:{adata.n_obs}")
sc.pp.filter_genes(adata, min_cells=10)"""
filtered out 44 cells that have less than 5000 counts
filtered out 130 cells that have more than 35000 counts
#cells after MT filter: 3861
filtered out 16916 genes that are detected in less than 10 cells
"""adata
"""
AnnData object with n_obs × n_vars = 3861 × 19685obs: 'in_tissue', 'array_row', 'array_col', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes', 'pct_counts_in_top_200_genes', 'pct_counts_in_top_500_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'n_counts'var: 'gene_ids', 'feature_types', 'genome', 'mt', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'n_cells'uns: 'spatial'obsm: 'spatial'
"""

我们继续使用scanpy的内置函数normaliza_total标准化Visium counts data（包含counts和空间转录信息的数据），并在此之后检测高变基因。

sc.pp.normalize_total(adata, inplace=True)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, flavor="seurat", n_top_genes=2000)"""
normalizing counts per cellfinished (0:00:00)
If you pass `n_top_genes`, all cutoffs are ignored.
extracting highly variable genesfinished (0:00:00)
--> added'highly_variable', boolean vector (adata.var)'means', float vector (adata.var)'dispersions', float vector (adata.var)'dispersions_norm', float vector (adata.var)
"""adata
"""
AnnData object with n_obs × n_vars = 3861 × 19685obs: 'in_tissue', 'array_row', 'array_col', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes', 'pct_counts_in_top_200_genes', 'pct_counts_in_top_500_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'n_counts'var: 'gene_ids', 'feature_types', 'genome', 'mt', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'n_cells', 'highly_variable', 'means', 'dispersions', 'dispersions_norm'uns: 'spatial', 'log1p', 'hvg'obsm: 'spatial'
"""

Manifold embedding and clustering based on transcriptional similarity

首先，我们按照标准聚类步骤进行（pca降维，计算neighbors graph，umap降维，leiden聚类）：

sc.pp.pca(adata)
sc.pp.neighbors(adata)
sc.tl.umap(adata)
sc.tl.leiden(adata, key_added="clusters")"""
computing PCAon highly variable geneswith n_comps=50finished (0:00:01)
computing neighborsusing 'X_pca' with n_pcs = 50finished: added to `.uns['neighbors']``.obsp['distances']`, distances for each pair of neighbors`.obsp['connectivities']`, weighted adjacency matrix (0:00:06)
computing UMAPfinished: added'X_umap', UMAP coordinates (adata.obsm) (0:00:06)
running Leiden clusteringfinished: found 10 clusters and added'clusters', the cluster labels (adata.obs, categorical) (0:00:00)
"""adata
"""
AnnData object with n_obs × n_vars = 3861 × 19685obs: 'in_tissue', 'array_row', 'array_col', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes', 'pct_counts_in_top_200_genes', 'pct_counts_in_top_500_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'n_counts', 'clusters'var: 'gene_ids', 'feature_types', 'genome', 'mt', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'n_cells', 'highly_variable', 'means', 'dispersions', 'dispersions_norm'uns: 'spatial', 'log1p', 'hvg', 'pca', 'neighbors', 'umap', 'leiden'obsm: 'spatial', 'X_pca', 'X_umap'varm: 'PCs'obsp: 'distances', 'connectivities'
"""

我们基于umap可视化一些相关变量，以检查是否存在与总计数（total_counts）和检测到的基因（n_genes_by_counts）相关的特定数据分布：

plt.rcParams["figure.figsize"] = (4, 4)
sc.pl.umap(adata, color=["total_counts", "n_genes_by_counts", "clusters"], wspace=0.4)

Visualization in spatial coordinates

我们可以先看看空间坐标：

adata.obsm['spatial'].shape
# (3861, 2)adata.obsm['spatial']
"""
array([[8346, 6982],[4270, 1363],...,[4822, 1840]], dtype=int64)
"""

现在让我们来看看 total_counts 和 n_genes_by_counts 在空间坐标中的分布。我们将使用 sc.pl.spatial 函数将圆形点（circular spots）覆盖在提供的 Hematoxylin and eosin stain（H＆E）图像上。

plt.rcParams["figure.figsize"] = (8, 8)
sc.pl.spatial(adata, img_key="hires", color=["total_counts", "n_genes_by_counts"])

函数sc.pl.spatial接受4个附加参数：

img_key：存储在adata.uns中images里的key；
crop_coord：用于裁剪的坐标（left, right, top, bottom）；
alpha_img：透明度；
bw：将图像转换为灰度图的标志；

此外，在sc.pl.spatial中，size参数会改变其行为：它会成为spots大小的比例因子。

之前，我们在基因表达空间中进行聚类，并使用UMAP将结果可视化。现在我们通过在空间维度上可视化样本，我们可以深入了解组织的结构，并可能深入了解细胞间的通信：

sc.pl.spatial(adata, img_key="hires", color="clusters", size=1.5)

在基因表达空间中属于同一簇的spots通常在空间维度上同时出现。例如，属于簇4的点通常被属于簇0的点包围。

我们可以放大特定的感兴趣区域，以获得特定的结论。此外，通过改变spot的alpha值，我们可以更好地从H&E图像中看到潜在的组织形态。

sc.pl.spatial(adata, img_key="hires", color="clusters", groups=["0", "4"], alpha=0.5, size=1.3)

簇的marker基因

我们进一步检查簇4。我们计算标记基因并绘制一个热图，图中显示了簇中前10个标记基因的表达水平。

sc.tl.rank_genes_groups(adata, "clusters", method="t-test")
sc.pl.rank_genes_groups_heatmap(adata, groups="4", n_genes=10, groupby="clusters")"""
ranking genesfinished: added to `.uns['rank_genes_groups']`'names', sorted np.recarray to be indexed by group ids'scores', sorted np.recarray to be indexed by group ids'logfoldchanges', sorted np.recarray to be indexed by group ids'pvals', sorted np.recarray to be indexed by group ids'pvals_adj', sorted np.recarray to be indexed by group ids (0:00:02)
WARNING: dendrogram data not found (using key=dendrogram_clusters). Running `sc.tl.dendrogram` with default parameters. For fine tuning it is recommended to run `sc.tl.dendrogram` independently.using 'X_pca' with n_pcs = 50
Storing dendrogram info using `.uns['dendrogram_clusters']`
WARNING: Groups are not reordered because the `groupby` categories and the `var_group_labels` are different.
categories: 0, 1, 2, etc.
var_group_labels: 4
"""

我们对基因CR2进行分析：

sc.pl.spatial(adata, img_key="hires", color=["clusters", "CR2"])

空间数据的高变基因

空间转录组学使研究人员能够研究基因表达趋势在空间中的变化，从而确定基因表达的空间模式。为此，我们使用SpatialDE，这是一种基于高斯过程的统计框架，旨在识别空间转录组的可变基因：

pip install spatialde

最近，还提出了其他用于识别空间变异基因的工具，例如：SPARK，trendsceek，HMRF

首先，我们将标准化的counts和coordinates转换为pandas dataframe，这是输入spatialDE所需的：

import SpatialDE%%time
counts = pd.DataFrame(adata.X.todense(), columns=adata.var_names, index=adata.obs_names)
coord = pd.DataFrame(adata.obsm['spatial'], columns=['x_coord', 'y_coord'], index=adata.obs_names)"""
CPU times: total: 93.8 ms
Wall time: 94.7 ms
"""

查看counts：

查看coord：

运行SpatialDE需要相当长的时间（本示例接近2小时）。

results = SpatialDE.run(coord, counts)

results中会保存基于空间转录组数据计算得到的可变基因。

MERFISH示例

如果我们使用基于FISH的技术生成空间数据，只需读cordinate表并将其分配给adata.obsm即可。

让我们看关于Xia等人2019年的例子。

首先，我们需要下载坐标和counts数据：

import urllib.requesturl_coord = "https://www.pnas.org/highwire/filestream/887973/field_highwire_adjunct_files/15/pnas.1912459116.sd15.xlsx"
filename_coord = "pnas.1912459116.sd15.xlsx"
urllib.request.urlretrieve(url_coord, filename_coord)url_counts = "https://www.pnas.org/highwire/filestream/887973/field_highwire_adjunct_files/12/pnas.1912459116.sd12.csv"
filename_counts = "pnas.1912459116.sd12.csv"
urllib.request.urlretrieve(url_counts, filename_counts)

如果不能访问以上链接，可以手动到github下载：spatial datasets，这是一个空间转录数据集的整理。

然后读取到adata：

coordinates = pd.read_excel("./pnas.1912459116.sd15.xlsx", index_col=0)
counts = sc.read_csv("./pnas.1912459116.sd12.csv").transpose()adata_merfish = counts[coordinates.index, :]
adata_merfish.obsm["spatial"] = coordinates.to_numpy()adata_merfish
"""
AnnData object with n_obs × n_vars = 645 × 12903obsm: 'spatial'
"""

我们将执行标准的预处理和降维：

sc.pp.normalize_per_cell(adata_merfish, counts_per_cell_after=1e6)
sc.pp.log1p(adata_merfish)
sc.pp.pca(adata_merfish, n_comps=15)
sc.pp.neighbors(adata_merfish)
sc.tl.umap(adata_merfish)
sc.tl.leiden(adata_merfish, key_added="clusters", resolution=0.5)"""
normalizing by total count per cellfinished (0:00:00): normalized adata.X and added    'n_counts', counts per cell before normalization (adata.obs)
computing PCAwith n_comps=15finished (0:00:00)
computing neighborsusing 'X_pca' with n_pcs = 15finished: added to `.uns['neighbors']``.obsp['distances']`, distances for each pair of neighbors`.obsp['connectivities']`, weighted adjacency matrix (0:00:05)
computing UMAPfinished: added'X_umap', UMAP coordinates (adata.obsm) (0:00:01)
running Leiden clusteringfinished: found 7 clusters and added'clusters', the cluster labels (adata.obs, categorical) (0:00:00)
"""adata_merfish
"""
AnnData object with n_obs × n_vars = 645 × 12903obs: 'n_counts', 'clusters'uns: 'log1p', 'pca', 'neighbors', 'umap', 'leiden'obsm: 'spatial', 'X_pca', 'X_umap'varm: 'PCs'obsp: 'distances', 'connectivities'
"""

我们可以在UMAP空间和spatial坐标中可视化Leiden得到的簇。

sc.pl.umap(adata_merfish, color="clusters")
sc.pl.embedding(adata_merfish, basis="spatial", color="clusters")

Scanpy(六)空间转录组数据的分析与可视化相关推荐

Python对阿里巴巴、谷歌、腾讯等六家公司股票数据进行分析与可视化实战（附源码超详细）
需要源码请点赞关注收藏后评论区留言私信~~~ 下面针对阿里巴巴.谷歌.亚马逊.Facebook.苹果和腾讯六家公司股票数据进行了分析与可视化描述,数据分析前需要安装互联数据获取包pandas-data ...
热点综述 | 单细胞+空间转录组的整合分析方法总结
目前scRNA-seq将每个转录物与单个细胞相关联,但关于这些转录物在组织中的位置信息丢失了:相反的,空间转录组学技术知道转录物的位置,却不知道是哪个细胞产生了转录物.因此,scRNA-seq与空间转 ...
单细胞转录组数据整合分析专题研讨会(2019.11)
2019年10月9日,单细胞转录组再等Nature.题为Decoding human fetal liver haematopoiesis的研究,对受孕后4周至17周的人胚胎肝脏.卵黄囊.肾脏和皮肤组 ...
爬取最好大学网数据、分析并可视化操作
爬虫爬取数据.分析并可视化操作本次对最好大学网进行爬虫示例. 1.获取网页响应 def getHTMLText(url):try:resp = request.urlopen(url)html_da ...
使用python生成词云——聆心云心理健康服务平台数据可视分析和可视化
实验题目:聆心云心理健康服务平台数据可视分析和可视化实验目的和要求:统计出在聆心云平台做沙盘游戏的次数.根据各次沙盘游戏所使用的沙具和进行的操作数据进行词云可视化,掌握Python词云制作方法实验 ...
【空间转录组】MIA分析
之前讲过一篇空间转录组的文献,里面首次提出了Multimodal intersection analysis(MIA)的空间转录组分析思路. 讲解视频在B站 MIA分析可以用来评估空间上某个regio ...
文本数据可视化_非结构化文本数据的分析和可视化
文本数据可视化 Stuck behind the paywall? Read this article with my friend link here. 卡在收费墙后面? 在这里与我的朋友链接阅读本 ...
对新闻数据的分析与可视化
该任务是在大学本科的时候的一个小作业,主要是对湾湾的新闻数据进行一个采集分析与可视化,任务比较简单.任务可分为两个部分,一个是爬取中时新闻网站的数据,二是对该数据进行处理与分析. 一.爬取数据数据: ...
波士顿犯罪数据时空分析及可视化
文章目录前言一.数据描述二.数据预处理三.犯罪时空分析 1.犯罪类型分析 2.犯罪时间特征分析 (1)基于年维度 (2)基于月维度 (3)基于日维度 3.犯罪空间特征分析四.基于犯罪空间理论 ...

Scanpy(六)空间转录组数据的分析与可视化

目录