scanpy教程：使用ingest和BBKNN整合多样本
男，
一个长大了才会遇到的帅哥，
稳健，潇洒，大方，靠谱。
一段生信缘，一棵技能树，
一枚大型测序工厂的螺丝钉，
一个随机森林中提灯觅食的津门旅客。
回顾
正文
随着单细胞技术的成熟，测序成本的降低，单细胞的数据量和样本量也日益增长。我们知道单细胞转录组的一个主要应用就是解释细胞的异质性，那么，不同器官，不同测序平台，不同物种之间的单细胞数据何如整合分析呢？特别是在单细胞的数据维度这么高的前提下，显然传统的基于回归的方法已经不适用了。于是出现了一批单细胞整合分析的工具，它们大多数是在R生态条件下的。如：
在我们理解单细胞数据的时候一张cell X gene 的大表不能离开我们的脑海。
1adata.to_df()
2Out[24]:
3 RP11-34P13.3 FAM138A ... AC213203.1 FAM231B
4AAACCCAAGCGTATGG-1 0.0 0.0 ... 0.0 0.0
5AAACCCAGTCCTACAA-1 0.0 0.0 ... 0.0 0.0
6AAACCCATCACCTCAC-1 0.0 0.0 ... 0.0 0.0
7AAACGCTAGGGCATGT-1 0.0 0.0 ... 0.0 0.0
8AAACGCTGTAGGTACG-1 0.0 0.0 ... 0.0 0.0
9 ... ... ... ... ...
10TTTGTTGCAGGTACGA-1 0.0 0.0 ... 0.0 0.0
11TTTGTTGCAGTCTCTC-1 0.0 0.0 ... 0.0 0.0
12TTTGTTGGTAATTAGG-1 0.0 0.0 ... 0.0 0.0
13TTTGTTGTCCTTGGAA-1 0.0 0.0 ... 0.0 0.0
14TTTGTTGTCGCACGAC-1 0.0 0.0 ... 0.0 0.0
当我们有多个样本的时候就是有多张这样的表，那让我们自己手动来整合这两张表的话，我们会怎么做呢？
肯定是行列分别对齐把它们拼在一起啊，就像拼积木一样的，但是这样的结果就是：
两个样本在图谱上完全的分开来了。我们不同平台的样本，相同的细胞类型应该是在一起的啊。于是我们开始思考如何完成这样的整合。
seurat提供了一套解决方案，就是在数据集中构建锚点，将不同数据集中相似的细胞锚在一起。
那么如何锚，选择哪些特征来锚定，又开发出不同的算法。不管算法如何，首先我们看看这种锚定可以为我们带来什么？相同的细胞类型mapping在一起，一个自然的作用就是用来mapping细胞类型未知的数据。
所以在scanpy中也如seurat一样在多样本分析中，分别给出reference的方法和整合的方法。目前在scanpy中分别是ingest和BBKNN（Batch balanced kNN），当然整合也是可以用来做reference的。scanpy.external.pp.mnn_correct(https://scanpy.readthedocs.io/en/latest/external/scanpy.external.pp.mnn_correct.html#scanpy-external-pp-mnn-correct)应该也是可以用的。
先来看ingest，通过投射到参考数据上的PCA(或备用模型)上，将一个adata的嵌入和注释与一个参考数据集adata_ref集成在一起。该函数使用knn分类器来映射标签，使用UMAP来映射嵌入。
再来看看bbknn(https://github.com/Teichlab/bbknn)是一个快速和直观的批处理效果去除工具，可以直接在scanpy工作流中使用。它是scanpy.api.pp.neighbors()的替代方法，这两个函数都创建了一个邻居图，以便后续在集群、伪时间和UMAP可视化中使用。标准方法首先确定整个数据结构中每个单元的k个最近邻，然后将候选单元转换为指数相关的连接，然后作为进一步分析的基础。
那么我们就来看一下在scanpy的实现吧。
1import scanpy as sc
2import pandas as pd
3import seaborn as sns
4import sklearn
5import sys
6import scipy
7import bbknn1sc.settings.verbosity = 1 # verbosity: errors (0), warnings (1), info (2), hints (3)
2sc.logging.print_versions()
3sc.settings.set_figure_params(dpi=80, frameon=False, figsize=(3, 3))1scanpy==1.4.5.1 anndata==0.7.1 umap==0.3.10 numpy==1.16.5 scipy==1.3.1 pandas==0.25.1 scikit-learn==0.21.3 statsmodels==0.10.1 python-igraph==0.8.0
ingest 注释
1adata_ref = sc.datasets.pbmc3k_processed() # this is an earlier version of the dataset from the pbmc3k tutorial
2
3adata_ref
4AnnData object with n_obs × n_vars = 2638 × 1838
5 obs: 'n_genes', 'percent_mito', 'n_counts', 'louvain'
6 var: 'n_cells'
7 uns: 'draw_graph', 'louvain', 'louvain_colors', 'neighbors', 'pca', 'rank_genes_groups'
8 obsm: 'X_pca', 'X_tsne', 'X_umap', 'X_draw_graph_fr'
9 varm: 'PCs'
我们一次看看以下参考数据集都有哪些内容：
1adata_ref.obs
2Out[9]:
3 n_genes percent_mito n_counts louvain
4index
5AAACATACAACCAC-1 781 0.030178 2419.0 CD4 T cells
6AAACATTGAGCTAC-1 1352 0.037936 4903.0 B cells
7AAACATTGATCAGC-1 1131 0.008897 3147.0 CD4 T cells
8AAACCGTGCTTCCG-1 960 0.017431 2639.0 CD14+ Monocytes
9AAACCGTGTATGCG-1 522 0.012245 980.0 NK cells
10 ... ... ... ...
11TTTCGAACTCTCAT-1 1155 0.021104 3459.0 CD14+ Monocytes
12TTTCTACTGAGGCA-1 1227 0.009294 3443.0 B cells
13TTTCTACTTCCTCG-1 622 0.021971 1684.0 B cells
14TTTGCATGAGAGGC-1 454 0.020548 1022.0 B cells
15TTTGCATGCCTCAC-1 724 0.008065 1984.0 CD4 T cells
16
17[2638 rows x 4 columns]1adata_ref.var
2Out[10]:
3 n_cells
4index
5TNFRSF4 155
6CPSF3L 202
7ATAD3C 9
8C1orf86 501
9RER1 608
10 ...
11ICOSLG 34
12SUMO3 570
13SLC19A1 31
14S100B 94
15PRMT2 588
16
17[1838 rows x 1 columns]1adata_ref.uns['louvain_colors']
2Out[14]:
3array(['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b',
4 '#e377c2', '#bcbd22'], dtype='<U7')
5
6 adata_ref.obsm
7Out[16]: AxisArrays with keys: X_pca, X_tsne, X_umap, X_draw_graph_fr
8
9adata_ref.obsm['X_umap']
10Out[17]:
11array([[ 1.35285574, 2.26612719],
12 [-0.47802448, 7.87730423],
13 [ 2.16588875, -0.24481226],
14 ...,
15 [ 0.34670979, 8.34967798],
16 [ 0.19864146, 9.56698797],
17 [ 2.62803322, 0.36722543]])
有没有再次理解AnnData 这个对象的数据结构呢？
可以看到在这个数据集中降维聚类都是做过的，我们可以画个图看看：
1sc.pl.umap(adata_ref, color='louvain')
接下来我们看看要预测的数据集是怎样的。
1adata = sc.datasets.pbmc68k_reduced()
2adata
3
4AnnData object with n_obs × n_vars = 700 × 765
5 obs: 'bulk_labels', 'n_genes', 'percent_mito', 'n_counts', 'S_score', 'G2M_score', 'phase', 'louvain'
6 var: 'n_counts', 'means', 'dispersions', 'dispersions_norm', 'highly_variable'
7 uns: 'bulk_labels_colors', 'louvain', 'louvain_colors', 'neighbors', 'pca', 'rank_genes_groups'
8 obsm: 'X_pca', 'X_umap'
9 varm: 'PCs'
可见它也是降维聚类过的了。
1sc.pl.umap(adata, color='louvain')
这个数据集并没有得到细胞类型的定义。
构建注释数据结构：
1var_names = adata_ref.var_names.intersection(adata.var_names) # 取交集
2adata_ref = adata_ref[:, var_names]
3adata = adata[:, var_names]1sc.pp.pca(adata_ref)
2sc.pp.neighbors(adata_ref)
3sc.tl.umap(adata_ref)
4sc.tl.leiden(adata_ref)# 新的聚类方法
5sc.pl.umap(adata_ref, color=['louvain','leiden'])
6adata_ref
7
8AnnData object with n_obs × n_vars = 2638 × 208
9 obs: 'n_genes', 'percent_mito', 'n_counts', 'louvain', 'leiden'
10 var: 'n_cells'
11 uns: 'draw_graph', 'louvain', 'louvain_colors', 'neighbors', 'pca', 'rank_genes_groups', 'umap', 'leiden', 'leiden_colors'
12 obsm: 'X_pca', 'X_tsne', 'X_umap', 'X_draw_graph_fr'
13 varm: 'PCs'
1sc.pp.pca(adata)
2sc.pp.neighbors(adata)
3sc.tl.umap(adata)
4sc.tl.leiden(adata)
5sc.pl.umap(adata, color=['louvain','leiden'])
6adata
7
8AnnData object with n_obs × n_vars = 700 × 208
9 obs: 'bulk_labels', 'n_genes', 'percent_mito', 'n_counts', 'S_score', 'G2M_score', 'phase', 'louvain', 'leiden'
10 var: 'n_counts', 'means', 'dispersions', 'dispersions_norm', 'highly_variable'
11 uns: 'bulk_labels_colors', 'louvain', 'louvain_colors', 'neighbors', 'pca', 'rank_genes_groups', 'umap', 'leiden', 'leiden_colors'
12 obsm: 'X_pca', 'X_umap'
13 varm: 'PCs'
用ingest来做细胞注释吧。
1sc.tl.ingest(adata, adata_ref, obs='louvain')
2adata.uns['louvain_colors'] = adata_ref.uns['louvain_colors'] # fix colors
我们来看看sc.tl.ingest的帮助文档：
1Help on function ingest in module scanpy.tools._ingest:
2
3ingest(adata: anndata._core.anndata.AnnData, adata_ref: anndata._core.anndata.AnnData, obs: Union[str, Iterable[str], NoneType] = None, embedding_method: Union[str, Iterable[str]] = ('umap', 'pca'), labeling_method: str = 'knn', inplace: bool = True, **kwargs)
4 Map labels and embeddings from reference data to new data.
5
6 :tutorial:`integrating-data-using-ingest`
7
8 Integrates embeddings and annotations of an `adata` with a reference dataset
9 `adata_ref` through projecting on a PCA (or alternate
10 model) that has been fitted on the reference data. The function uses a knn
11 classifier for mapping labels and the UMAP package [McInnes18]_ for mapping
12 the embeddings.
13
14 .. note::
15
16 We refer to this *asymmetric* dataset integration as *ingesting*
17 annotations from reference data to new data. This is different from
18 learning a joint representation that integrates both datasets in an
19 unbiased way, as CCA (e.g. in Seurat) or a conditional VAE (e.g. in
20 scVI) would do.
21
22 You need to run :func:`~scanpy.pp.neighbors` on `adata_ref` before
23 passing it.
24
25 Parameters
26 ----------
27 adata
28 The annotated data matrix of shape `n_obs` × `n_vars`. Rows correspond
29 to cells and columns to genes. This is the dataset without labels and
30 embeddings.
31 adata_ref
32 The annotated data matrix of shape `n_obs` × `n_vars`. Rows correspond
33 to cells and columns to genes.
34 Variables (`n_vars` and `var_names`) of `adata_ref` should be the same
35 as in `adata`.
36 This is the dataset with labels and embeddings
37 which need to be mapped to `adata`.
38 obs
39 Labels' keys in `adata_ref.obs` which need to be mapped to `adata.obs`
40 (inferred for observation of `adata`).
41 embedding_method
42 Embeddings in `adata_ref` which need to be mapped to `adata`.
43 The only supported values are 'umap' and 'pca'.
44 labeling_method
45 The method to map labels in `adata_ref.obs` to `adata.obs`.
46 The only supported value is 'knn'.
47 inplace
48 Only works if `return_joint=False`.
49 Add labels and embeddings to the passed `adata` (if `True`)
50 or return a copy of `adata` with mapped embeddings and labels.
51
52 Returns
53 -------
54 * if `inplace=False` returns a copy of `adata`
55 with mapped embeddings and labels in `obsm` and `obs` correspondingly
56 * if `inplace=True` returns `None` and updates `adata.obsm` and `adata.obs`
57 with mapped embeddings and labels
58
59 Example
60 -------
61 Call sequence:
62
63import scanpy as sc
64sc.pp.neighbors(adata_ref)
65sc.tl.umap(adata_ref)
66sc.tl.ingest(adata, adata_ref, obs='cell_type')
67
68 .. _ingest PBMC tutorial: https://scanpy-tutorials.readthedocs.io/en/latest/integrating-pbmcs-using-ingest.html
69 .. _ingest Pancreas tutorial: https://scanpy-tutorials.readthedocs.io/en/latest/integrating-pancreas-using-ingest.html
70
通过比较'bulk_label’注释和'louvain’注释，我们发现数据被合理地映射，只有树突细胞的注释似乎是含糊不清的，在adata中可能已经是模糊的了。我们来对adata做进一步的处理。
1adata_concat = adata_ref.concatenate(adata, batch_categories=['ref', 'new'])
2adata_concat.obs.louvain = adata_concat.obs.louvain.astype('category')
3adata_concat.obs.louvain.cat.reorder_categories(adata_ref.obs.louvain.cat.categories, inplace=True) # fix category ordering
4adata_concat.uns['louvain_colors'] = adata_ref.uns['louvain_colors'] # fix category colors
5adata_concat
6sc.pl.umap(adata_concat, color=['batch', 'louvain'])
7
8AnnData object with n_obs × n_vars = 3338 × 208
9 obs: 'G2M_score', 'S_score', 'batch', 'bulk_labels', 'leiden', 'louvain', 'n_counts', 'n_genes', 'percent_mito', 'phase'
10 var: 'n_cells-ref', 'n_counts-new', 'means-new', 'dispersions-new', 'dispersions_norm-new', 'highly_variable-new'
11 obsm: 'X_pca', 'X_umap'
虽然在单核细胞和树突状细胞簇中似乎存在一些批处理效应，但在其他方面，新数据被绘制得相对均匀。
巨核细胞只存在于adata_ref中，没有来自adata映射的单元格。如果交换参考数据和查询数据，巨核细胞不再作为单独的集群出现。这是一个极端的情况，因为参考数据非常小;但是，人们应该始终质疑参考数据是否包含足够的生物变异，以便有意义地容纳查询数据。
使用BBKNN整合
1sc.tl.pca(adata_concat)
2sc.external.pp.bbknn(adata_concat, batch_key='batch') # running bbknn 1.3.6
3sc.tl.umap(adata_concat)
4sc.pl.umap(adata_concat, color=['batch', 'louvain'])
1adata_concat
2Out[45]:
3AnnData object with n_obs × n_vars = 3338 × 208
4 obs: 'G2M_score', 'S_score', 'batch', 'bulk_labels', 'leiden', 'louvain', 'n_counts', 'n_genes', 'percent_mito', 'phase'
5 var: 'n_cells-ref', 'n_counts-new', 'means-new', 'dispersions-new', 'dispersions_norm-new', 'highly_variable-new'
6 uns: 'batch_colors', 'louvain_colors', 'pca', 'neighbors', 'umap'
7 obsm: 'X_pca', 'X_umap'
8 varm: 'PCs'
BBKNN并不维持巨核细胞簇。然而，它似乎更均匀地混合细胞。
一个使用BBKNN整合数据的例子
以下数据已在scGen论文[Lotfollahi19]中使用。点击pancreas(http://ftp//ngs.sanger.ac.uk/production/teichmann/BBKNN/objects-pancreas.zip)下载数据。
它包含了来自4个不同研究(Segerstolpe16, Baron16, Wang16, Muraro16)的人类胰腺数据，这些数据在单细胞数据集集成的开创性论文(Butler18, Haghverdi18)中被使用过，并在此后多次被使用。
1h5ad = 'E:\\learnscanpy\\data\\objects-pancreas\\pancreas.h5ad'
2adata_all = sc.read_h5ad(h5ad)
3adata_all
4
5AnnData object with n_obs × n_vars = 14693 × 2448
6 obs: 'celltype', 'sample', 'n_genes', 'batch', 'n_counts', 'louvain'
7 var: 'n_cells-0', 'n_cells-1', 'n_cells-2', 'n_cells-3'
8 uns: 'celltype_colors', 'louvain', 'neighbors', 'pca', 'sample_colors'
9 obsm: 'X_pca', 'X_umap'
10 varm: 'PCs'1 counts = adata_all.obs.celltype.value_counts()
2counts
3Out[173]:
4alpha 4214
5beta 3354
6ductal 1804
7acinar 1368
8not applicable 1154
9delta 917
10gamma 571
11endothelial 289
12activated_stellate 284
13dropped 178
14quiescent_stellate 173
15mesenchymal 80
16macrophage 55
17PSC 54
18unclassified endocrine 41
19co-expression 39
20mast 32
21epsilon 28
22mesenchyme 27
23schwann 13
24t_cell 7
25MHC class II 5
26unclear 4
27unclassified 2
28Name: celltype, dtype: int641adata_all.obs
2Out[171]:
3 celltype sample ... n_counts louvain
4index ...
5human1_lib1.final_cell_0001-0 acinar Baron ... 2.241100e+04 2
6human1_lib1.final_cell_0002-0 acinar Baron ... 2.794900e+04 2
7human1_lib1.final_cell_0003-0 acinar Baron ... 1.689200e+04 2
8human1_lib1.final_cell_0004-0 acinar Baron ... 1.929900e+04 2
9human1_lib1.final_cell_0005-0 acinar Baron ... 1.506700e+04 2
10 ... ... ... ... ...
11reads.29499-3 ductal Wang ... 1.056558e+06 10
12reads.29500-3 ductal Wang ... 9.926309e+05 10
13reads.29501-3 beta Wang ... 1.751338e+06 10
14reads.29502-3 dropped Wang ... 2.163764e+06 10
15reads.29503-3 beta Wang ... 2.038979e+06 10
16
17[14693 rows x 6 columns]
可以看出这个数据集已经降维聚类好了，所以我们可以可视化一下：
1sc.pl.umap(adata_all,color=['sample', 'celltype','louvain'])
样本之间的批次很严重啊。
去掉细胞数较小的小群，
1minority_classes = counts.index[-5:].tolist() # get the minority classes
2
3# ['schwann', 't_cell', 'MHC class II', 'unclear', 'unclassified']
4
5adata_all = adata_all[ # actually subset
6 ~adata_all.obs.celltype.isin(minority_classes)]
7adata_all.obs.celltype.cat.reorder_categories( # reorder according to abundance
8 counts.index[:-5].tolist(), inplace=True)
9
10adata_all.obs.celltype.value_counts()
11Out[182]:
12alpha 4214
13beta 3354
14ductal 1804
15acinar 1368
16not applicable 1154
17delta 917
18gamma 571
19endothelial 289
20activated_stellate 284
21dropped 178
22quiescent_stellate 173
23mesenchymal 80
24macrophage 55
25PSC 54
26unclassified endocrine 41
27co-expression 39
28mast 32
29epsilon 28
30mesenchyme 27
进行pca降维和umap降维：
1sc.pp.pca(adata_all)
2sc.pp.neighbors(adata_all)
3sc.tl.umap(adata_all)
4sc.pl.umap(adata_all, color=['batch', 'celltype'], palette=sc.pl.palettes.vega_20_scanpy)
下面我们使用BBKNN来整合数据：
1sc.external.pp.bbknn(adata_all, batch_key='batch')
2sc.tl.umap(adata_all)
3adata_all
4sc.pl.umap(adata_all, color=['sample','batch', 'celltype'])
果然要比原始的数据好多了。但是改变的是什么？
1AnnData object with n_obs × n_vars = 14662 × 2448
2 obs: 'celltype', 'sample', 'n_genes', 'batch', 'n_counts', 'louvain'
3 var: 'n_cells-0', 'n_cells-1', 'n_cells-2', 'n_cells-3'
4 uns: 'celltype_colors', 'louvain', 'neighbors', 'pca', 'sample_colors', 'louvain_colors', 'umap', 'batch_colors'
5 obsm: 'X_pca', 'X_umap'
6 varm: 'PCs'
如果想对其中某个样本进行单独的注释，可以用上面提到的ingest。选择一个参考批次来训练模型和建立邻域图(这里是一个PCA)，并分离出所有其他批次。
1adata_ref = adata_all[adata_all.obs.batch == '0']
2adata_ref
3Out[191]:
4View of AnnData object with n_obs × n_vars = 8549 × 2448
5 obs: 'celltype', 'sample', 'n_genes', 'batch', 'n_counts', 'louvain'
6 var: 'n_cells-0', 'n_cells-1', 'n_cells-2', 'n_cells-3'
7 uns: 'celltype_colors', 'louvain', 'neighbors', 'pca', 'sample_colors', 'louvain_colors', 'umap', 'batch_colors'
8 obsm: 'X_pca', 'X_umap'
9 varm: 'PCs'1sc.pp.pca(adata_ref)
2sc.pp.neighbors(adata_ref)
3sc.tl.umap(adata_ref)
4
5adata_ref
6Out[197]:
7AnnData object with n_obs × n_vars = 8549 × 2448
8 obs: 'celltype', 'sample', 'n_genes', 'batch', 'n_counts', 'louvain'
9 var: 'n_cells-0', 'n_cells-1', 'n_cells-2', 'n_cells-3'
10 uns: 'celltype_colors', 'louvain', 'neighbors', 'pca', 'sample_colors', 'louvain_colors', 'umap', 'batch_colors'
11 obsm: 'X_pca', 'X_umap'
12 varm: 'PCs'
13
14sc.pl.umap(adata_ref, color='celltype')
选取数据集用ingest于adata_ref进行mapping：
1adatas = [adata_all[adata_all.obs.batch == i].copy() for i in ['1', '2', '3']]
2sc.settings.verbosity = 2 # a bit more logging
3for iadata, adata in enumerate(adatas):
4 print(f'... integrating batch {iadata+1}')
5 adata.obs['celltype_orig'] = adata.obs.celltype # save the original cell type
6 sc.tl.ingest(adata, adata_ref, obs='celltype')
7
8integrating batch 1
9running ingest
10 finished (0:00:08)
11integrating batch 2
12running ingest
13 finished (0:00:06)
14integrating batch 3
15running ingest
16 finished (0:00:03)
17
18adata_concat = adata_ref.concatenate(adatas)
19adata_concat
20
21Out[200]:
22AnnData object with n_obs × n_vars = 14662 × 2448
23 obs: 'batch', 'celltype', 'celltype_orig', 'louvain', 'n_counts', 'n_genes', 'sample'
24 var: 'n_cells-0', 'n_cells-1', 'n_cells-2', 'n_cells-3'
25 obsm: 'X_pca', 'X_umap'1adata_concat.obs.celltype = adata_concat.obs.celltype.astype('category')
2adata_concat.obs.celltype.cat.reorder_categories(adata_ref.obs.celltype.cat.categories, inplace=True) # fix category ordering
3adata_concat.uns['celltype_colors'] = adata_ref.uns['celltype_colors'] # fix category coloring
4
5sc.pl.umap(adata_concat, color=['celltype_orig','batch', 'celltype'])
与BBKNN的结果相比，这是以一种更加明显的方式保持分群。如果已经观察到一个想要的连续结构(例如在造血数据集中)，摄取允许容易地维持这个结构。
一致性评估
1adata_query = adata_concat[adata_concat.obs.batch.isin(['1', '2', '3'])]
2
3View of AnnData object with n_obs × n_vars = 6113 × 2448
4 obs: 'batch', 'celltype', 'celltype_orig', 'louvain', 'n_counts', 'n_genes', 'sample'
5 var: 'n_cells-0', 'n_cells-1', 'n_cells-2', 'n_cells-3'
6 uns: 'celltype_colors', 'celltype_orig_colors', 'batch_colors'
7 obsm: 'X_pca', 'X_umap'1sc.pl.umap(
2 adata_query, color=['batch', 'celltype', 'celltype_orig'], wspace=0.4)
这个结果依然不能很好的反映一致性，让我们首先关注与参考保守的细胞类型，以简化混淆矩阵的reads。
1obs_query = adata_query.obs
2conserved_categories = obs_query.celltype.cat.categories.intersection(obs_query.celltype_orig.cat.categories) # intersected categories
3obs_query_conserved = obs_query.loc[obs_query.celltype.isin(conserved_categories) & obs_query.celltype_orig.isin(conserved_categories)] # intersect categories
4obs_query_conserved.celltype.cat.remove_unused_categories(inplace=True) # remove unused categoriyes
5obs_query_conserved.celltype_orig.cat.remove_unused_categories(inplace=True) # remove unused categoriyes
6obs_query_conserved.celltype_orig.cat.reorder_categories(obs_query_conserved.celltype.cat.categories, inplace=True) # fix category ordering
7
8obs_query_conserved
9Out[214]:
10 batch celltype celltype_orig ... n_counts n_genes sample
11D28.1_1-1-1 1 alpha alpha ... 2.322583e+04 5448 Muraro
12D28.1_13-1-1 1 ductal ductal ... 2.334263e+04 5911 Muraro
13D28.1_15-1-1 1 alpha alpha ... 2.713471e+04 5918 Muraro
14D28.1_17-1-1 1 alpha alpha ... 1.581207e+04 4522 Muraro
15D28.1_2-1-1 1 endothelial endothelial ... 3.173151e+04 6464 Muraro
16 ... ... ... ... ... ... ...
17reads.29498-3-3 3 ductal ductal ... 1.362606e+06 19950 Wang
18reads.29499-3-3 3 ductal ductal ... 1.056558e+06 19950 Wang
19reads.29500-3-3 3 ductal ductal ... 9.926309e+05 19950 Wang
20reads.29501-3-3 3 beta beta ... 1.751338e+06 19950 Wang
21reads.29503-3-3 3 beta beta ... 2.038979e+06 19950 Wang1pd.crosstab(obs_query_conserved.celltype, obs_query_conserved.celltype_orig)
2Out[215]:
3celltype_orig alpha beta ductal acinar delta gamma endothelial mast
4celltype
5alpha 1819 3 7 0 1 20 0 5
6beta 49 803 3 1 10 26 0 0
7ductal 8 5 693 263 0 0 0 0
8acinar 1 3 2 145 0 3 0 0
9delta 5 4 0 0 305 73 0 0
10gamma 1 5 0 0 0 194 0 0
11endothelial 2 0 0 0 0 0 36 0
12mast 0 0 1 0 0 0 0 2
总的来说，保守的细胞类型也如预期的那样被映射。主要的例外是原始注释中出现的一些腺泡细胞。然而，已经观察到参考数据同时具有腺泡和导管细胞，这就解释了差异，并指出了初始注释中潜在的不一致性。
现在让我们继续看看所有的细胞类型。
1 pd.crosstab(adata_query.obs.celltype, adata_query.obs.celltype_orig)
2Out[216]:
3celltype_orig PSC acinar ... not applicable unclassified endocrine
4celltype ...
5alpha 0 0 ... 304 11
6beta 0 1 ... 522 24
7ductal 0 263 ... 106 1
8acinar 0 145 ... 86 0
9delta 0 0 ... 95 5
10gamma 0 0 ... 14 0
11endothelial 1 0 ... 7 0
12activated_stellate 49 1 ... 17 0
13quiescent_stellate 4 0 ... 1 0
14macrophage 0 0 ... 1 0
15mast 0 0 ... 1 0
16
17[11 rows x 16 columns]
我们观察到PSC(胰腺星状细胞)细胞实际上只是不一致地注释并正确地映射到“激活的星状细胞”上。
此外，很高兴看到“间充质”和“间充质”细胞都属于同一类别。但是，这个类别又是“activated_stellate”，而且可能是错误的。这就是我们说的，算法只能接近真相，而不能定义真相。
可视化分布的批次
通常，批量对应的是想要比较的实验。Scanpy提供了方便的可视化可能性，主要有
a density plot
a partial visualization of a subset of categories/groups in an emnbedding
1sc.tl.embedding_density(adata_concat, groupby='batch')
2sc.pl.embedding_density(adata_concat, groupby='batch')
1for batch in ['1', '2', '3']:
2 sc.pl.umap(adata_concat, color='batch', groups=[batch])
References
[1] BBKNN: fast batch alignment of single cell transcriptomes : https://academic.oup.com/bioinformatics/article/36/3/964/5545955
[2] integrating-data-using-ingest : https://scanpy-tutorials.readthedocs.io/en/latest/integrating-data-using-ingest.html