美文网首页单细胞分析单细胞分析之Scanpypython单细胞分析
单细胞转录组数据分析|| scanpy教程:使用ingest和B

单细胞转录组数据分析|| scanpy教程:使用ingest和B

作者: 周运来就是我 | 来源:发表于2020-03-27 12:14 被阅读0次

    随着单细胞技术的成熟,测序成本的降低,单细胞的数据量和样本量也日益增长。我们知道单细胞转录组的一个主要应用就是解释细胞的异质性,那么,不同器官,不同测序平台,不同物种之间的单细胞数据何如整合分析呢?特别是在单细胞的数据维度这么高的前提下,显然传统的基于回归的方法已经不适用了。于是出现了一批单细胞整合分析的工具,它们大多数是在R生态条件下的。如:

    在我们理解单细胞数据的时候一张cell X gene 的大表不能离开我们的脑海。

    adata.to_df()
    Out[24]: 
                        RP11-34P13.3  FAM138A  ...  AC213203.1  FAM231B
    AAACCCAAGCGTATGG-1           0.0      0.0  ...         0.0      0.0
    AAACCCAGTCCTACAA-1           0.0      0.0  ...         0.0      0.0
    AAACCCATCACCTCAC-1           0.0      0.0  ...         0.0      0.0
    AAACGCTAGGGCATGT-1           0.0      0.0  ...         0.0      0.0
    AAACGCTGTAGGTACG-1           0.0      0.0  ...         0.0      0.0
                             ...      ...  ...         ...      ...
    TTTGTTGCAGGTACGA-1           0.0      0.0  ...         0.0      0.0
    TTTGTTGCAGTCTCTC-1           0.0      0.0  ...         0.0      0.0
    TTTGTTGGTAATTAGG-1           0.0      0.0  ...         0.0      0.0
    TTTGTTGTCCTTGGAA-1           0.0      0.0  ...         0.0      0.0
    TTTGTTGTCGCACGAC-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应该也是可以用的。

    先来看ingest,通过投射到参考数据上的PCA(或备用模型)上,将一个adata的嵌入和注释与一个参考数据集adata_ref集成在一起。该函数使用knn分类器来映射标签,使用UMAP来映射嵌入。

    再来看看bbknn是一个快速和直观的批处理效果去除工具,可以直接在scanpy工作流中使用。它是scanpy.api.pp.neighbors()的替代方法,这两个函数都创建了一个邻居图,以便后续在集群、伪时间和UMAP可视化中使用。标准方法首先确定整个数据结构中每个单元的k个最近邻,然后将候选单元转换为指数相关的连接,然后作为进一步分析的基础。

    那么我们就来看一下在scanpy的实现吧。

    import scanpy as sc
    import pandas as pd
    import seaborn as sns
    import sklearn
    import sys
    import scipy
    import bbknn
    
    sc.settings.verbosity = 1             # verbosity: errors (0), warnings (1), info (2), hints (3)
    sc.logging.print_versions()
    sc.settings.set_figure_params(dpi=80, frameon=False, figsize=(3, 3))
    
    scanpy==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 注释
    adata_ref = sc.datasets.pbmc3k_processed()  # this is an earlier version of the dataset from the pbmc3k tutorial
    
    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'
    

    我们一次看看以下参考数据集都有哪些内容:

    adata_ref.obs
    Out[9]: 
                      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_ref.var
    Out[10]: 
             n_cells
    index           
    TNFRSF4      155
    CPSF3L       202
    ATAD3C         9
    C1orf86      501
    RER1         608
             ...
    ICOSLG        34
    SUMO3        570
    SLC19A1       31
    S100B         94
    PRMT2        588
    
    [1838 rows x 1 columns]
    
    adata_ref.uns['louvain_colors']
    Out[14]: 
    array(['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b',
           '#e377c2', '#bcbd22'], dtype='<U7')
    
     adata_ref.obsm
    Out[16]: AxisArrays with keys: X_pca, X_tsne, X_umap, X_draw_graph_fr
    
    adata_ref.obsm['X_umap']
    Out[17]: 
    array([[ 1.35285574,  2.26612719],
           [-0.47802448,  7.87730423],
           [ 2.16588875, -0.24481226],
           ...,
           [ 0.34670979,  8.34967798],
           [ 0.19864146,  9.56698797],
           [ 2.62803322,  0.36722543]])
    

    有没有再次理解AnnData 这个对象的数据结构呢?

    可以看到在这个数据集中降维聚类都是做过的,我们可以画个图看看:

    sc.pl.umap(adata_ref, color='louvain')
    

    接下来我们看看要预测的数据集是怎样的。

    adata = sc.datasets.pbmc68k_reduced()
    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'
    

    可见它也是降维聚类过的了。

    sc.pl.umap(adata, color='louvain')
    

    这个数据集并没有得到细胞类型的定义。

    构建注释数据结构:

    var_names = adata_ref.var_names.intersection(adata.var_names) # 取交集
    adata_ref = adata_ref[:, var_names]
    adata = adata[:, var_names]
    
    sc.pp.pca(adata_ref)
    sc.pp.neighbors(adata_ref)
    sc.tl.umap(adata_ref)
    sc.tl.leiden(adata_ref)# 新的聚类方法
    sc.pl.umap(adata_ref, color=['louvain','leiden'])
    adata_ref
    
    AnnData object with n_obs × n_vars = 2638 × 208 
        obs: 'n_genes', 'percent_mito', 'n_counts', 'louvain', 'leiden'
        var: 'n_cells'
        uns: 'draw_graph', 'louvain', 'louvain_colors', 'neighbors', 'pca', 'rank_genes_groups', 'umap', 'leiden', 'leiden_colors'
        obsm: 'X_pca', 'X_tsne', 'X_umap', 'X_draw_graph_fr'
        varm: 'PCs'
    
    sc.pp.pca(adata)
    sc.pp.neighbors(adata)
    sc.tl.umap(adata)
    sc.tl.leiden(adata)
    sc.pl.umap(adata, color=['louvain','leiden'])
    adata
    
    AnnData object with n_obs × n_vars = 700 × 208 
        obs: 'bulk_labels', 'n_genes', 'percent_mito', 'n_counts', 'S_score', 'G2M_score', 'phase', 'louvain', 'leiden'
        var: 'n_counts', 'means', 'dispersions', 'dispersions_norm', 'highly_variable'
        uns: 'bulk_labels_colors', 'louvain', 'louvain_colors', 'neighbors', 'pca', 'rank_genes_groups', 'umap', 'leiden', 'leiden_colors'
        obsm: 'X_pca', 'X_umap'
        varm: 'PCs'
    

    用ingest来做细胞注释吧。

    sc.tl.ingest(adata, adata_ref, obs='louvain')
    adata.uns['louvain_colors'] = adata_ref.uns['louvain_colors']  # fix colors
    

    我们来看看sc.tl.ingest的帮助文档:

    Help on function ingest in module scanpy.tools._ingest:
    
    ingest(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)
        Map labels and embeddings from reference data to new data.
        
        :tutorial:`integrating-data-using-ingest`
        
        Integrates embeddings and annotations of an `adata` with a reference dataset
        `adata_ref` through projecting on a PCA (or alternate
        model) that has been fitted on the reference data. The function uses a knn
        classifier for mapping labels and the UMAP package [McInnes18]_ for mapping
        the embeddings.
        
        .. note::
        
            We refer to this *asymmetric* dataset integration as *ingesting*
            annotations from reference data to new data. This is different from
            learning a joint representation that integrates both datasets in an
            unbiased way, as CCA (e.g. in Seurat) or a conditional VAE (e.g. in
            scVI) would do.
        
        You need to run :func:`~scanpy.pp.neighbors` on `adata_ref` before
        passing it.
        
        Parameters
        ----------
        adata
            The annotated data matrix of shape `n_obs` × `n_vars`. Rows correspond
            to cells and columns to genes. This is the dataset without labels and
            embeddings.
        adata_ref
            The annotated data matrix of shape `n_obs` × `n_vars`. Rows correspond
            to cells and columns to genes.
            Variables (`n_vars` and `var_names`) of `adata_ref` should be the same
            as in `adata`.
            This is the dataset with labels and embeddings
            which need to be mapped to `adata`.
        obs
            Labels' keys in `adata_ref.obs` which need to be mapped to `adata.obs`
            (inferred for observation of `adata`).
        embedding_method
            Embeddings in `adata_ref` which need to be mapped to `adata`.
            The only supported values are 'umap' and 'pca'.
        labeling_method
            The method to map labels in `adata_ref.obs` to `adata.obs`.
            The only supported value is 'knn'.
        inplace
            Only works if `return_joint=False`.
            Add labels and embeddings to the passed `adata` (if `True`)
            or return a copy of `adata` with mapped embeddings and labels.
        
        Returns
        -------
        * if `inplace=False` returns a copy of `adata`
          with mapped embeddings and labels in `obsm` and `obs` correspondingly
        * if `inplace=True` returns `None` and updates `adata.obsm` and `adata.obs`
          with mapped embeddings and labels
        
        Example
        -------
        Call sequence:
        
    import scanpy as sc
    sc.pp.neighbors(adata_ref)
    sc.tl.umap(adata_ref)
    sc.tl.ingest(adata, adata_ref, obs='cell_type')
        
        .. _ingest PBMC tutorial: https://scanpy-tutorials.readthedocs.io/en/latest/integrating-pbmcs-using-ingest.html
        .. _ingest Pancreas tutorial: https://scanpy-tutorials.readthedocs.io/en/latest/integrating-pancreas-using-ingest.html
    

    通过比较‘bulk_label’注释和‘louvain’注释,我们发现数据被合理地映射,只有树突细胞的注释似乎是含糊不清的,在adata中可能已经是模糊的了。我们来对adata做进一步的处理。

    adata_concat = adata_ref.concatenate(adata, batch_categories=['ref', 'new'])
    adata_concat.obs.louvain  = adata_concat.obs.louvain.astype('category')
    adata_concat.obs.louvain.cat.reorder_categories(adata_ref.obs.louvain.cat.categories, inplace=True)  # fix category ordering
    adata_concat.uns['louvain_colors'] = adata_ref.uns['louvain_colors']  # fix category colors
    adata_concat
    sc.pl.umap(adata_concat, color=['batch', 'louvain'])
    
    
    AnnData object with n_obs × n_vars = 3338 × 208 
        obs: 'G2M_score', 'S_score', 'batch', 'bulk_labels', 'leiden', 'louvain', 'n_counts', 'n_genes', 'percent_mito', 'phase'
        var: 'n_cells-ref', 'n_counts-new', 'means-new', 'dispersions-new', 'dispersions_norm-new', 'highly_variable-new'
        obsm: 'X_pca', 'X_umap'
    

    虽然在单核细胞和树突状细胞簇中似乎存在一些批处理效应,但在其他方面,新数据被绘制得相对均匀。
    巨核细胞只存在于adata_ref中,没有来自adata映射的单元格。如果交换参考数据和查询数据,巨核细胞不再作为单独的集群出现。这是一个极端的情况,因为参考数据非常小;但是,人们应该始终质疑参考数据是否包含足够的生物变异,以便有意义地容纳查询数据。

    使用BBKNN整合
    sc.tl.pca(adata_concat)
    sc.external.pp.bbknn(adata_concat, batch_key='batch')  # running bbknn 1.3.6
    sc.tl.umap(adata_concat)
    sc.pl.umap(adata_concat, color=['batch', 'louvain'])
    
    adata_concat
    Out[45]: 
    AnnData object with n_obs × n_vars = 3338 × 208 
        obs: 'G2M_score', 'S_score', 'batch', 'bulk_labels', 'leiden', 'louvain', 'n_counts', 'n_genes', 'percent_mito', 'phase'
        var: 'n_cells-ref', 'n_counts-new', 'means-new', 'dispersions-new', 'dispersions_norm-new', 'highly_variable-new'
        uns: 'batch_colors', 'louvain_colors', 'pca', 'neighbors', 'umap'
        obsm: 'X_pca', 'X_umap'
        varm: 'PCs'
    

    BBKNN并不维持巨核细胞簇。然而,它似乎更均匀地混合细胞。

    一个例子使用BBKNN整合数据的例子

    以下数据已在scGen论文[Lotfollahi19]中使用。点击pancreas下载数据。

    它包含了来自4个不同研究(Segerstolpe16, Baron16, Wang16, Muraro16)的人类胰腺数据,这些数据在单细胞数据集集成的开创性论文(Butler18, Haghverdi18)中被使用过,并在此后多次被使用。

    h5ad = 'E:\\learnscanpy\\data\\objects-pancreas\\pancreas.h5ad'
    adata_all = sc.read_h5ad(h5ad)
    adata_all
    
    AnnData object with n_obs × n_vars = 14693 × 2448 
        obs: 'celltype', 'sample', 'n_genes', 'batch', 'n_counts', 'louvain'
        var: 'n_cells-0', 'n_cells-1', 'n_cells-2', 'n_cells-3'
        uns: 'celltype_colors', 'louvain', 'neighbors', 'pca', 'sample_colors'
        obsm: 'X_pca', 'X_umap'
        varm: 'PCs'
    
     counts = adata_all.obs.celltype.value_counts()
    counts
    Out[173]: 
    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
    
    adata_all.obs
    Out[171]: 
                                  celltype sample  ...      n_counts louvain
    index                                          ...                      
    human1_lib1.final_cell_0001-0   acinar  Baron  ...  2.241100e+04       2
    human1_lib1.final_cell_0002-0   acinar  Baron  ...  2.794900e+04       2
    human1_lib1.final_cell_0003-0   acinar  Baron  ...  1.689200e+04       2
    human1_lib1.final_cell_0004-0   acinar  Baron  ...  1.929900e+04       2
    human1_lib1.final_cell_0005-0   acinar  Baron  ...  1.506700e+04       2
                                   ...    ...  ...           ...     ...
    reads.29499-3                   ductal   Wang  ...  1.056558e+06      10
    reads.29500-3                   ductal   Wang  ...  9.926309e+05      10
    reads.29501-3                     beta   Wang  ...  1.751338e+06      10
    reads.29502-3                  dropped   Wang  ...  2.163764e+06      10
    reads.29503-3                     beta   Wang  ...  2.038979e+06      10
    
    [14693 rows x 6 columns]
    

    可以看出这个数据集已经降维聚类好了,所以我们可以可视化一下:

    sc.pl.umap(adata_all,color=['sample', 'celltype','louvain'])
    

    样本之间的批次很严重啊。

    去掉细胞数较小的小群,

    minority_classes = counts.index[-5:].tolist()        # get the minority classes
    
    # ['schwann', 't_cell', 'MHC class II', 'unclear', 'unclassified']
    
    adata_all = adata_all[                               # actually subset
        ~adata_all.obs.celltype.isin(minority_classes)]
    adata_all.obs.celltype.cat.reorder_categories(       # reorder according to abundance
        counts.index[:-5].tolist(), inplace=True)
    
    adata_all.obs.celltype.value_counts()
    Out[182]: 
    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
    

    进行pca降维和umap降维:

    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)
    

    下面我们使用BBKNN来整合数据:

    sc.external.pp.bbknn(adata_all, batch_key='batch')
    sc.tl.umap(adata_all)
    adata_all
    sc.pl.umap(adata_all, color=['sample','batch', 'celltype'])
    

    果然要比原始的数据好多了。但是改变的是什么?

    AnnData object with n_obs × n_vars = 14662 × 2448 
        obs: 'celltype', 'sample', 'n_genes', 'batch', 'n_counts', 'louvain'
        var: 'n_cells-0', 'n_cells-1', 'n_cells-2', 'n_cells-3'
        uns: 'celltype_colors', 'louvain', 'neighbors', 'pca', 'sample_colors', 'louvain_colors', 'umap', 'batch_colors'
        obsm: 'X_pca', 'X_umap'
        varm: 'PCs'
    

    如果想对其中某个样本进行单独的注释,可以用上面提到的ingest。选择一个参考批次来训练模型和建立邻域图(这里是一个PCA),并分离出所有其他批次。

    adata_ref = adata_all[adata_all.obs.batch == '0']
    adata_ref
    Out[191]: 
    View of AnnData object with n_obs × n_vars = 8549 × 2448 
        obs: 'celltype', 'sample', 'n_genes', 'batch', 'n_counts', 'louvain'
        var: 'n_cells-0', 'n_cells-1', 'n_cells-2', 'n_cells-3'
        uns: 'celltype_colors', 'louvain', 'neighbors', 'pca', 'sample_colors', 'louvain_colors', 'umap', 'batch_colors'
        obsm: 'X_pca', 'X_umap'
        varm: 'PCs'
    
    sc.pp.pca(adata_ref)
    sc.pp.neighbors(adata_ref)
    sc.tl.umap(adata_ref)
    
    adata_ref
    Out[197]: 
    AnnData object with n_obs × n_vars = 8549 × 2448 
        obs: 'celltype', 'sample', 'n_genes', 'batch', 'n_counts', 'louvain'
        var: 'n_cells-0', 'n_cells-1', 'n_cells-2', 'n_cells-3'
        uns: 'celltype_colors', 'louvain', 'neighbors', 'pca', 'sample_colors', 'louvain_colors', 'umap', 'batch_colors'
        obsm: 'X_pca', 'X_umap'
        varm: 'PCs'
    
    sc.pl.umap(adata_ref, color='celltype')
    

    选取数据集用ingest于adata_ref进行mapping:

    adatas = [adata_all[adata_all.obs.batch == i].copy() for i in ['1', '2', '3']]
    sc.settings.verbosity = 2  # a bit more logging
    for iadata, adata in enumerate(adatas):
        print(f'... integrating batch {iadata+1}')
        adata.obs['celltype_orig'] = adata.obs.celltype  # save the original cell type
        sc.tl.ingest(adata, adata_ref, obs='celltype')
    
    integrating batch 1
    running ingest
        finished (0:00:08)
    integrating batch 2
    running ingest
        finished (0:00:06)
    integrating batch 3
    running ingest
        finished (0:00:03)
    
    adata_concat = adata_ref.concatenate(adatas)
    adata_concat
    
    Out[200]: 
    AnnData object with n_obs × n_vars = 14662 × 2448 
        obs: 'batch', 'celltype', 'celltype_orig', 'louvain', 'n_counts', 'n_genes', 'sample'
        var: 'n_cells-0', 'n_cells-1', 'n_cells-2', 'n_cells-3'
        obsm: 'X_pca', 'X_umap'
    
    adata_concat.obs.celltype = adata_concat.obs.celltype.astype('category')
    adata_concat.obs.celltype.cat.reorder_categories(adata_ref.obs.celltype.cat.categories, inplace=True)  # fix category ordering
    adata_concat.uns['celltype_colors'] = adata_ref.uns['celltype_colors']  # fix category coloring
    
    sc.pl.umap(adata_concat, color=['celltype_orig','batch', 'celltype'])
    

    与BBKNN的结果相比,这是以一种更加明显的方式保持分群。如果已经观察到一个想要的连续结构(例如在造血数据集中),摄取允许容易地维持这个结构。

    一致性评估
    adata_query = adata_concat[adata_concat.obs.batch.isin(['1', '2', '3'])]
    
    View of AnnData object with n_obs × n_vars = 6113 × 2448 
        obs: 'batch', 'celltype', 'celltype_orig', 'louvain', 'n_counts', 'n_genes', 'sample'
        var: 'n_cells-0', 'n_cells-1', 'n_cells-2', 'n_cells-3'
        uns: 'celltype_colors', 'celltype_orig_colors', 'batch_colors'
        obsm: 'X_pca', 'X_umap'
    
    sc.pl.umap(
        adata_query, color=['batch', 'celltype', 'celltype_orig'], wspace=0.4)
    

    这个结果依然不能很好的反映一致性,让我们首先关注与参考保守的细胞类型,以简化混淆矩阵的reads。

    obs_query = adata_query.obs
    conserved_categories = obs_query.celltype.cat.categories.intersection(obs_query.celltype_orig.cat.categories)  # intersected categories
    obs_query_conserved = obs_query.loc[obs_query.celltype.isin(conserved_categories) & obs_query.celltype_orig.isin(conserved_categories)]  # intersect categories
    obs_query_conserved.celltype.cat.remove_unused_categories(inplace=True)  # remove unused categoriyes
    obs_query_conserved.celltype_orig.cat.remove_unused_categories(inplace=True)  # remove unused categoriyes
    obs_query_conserved.celltype_orig.cat.reorder_categories(obs_query_conserved.celltype.cat.categories, inplace=True)  # fix category ordering
    
    obs_query_conserved
    Out[214]: 
                    batch     celltype celltype_orig  ...      n_counts  n_genes  sample
    D28.1_1-1-1         1        alpha         alpha  ...  2.322583e+04     5448  Muraro
    D28.1_13-1-1        1       ductal        ductal  ...  2.334263e+04     5911  Muraro
    D28.1_15-1-1        1        alpha         alpha  ...  2.713471e+04     5918  Muraro
    D28.1_17-1-1        1        alpha         alpha  ...  1.581207e+04     4522  Muraro
    D28.1_2-1-1         1  endothelial   endothelial  ...  3.173151e+04     6464  Muraro
                  ...          ...           ...  ...           ...      ...     ...
    reads.29498-3-3     3       ductal        ductal  ...  1.362606e+06    19950    Wang
    reads.29499-3-3     3       ductal        ductal  ...  1.056558e+06    19950    Wang
    reads.29500-3-3     3       ductal        ductal  ...  9.926309e+05    19950    Wang
    reads.29501-3-3     3         beta          beta  ...  1.751338e+06    19950    Wang
    reads.29503-3-3     3         beta          beta  ...  2.038979e+06    19950    Wang
    
    pd.crosstab(obs_query_conserved.celltype, obs_query_conserved.celltype_orig)
    Out[215]: 
    celltype_orig  alpha  beta  ductal  acinar  delta  gamma  endothelial  mast
    celltype                                                                   
    alpha           1819     3       7       0      1     20            0     5
    beta              49   803       3       1     10     26            0     0
    ductal             8     5     693     263      0      0            0     0
    acinar             1     3       2     145      0      3            0     0
    delta              5     4       0       0    305     73            0     0
    gamma              1     5       0       0      0    194            0     0
    endothelial        2     0       0       0      0      0           36     0
    mast               0     0       1       0      0      0            0     2
    

    总的来说,保守的细胞类型也如预期的那样被映射。主要的例外是原始注释中出现的一些腺泡细胞。然而,已经观察到参考数据同时具有腺泡和导管细胞,这就解释了差异,并指出了初始注释中潜在的不一致性。

    现在让我们继续看看所有的细胞类型。

     pd.crosstab(adata_query.obs.celltype, adata_query.obs.celltype_orig)
    Out[216]: 
    celltype_orig       PSC  acinar  ...  not applicable  unclassified endocrine
    celltype                         ...                                        
    alpha                 0       0  ...             304                      11
    beta                  0       1  ...             522                      24
    ductal                0     263  ...             106                       1
    acinar                0     145  ...              86                       0
    delta                 0       0  ...              95                       5
    gamma                 0       0  ...              14                       0
    endothelial           1       0  ...               7                       0
    activated_stellate   49       1  ...              17                       0
    quiescent_stellate    4       0  ...               1                       0
    macrophage            0       0  ...               1                       0
    mast                  0       0  ...               1                       0
    
    [11 rows x 16 columns]
    

    我们观察到PSC(胰腺星状细胞)细胞实际上只是不一致地注释并正确地映射到“激活的星状细胞”上。
    此外,很高兴看到“间充质”和“间充质”细胞都属于同一类别。但是,这个类别又是“activated_stellate”,而且可能是错误的。这就是我们说的,算法只能接近真相,而不能定义真相。

    可视化分布的批次

    通常,批量对应的是想要比较的实验。Scanpy提供了方便的可视化可能性,主要有

    • a density plot
    • a partial visualization of a subset of categories/groups in an emnbedding
    sc.tl.embedding_density(adata_concat, groupby='batch')
    sc.pl.embedding_density(adata_concat, groupby='batch')
    
    for batch in ['1', '2', '3']:
        sc.pl.umap(adata_concat, color='batch', groups=[batch])
    


    BBKNN: fast batch alignment of single cell transcriptomes
    integrating-data-using-ingest

    相关文章

      网友评论

        本文标题:单细胞转录组数据分析|| scanpy教程:使用ingest和B

        本文链接:https://www.haomeiwen.com/subject/gkxkuhtx.html