scanpy官方教程2022||04-数据整合：ingest a

作者: 信你个鬼 | 来源:发表于2022-08-28 18:24 被阅读0次

学习资料来源：

scanpy主页：https://scanpy.readthedocs.io/en/stable/
官网：https://scanpy-tutorials.readthedocs.io/en/latest/integrating-data-using-ingest.html【注意教程有两个版本，这里是latest版本的学习笔记】

此教程描述了一种简单的基于pca的方法来整合数据，我们称之为ingest，并将其与BBKNN [Polanski19]进行了比较。BBKNN与Scanpy工作流整合得很好，可以通过BBKNN功能访问。

ingest函数假设有一个带注释的参考数据集，该数据集捕获了感兴趣的生物变异。合理的做法是在参考数据上拟合一个模型，并使用它来投射新的数据。

目前，这个模型是一个结合了邻居查找搜索树的PCA，我们将使用UMAP来实现这个功能。类似的基于pca的整合之前已经被使用过，例如，[Weinreb18]。ingest的优点：

ingest简单，流程清晰且运行快
与BBKNN一样，ingest不改变数据matrix本身
不同于BBKNN的是，ingest解决了标签映射问题(如scmap)，并保持了可能具有特定cluster或轨迹等所需属性的嵌入

我们将这种非对称数据集整合称为将注释从带注释的参考adata_ref摄取到没有注释的数据中。它不同于学习以对称方式整合数据集的joint，如BBKNN, Scanorma, Conos, CCA(如Seurat)或有条件的VAE(如scVI, trVAE)，但可与scran中的初始MNN实现相比较。

01 PBMCs数据

用到两个数据集：一个带注释的参考数据集adata_ref和一个需要被注释和嵌入数据的数据集。

import scanpy as sc
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as pl
# verbosity: errors (0), warnings (1), info (2), hints (3)
sc.settings.verbosity = 1  
sc.logging.print_versions()
sc.settings.set_figure_params(dpi=80, frameon=False, figsize=(3, 3), facecolor='white')
outdir = '/Pub/Users/zhangjuan/project/scanpy/Intergrating/'

# this is an earlier version of the dataset from the pbmc3k tutorial
adata_ref = sc.datasets.pbmc3k_processed()  
adata = sc.datasets.pbmc68k_reduced()

环境中软件版本如下：

-----
anndata     0.8.0
scanpy      1.9.1
-----
PIL                 9.2.0
beta_ufunc          NA
binom_ufunc         NA
colorama            0.4.5
cycler              0.10.0
cython_runtime      NA
dateutil            2.8.2
h5py                3.7.0
hypergeom_ufunc     NA
igraph              0.9.11
joblib              1.1.0
kiwisolver          1.4.4
leidenalg           0.8.10
llvmlite            0.38.1
louvain             0.7.1
matplotlib          3.4.3
mpl_toolkits        NA
natsort             8.1.0
nbinom_ufunc        NA
numba               0.55.2
numexpr             2.8.3
numpy               1.22.4
packaging           21.3
pandas              1.4.3
pkg_resources       NA
pyparsing           3.0.9
pytz                2022.1
scipy               1.8.1
seaborn             0.11.2
session_info        1.0.0
setuptools          63.2.0
six                 1.16.0
sklearn             1.1.1
statsmodels         0.13.2
texttable           1.6.4
threadpoolctl       3.1.0
-----
Python 3.9.13 | packaged by conda-forge | (main, May 27 2022, 16:58:50) [GCC 10.3.0]
Linux-5.10.0-9-amd64-x86_64-with-glibc2.31
-----
Session information updated at 2022-08-03 23:55

参考数据以及查询数据的情况：

adata_ref
AnnData object with n_obs × n_vars = 2638 × 1838
    obs: 'n_genes', 'percent_mito', 'n_counts', 'louvain'
    var: 'n_cells'
    uns: 'draw_graph', 'louvain', 'louvain_colors', 'neighbors', 'pca', 'rank_genes_groups'
    obsm: 'X_pca', 'X_tsne', 'X_umap', 'X_draw_graph_fr'
    varm: 'PCs'
    obsp: 'distances', 'connectivities'

# 注释：
adata_ref.obs
                  n_genes  percent_mito  n_counts          louvain
index                                                             
AAACATACAACCAC-1      781      0.030178    2419.0      CD4 T cells
AAACATTGAGCTAC-1     1352      0.037936    4903.0          B cells
AAACATTGATCAGC-1     1131      0.008897    3147.0      CD4 T cells
AAACCGTGCTTCCG-1      960      0.017431    2639.0  CD14+ Monocytes
AAACCGTGTATGCG-1      522      0.012245     980.0         NK cells
...                   ...           ...       ...              ...
TTTCGAACTCTCAT-1     1155      0.021104    3459.0  CD14+ Monocytes
TTTCTACTGAGGCA-1     1227      0.009294    3443.0          B cells
TTTCTACTTCCTCG-1      622      0.021971    1684.0          B cells
TTTGCATGAGAGGC-1      454      0.020548    1022.0          B cells
TTTGCATGCCTCAC-1      724      0.008065    1984.0      CD4 T cells

[2638 rows x 4 columns]


adata
AnnData object with n_obs × n_vars = 700 × 765
    obs: 'bulk_labels', 'n_genes', 'percent_mito', 'n_counts', 'S_score', 'G2M_score', 'phase', 'louvain'
    var: 'n_counts', 'means', 'dispersions', 'dispersions_norm', 'highly_variable'
    uns: 'bulk_labels_colors', 'louvain', 'louvain_colors', 'neighbors', 'pca', 'rank_genes_groups'
    obsm: 'X_pca', 'X_umap'
    varm: 'PCs'
    obsp: 'distances', 'connectivities'

使用sc.tl函数，需要在相同的变量（即需要提共有基因数据）上定义数据集

var_names = adata_ref.var_names.intersection(adata.var_names)
adata_ref = adata_ref[:, var_names]
adata = adata[:, var_names]

只有208个共同基因！

在参考数据上训练的模型和图形(这里是PCA，邻居，UMAP)将解释其中观察到的生物变异。

sc.pp.pca(adata_ref)
sc.pp.neighbors(adata_ref)
sc.tl.umap(adata_ref)

该流形看起来与聚类教程中的流形基本相同。

sc.pl.umap(adata_ref, color='louvain')
pl.savefig(outdir + "01-adata_ref-UMAP.png")

结果图：

image-20220726093757581.png

使用ingest整合pbmc

基于选择的特征，我们使用adata_ref数据集合的labels and embeddings映射到adata。在这里，我们使用adata_ref.obsm['X_pca']来map 细胞亚群标签 and the UMAP 坐标。

sc.tl.ingest函数：

obs：Labels’ keys in adata_ref.obs which need to be mapped to adata.obs (inferred for observation of adata)

sc.tl.ingest(adata, adata_ref, obs='louvain')
# fix colors
adata.uns['louvain_colors'] = adata_ref.uns['louvain_colors'] 
sc.pl.umap(adata, color=['louvain', 'bulk_labels'], wspace=0.5)
## 修改图片边距，之前图片保存字体显示不全
pl.margins(0, 0)
pl.savefig(outdir + "02-adata-UMAP.png", bbox_inches='tight', dpi=300, pad_inches=0.1)

结果图：这个图片跟官方的结果是个镜像，180度反的，左边louvain注释是通过ingest整合的结果，右边bulk_labels是数据原有的注释标签：

image-20220727102219103.png

通过比较louvain注释结果与bulk_labels注释结果，我们发现基本上都能对上，除了DC细胞。但是DC细胞可能在adata中就存在注释问题。

adata_concat = adata_ref.concatenate(adata, batch_categories=['ref', 'new'])
adata_concat.obs.louvain = adata_concat.obs.louvain.astype('category')
# fix category ordering
adata_concat.obs.louvain.cat.reorder_categories(adata_ref.obs.louvain.cat.categories, inplace=True)  
# fix category colors
adata_concat.uns['louvain_colors'] = adata_ref.uns['louvain_colors']  
sc.pl.umap(adata_concat, color=['batch', 'louvain'])
pl.margins(0, 0)
pl.savefig(outdir + "03-adata_ref-UMAP-new.png", bbox_inches='tight', dpi=300, pad_inches=0.1)

矫正之后的结果：虽然在单核细胞和树突状细胞簇中似乎存在一些批量效应，但是新的数据在其他方面被映射得相对均匀。

巨核细胞megakaryoctes只存在于 adata _ ref 中，adata中没有cell映射到它们

如果交换adata _ ref 和adata，巨核细胞megakaryoctes不再作为一个单独的cluster出现。这是一个极端的情况，因为参考数据非常小; 但是人们应该始终质疑参考数据是否包含足够的生物变异来有意义地容纳查询数据。

image-20220801084713388.png

使用BBKNN整合数据

adata_concat是adata_ref（2638个细胞）与adata（700个细胞）两个矩阵合并后的数据

参数join默认采用inner方式即交集合并数据 : str (default: 'inner')，Use intersection ('inner') or union ('outer') of variables

合并后，adata_concat包含：3338个细胞和208个共同基因

adata_concat
AnnData object with n_obs × n_vars = 3338 × 208
    obs: 'n_genes', 'percent_mito', 'n_counts', 'louvain', 'bulk_labels', 'S_score', 'G2M_score', 'phase', 'batch'
    var: 'n_counts-new', 'means-new', 'dispersions-new', 'dispersions_norm-new', 'highly_variable-new', 'n_cells-ref'
    uns: 'louvain_colors', 'batch_colors', 'pca'
    obsm: 'X_pca', 'X_umap'
    varm: 'PCs'

整合数据：

sc.tl.pca(adata_concat)
# running bbknn 1.3.6
sc.external.pp.bbknn(adata_concat, batch_key='batch')  
sc.tl.umap(adata_concat)
sc.pl.umap(adata_concat, color=['batch', 'louvain'])
pl.savefig(outdir + "04-adata_concat-UMAP-BBKNN.png",bbox_inches='tight', dpi=300, pad_inches=0.1)

结果图：左图为adata_ref数据与adata数据整合在一起的情况，右边为注释后的结果

image-20220801091835116.png

02 Pancreas

本数据集在the scGen paper [Lotfollahi19]，here中使用过，在here被核实过，可以在here (the BBKNN paper)这里进行下载。

它包含了来自4个不同研究（Segerstolpe16, Baron16, Wang16, Muraro16）的人类胰腺的数据，它已经在单细胞数据集整合的开创性论文中被使用(Butler18，Haghverdi18) ，并且从那时起被多次使用。

将数据下载下来并存放在data目录下：

# note that this collection of batches is already intersected on the genes
adata_all = sc.read('data/pancreas.h5ad', backup_url='https://www.dropbox.com/s/qj1jlm9w10wmt0u/pancreas.h5ad?dl=1')
adata_all.shape

(14693, 2448)

检查在这些研究中观察到的细胞类型：

counts = adata_all.obs.celltype.value_counts()
counts

alpha                     4214
beta                      3354
ductal                    1804
acinar                    1368
not applicable            1154
delta                      917
gamma                      571
endothelial                289
activated_stellate         284
dropped                    178
quiescent_stellate         173
mesenchymal                 80
macrophage                  55
PSC                         54
unclassified endocrine      41
co-expression               39
mast                        32
epsilon                     28
mesenchyme                  27
schwann                     13
t_cell                       7
MHC class II                 5
unclear                      4
unclassified                 2
Name: celltype, dtype: int64

为了简化可视化，让我们删除5个细胞数量比较少的cluster：

# get the minority classes
minority_classes = counts.index[-5:].tolist() 
# actually subset
adata_all = adata_all[~adata_all.obs.celltype.isin(minority_classes)]
# reorder according to abundance
adata_all.obs.celltype.cat.reorder_categories(counts.index[:-5].tolist(), inplace=True)

查看批次效应：

sc.pp.pca(adata_all)
sc.pp.neighbors(adata_all)
sc.tl.umap(adata_all)
sc.pl.umap(adata_all, color=['batch', 'celltype'], palette=sc.pl.palettes.vega_20_scanpy)
pl.savefig(outdir + "05-adata_all-UMAP.png",bbox_inches='tight', dpi=300, pad_inches=0.1)

可以看到很明显的批次效应：这里跟官网的图有点不一样，聚类形状差不多，但是就是有些cluster旋转了某个角度

image-20220804002237472.png

使用BBKNN整合数据

上面数据中的批次可以很好地使用 BBKNN [Polanski19]来解决。

sc.external.pp.bbknn(adata_all, batch_key='batch')
sc.tl.umap(adata_all)
sc.pl.umap(adata_all, color=['batch', 'celltype'])
pl.savefig(outdir + "06-adata_all-UMAP_batchAdjust.png",bbox_inches='tight', dpi=300, pad_inches=0.1)

结果如下：

image-20220805001043795.png

使用ingest整合数据

选择一个参考批处理来训练模型和建立邻域图(这里是PCA) ，并分离出所有其他批处理。

和以前一样，对参考批次进行训练的模型将解释在其中观察到的生物学变异。

这里数据有4个batch，选择batch 0 作为adata_ref，其余为adata

adata_ref = adata_all[adata_all.obs.batch == '0']

# Compute the PCA, neighbors and UMAP on the reference data
sc.pp.pca(adata_ref)
sc.pp.neighbors(adata_ref)
sc.tl.umap(adata_ref)

# The reference batch contains 12 of the 19 cell types across all batches.
sc.pl.umap(adata_ref, color='celltype')
pl.savefig(outdir + "07-adata_all-UMAP_ingest.png",bbox_inches='tight', dpi=300, pad_inches=0.1)

结果图如下：参考batch中包含所有19中细胞类型中的12种细胞

image-20220805082551546.png

迭代地将label(例如‘ 细胞类型’)和embeddings(例如‘ X _ pca’和‘ X _ umap’)从参考数据映射到需要注释的数据中：

adatas = [adata_all[adata_all.obs.batch == i].copy() for i in ['1', '2', '3']]
# a bit more logging
sc.settings.verbosity = 2  
for iadata, adata in enumerate(adatas):
    print(f'... integrating batch {iadata+1}')
    # save the original cell type
    adata.obs['celltype_orig'] = adata.obs.celltype  
    sc.tl.ingest(adata, adata_ref, obs='celltype')
    

... integrating batch 1
running ingest
    finished (0:00:15)
... integrating batch 2
running ingest
    finished (0:00:06)
... integrating batch 3
running ingest
    finished (0:00:02)

现在，每个查询的batch都带有已经与 adata _ ref 进行上下文关联的注释

使用concatenate，我们可以一起可视化看看结果：

adata_concat = adata_ref.concatenate(adatas)
adata_concat.obs.celltype = adata_concat.obs.celltype.astype('category')
# fix category ordering
adata_concat.obs.celltype.cat.reorder_categories(adata_ref.obs.celltype.cat.categories, inplace=True)  
# fix category coloring 
adata_concat.uns['celltype_colors'] = adata_ref.uns['celltype_colors']  
sc.pl.umap(adata_concat, color=['batch', 'celltype'])
pl.savefig(outdir + "08-adata_concat-UMAP_ingest.png",bbox_inches='tight', dpi=300, pad_inches=0.1)

结果图：与BBKNN的结果相比，ingest 以一种更明显的方式维持cluster。如果已经观察到想要的连续结构(例如在造血数据集中)，ingest 可以轻松地维持这个结构。

image-20220813005417132-16611794386863.png

评估一致性

提取数据中刚刚被当做查询的部分，即batch1,2,3，包含6113个细胞

## 评估一致性
adata_query = adata_concat[adata_concat.obs.batch.isin(['1', '2', '3'])]
adata_query
View of AnnData object with n_obs × n_vars = 6113 × 2448

# The following plot is a bit hard to read, hence, move on to confusion matrices below
sc.pl.umap(adata_query, color=['batch', 'celltype', 'celltype_orig'], wspace=0.4)
pl.savefig(outdir + "09-adata_query-consistency.png",bbox_inches='tight', dpi=300, pad_inches=0.1)

结果图：整合后数据的batch，整合后数据的注释，数据原有的注释label

image-20220813010127651.png

不同batches中保守的细胞类型

我们首先来查看在reference batch与query batch中保守的细胞类型：

obs_query = adata_query.obs

# intersected categories
conserved_categories = obs_query.celltype.cat.categories.intersection(obs_query.celltype_orig.cat.categories)  

# intersect categories
obs_query_conserved = obs_query.loc[obs_query.celltype.isin(conserved_categories) & obs_query.celltype_orig.isin(conserved_categories)]  

# remove unused categoriyes
obs_query_conserved.celltype.cat.remove_unused_categories(inplace=True)  

# remove unused categoriyes
obs_query_conserved.celltype_orig.cat.remove_unused_categories(inplace=True)  

# fix category ordering
obs_query_conserved.celltype_orig.cat.reorder_categories(obs_query_conserved.celltype.cat.categories, inplace=True)  

pd.crosstab(obs_query_conserved.celltype, obs_query_conserved.celltype_orig)

总的来说，保守的细胞类型按预期那样映射。意外的是原注释中一些腺泡细胞表现为腺泡细胞。然而，已经观察到的参考数据显示腺泡和导管细胞聚集，这解释了差异，并表明初始注释中可能存在不一致。

image-20220823000738415.png

所有的细胞类型

现在，我们来看看所有的细胞类型：

pd.crosstab(adata_query.obs.celltype, adata_query.obs.celltype_orig)

结果：

我们观察到PSC(胰腺星状细胞)细胞实际上只是不一致地注释和正确地映射在“激活星状细胞”上；
此外，很高兴看到“间充质”和“间充质”细胞都属于同一类别。然而，这个类别仍然是'activated_stellate'，而且很可能是错误的

image-20220823000658941.png

各批次的可视化分布

通常，批量与人们想要比较的实验相对应。Scanpy为此提供了方便的可视化可能性

Density plot

sc.tl.embedding_density(adata_concat, groupby='batch')
sc.pl.embedding_density(adata_concat, groupby='batch')
pl.savefig(outdir + "10-adata_concat-density_plot.png",bbox_inches='tight', dpi=300, pad_inches=0.1)

结果图：

image-20220823001148726.png

嵌入群子集的部分可视化

for batch in ['1', '2', '3']:
    sc.pl.umap(adata_concat, color='batch', groups=[batch])
    out = outdir + "11-adata_concat-Partial_visualizaton"+batch+".png"
    pl.savefig(out, bbox_inches='tight', dpi=300, pad_inches=0.1)

image-20220823002643460.png

下回见~

scanpy官方教程2022||04-数据整合：ingest a

学习资料来源：

01 PBMCs数据

使用ingest整合pbmc

使用BBKNN整合数据

02 Pancreas

使用BBKNN整合数据

使用ingest整合数据

评估一致性

不同batches中保守的细胞类型

所有的细胞类型

各批次的可视化分布

Density plot

嵌入群子集的部分可视化

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

单细胞测序

scanpy