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

相关文章推荐
跑龙套的圣诞树 · LocalDate - 抱糖果彡 - 博客园· 1 月前 ·
欢乐的篮球 · Array.prototype.splice ...· 5 月前 ·
任性的沙发 · Android项目中引入aar包的正确方法介 ...· 6 月前 ·
烦恼的紫菜汤 · java 校验字符串是否是 json ...· 1 年前 ·
谈吐大方的小蝌蚪 · 转发 ...· 1 年前 ·
Reading data
QC and preprocessing
Manifold embedding and clustering based on transcriptional similarity
Visualization in spatial coordinates
簇的marker基因
空间数据的高变基因
MERFISH示例
本篇内容演示如何在Scanpy中使用空间转录组数据（spatial transcriptomics data），首先，我们专注于10x Genomics Visium data，最后，我们为MERFISH技术的空间转录数据分析提供了示例。
目前空间转录组技术基本上还未到单细胞水平，一般一个空间位置spots上有多个细胞，10X Genomics Visum一般为1-10个细胞，空间转录组学能够提供空间位置，但是分辨率达不到单细胞水平，单细胞转录组学分辨率能够达到单细胞分辨率，但是没有空间位置信息，因此将空间转录组与单细胞转录组结合的话，那就既有空间位置信息，又达到了单细胞分辨率，这为科学界提供了有力的工具，大大的提高了生物体的组织、器官等研究的分辨率和准确性。
我们导入scanpy工具：
import scanpy as sc
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
sc.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 × 36601
    obs: '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 × 36601
    obs: '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 × 19685
    obs: '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 cell
    finished (0:00:00)
If you pass `n_top_genes`, all cutoffs are ignored.
extracting highly variable genes
    finished (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 × 19685
    obs: '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 PCA
    on highly variable genes
    with n_comps=50
    finished (0:00:01)
computing neighbors
    using 'X_pca' with n_pcs = 50
    finished: added to `.uns['neighbors']`
    `.obsp['distances']`, distances for each pair of neighbors
    `.obsp['connectivities']`, weighted adjacency matrix (0:00:06)
computing UMAP
    finished: added
    'X_umap', UMAP coordinates (adata.obsm) (0:00:06)
running Leiden clustering
    finished: found 10 clusters and added
    'clusters', the cluster labels (adata.obs, categorical) (0:00:00)
adata
AnnData object with n_obs × n_vars = 3861 × 19685
    obs: '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 genes
    finished: 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.request
url_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 × 12903
    obsm: '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 cell
    finished (0:00:00): normalized adata.X and added    'n_counts', counts per cell before normalization (adata.obs)
computing PCA
    with n_comps=15
    finished (0:00:00)
computing neighbors
    using 'X_pca' with n_pcs = 15
    finished: added to `.uns['neighbors']`
    `.obsp['distances']`, distances for each pair of neighbors
    `.obsp['connectivities']`, weighted adjacency matrix (0:00:05)
computing UMAP
    finished: added
    'X_umap', UMAP coordinates (adata.obsm) (0:00:01)
running Leiden clustering
    finished: 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 × 12903
    obs: '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中使用空间转录组数据（spatial transcriptomics data），首先，我们专注于10x Genomics Visium data，最后，我们为MERFISH技术的空间转录数据分析提供了示例。目前空间转录组技术基本上还未到单细胞水平，一般一个空间位置spots上有多个细胞，10X Genomics Visum一般为1-10个细胞，空间转录组学能够提供空间位置，但是分辨率达不到单细胞水平，单细胞转录组学分辨率能够达到单细胞分辨率，但是没有空间位置信息，因此将空间转
				CytoSPACE：scRNA-seq数据到空间转录组学数据的最佳映射
CytoSPACE是一种新颖的计算策略，用于在空间转录组（ST）测量可能包含多个细胞的贡献的情况下，将单细胞转录组分配给原位空间转录组数据。 我们的方法通过基于线性编程的优化例程将基于相关的成本函数最小化，从而解决了单个像元/点分配问题。
 该存储库包含用于实现和评估我们的方法的代码以及一个应用该方法的案例研究。
 我们方法的关键创新是：
 与常规方法相比，CytoSPACE在单个细胞水平上解剖给定组织中细胞的空间组织。
 由于我们的方法从scRNA测序数据中绘制了单个细胞，与可用的空间转录组学技术相比，每个细胞中都有大量的基因被测序，因此我们的方法显着改善了重建组织的基因覆盖率。
 我们的方法不需要有关细胞类型和细胞状态的先验知识。
 主要实现是作为Python 3软件包。 要查看SpatialDE的用法示例，请继续




    

				使用Seurat对空间数据集进行分析，可视化和集成
本教程演示了如何使用Seurat（> = 3.2）分析空间分布的RNA-seq数据。尽管分析流程类似于单细胞RNA序列分析的Seurat工作流程，但我们引入了更新的交互作用和可视化工具，特别着重于空间和分子信息的整合。本教程将涵盖以下任务，我们认为这些任务在许多空间分析中都是常见的：
Normalization标准化
Dimensional reduction and clustering降维和聚类
Detecting spatially
我把常用的函数写成了几个包，方便之后使用；
bioquest 包括三个子包 tl、pl、st 分别是常用的工具包 DdatafFame 的处理、画图、字符串处理；
genekit 包括基因名转换、格式转换、差异分析、提取 TCGA 数据等的函数；
sckit 包括单细胞分析的一些函数；
可以在https://jihulab.com/BioQuest找到这些函数。
读入 spaceranger 的输出文件
在读入和后续分析的时候会出现报错，解决方法如下
				之前讲过一篇空间转录组的文献，里面首次提出了Multimodal intersection analysis(MIA)的空间转录组分析思路。
讲解视频在B站
MIA分析可以用来评估空间上某个region或者cluster中富集的细胞类型。需要单细胞和空间转录组两种组学数据，数据最好配对。
MIA原理
上图是示例图，一个region是否富含某一种细胞类型，看的是一个region高表达的基因和一个celltype高表达的基因是不是有足够多的重叠。
原理不算难，跟常规的「差异基因」做「富集分析」原理是一样的，详细
				01、空间转录组技术的发展
近年来单细胞转录组测序技术的应用大大拓宽了人们的视野，使人们能够深入了解组织中细胞的构成的多样性和基因表达状态。众所周知，基因表达具有时间和空间的特异性，通过对不同时间点的样本取材，使用单细胞转录组测序技术能够解析时间维度上细胞类型和基因表达的变化过程。
图1. 早期胚胎发育中基因表达的时间特异性【1】
然而单细胞测序实验的前提是组织必须通过机械分离或酶解消化成单细胞悬液，此过程不可避免的丢失了组织中细胞所处的原始位置信息，也导致了细胞间的通讯网络被打破，这使我们难以获得组织中不
				目录PBMCsMapping PBMCs using ingest
下面的tutorial描述了一种基于PCA的简单方法用于集成数据（integrating data），我们称之为ingest，并将其与BBKNN进行了比较。此外，BBKNN可以与Scanpy工作流良好结合，可通过bbknn接口访问。
ingest函数假设有一个带注释的参考数据集，该数据集捕获了感兴趣的生物变异性。rational是在参考数据上拟合模型，并使用它来投影新数据。目前，该模型是一个结合了邻居查找搜索树（neighbor looku
				本篇内容介绍如何使用多个Visium数据集，以及如何用scRNA-seq数据集进行整合。我们将使用Scanorama来进行整合。首先我们要安装Scanorama：
加载相关框架：
建议先阅读论文阅读笔记-利用Scanorama高效整合异质单细胞转录组我们将使用两个小鼠大脑的Visium空间转录组数据集，该数据集可从10x genomics website获取。函数从10x genomics下载数据集并返回adata对象（包含counts，images和spatial coordinates），我们使用计算
				文章目录一、安装二、使用1、准备工作2、预处理过滤低质量细胞样本3、检测特异性基因4、主成分分析（Principal component analysis）5、领域图，聚类图（Neighborhood graph）6、检索标记基因7、保存数据8、番外
如果没有conda 基础，参考： Conda 安装使用图文详解（2021版）
pip install scanpy
conda install -y -c conda-forge leidenalg
1、准备工作
# 载入包
import