Tutorial:Single-cell RNA-seq: Annotation: Celltype annotation migration(mapping) with TOSICA
0
0
Entering edit mode
3 months ago
Julia Ma ▴ 120

We know that when all samples cannot be obtained at the same time, it would be desirable to classify the cell types on the first batch of data and use them to annotate the data obtained later or to be obtained in the future with the same standard, without the need to processing and mapping them together again.

So migration(mapping) the reference cell annotation is necessary. This tutorial focuses on how to migration(mapping) the cell annotation from reference scRNA-seq atlas to new scRNA-seq data.

Paper: Transformer for one stop interpretable cell type annotation

Code: https://github.com/JackieHanLab/TOSICA

Colab_Reproducibility: https://colab.research.google.com/drive/1BjPEG-kLAgicP8iQvtVtpzzbIOmk1X23?usp=sharing

import omicverse as ov
import scanpy as sc
ov.utils.ov_plot_set()

Loading data

  • demo_train.h5ad: Braon(GSE84133) and Muraro(GSE85241)
  • demo_test.h5ad : xin(GSE81608), segerstolpe(E-MTAB-5061) and Lawlor(GSE86473)

They can be downloaded at https://figshare.com/projects/TOSICA_demo/158489.

ref_adata = sc.read('demo_train.h5ad')
ref_adata = ref_adata[:,ref_adata.var_names]
print(ref_adata)
print(ref_adata.obs.Celltype.value_counts())

View of AnnData object with n_obs × n_vars = 10600 × 3000
    obs: 'Celltype'
    var: 'Gene Symbol'
alpha          3136
beta           2966
ductal         1290
acinar         1144
delta           793
PSC             524
PP              356
endothelial     273
macrophage       52
mast             25
epsilon          21
schwann          13
t_cell            7
Name: Celltype, dtype: int64



query_adata = sc.read('demo_test.h5ad')
query_adata = query_adata[:,ref_adata.var_names]
print(query_adata)
print(query_adata.obs.Celltype.value_counts())

View of AnnData object with n_obs × n_vars = 4218 × 3000
    obs: 'Celltype'
    var: 'Gene Symbol'
alpha           2011
beta            1006
ductal           414
PP               282
acinar           209
delta            188
PSC               73
endothelial       16
epsilon            7
mast               7
MHC class II       5
Name: Celltype, dtype: int64

We need to select the same gene training and predicting the celltype

ref_adata.var_names_make_unique()
query_adata.var_names_make_unique()
ret_gene=list(set(query_adata.var_names) & set(ref_adata.var_names))
len(ret_gene)

query_adata=query_adata[:,ret_gene]
ref_adata=ref_adata[:,ret_gene]

print(f"The max of ref_adata is {ref_adata.X.max()}, query_data is {query_adata.X.max()}",)

The max of ref_adata is 8.72524356842041, query_data is 9.2676362991333

By comparing the maximum values of the two data, we can see that the data has been normalised sc.pp.normalize_total and logarithmised sc.pp.log1p. The same treatment is applied to the data when we use our own data for analysis.

Download Genesets

Here, we need to download the genesets as pathway at first. You can use ov.utils.download_tosica_gmt() to download automatically or manual download from:

Initialisation the TOSICA model

We first need to train the TOSICA model on the REFERENCE dataset, where omicverse provides a simple class pyTOSICA, and all subsequent operations can be done with pyTOSICA. We need to set the parameters for model initialisation.

  • adata: the reference adata object
  • gmt_path: default pre-prepared mask or path to .gmt files. you can use ov.utils.download_tosica_gmt() to obtain the genesets
  • depth: the depth of transformer model, When it is set to 2, a memory leak may occur
  • label_name: the reference key of celltype in adata.obs
  • project_path: the save path of TOSICA model
  • batch_size: indicates the number of cells passed to the programme for training in a single pass

    tosica_obj=ov.single.pyTOSICA(adata=ref_adata,

                                gmt_path='genesets/GO_bp.gmt', depth=1,
                                label_name='Celltype',
                                project_path='hGOBP_demo',
                                batch_size=8)
    

    cuda:0 Mask loaded!

Training the TOSICA model

There're 4 arguments to set when training the TOSICA model.

  • pre_weights: The path of the pre-trained weights.
  • lr: The learning rate.
  • epochs: The number of epochs.
  • lrf: The learning rate of the last layer.

    tosica_obj.train(epochs=5)

The model can be loaded from project_path automatically.

tosica_obj.load()

Model loaded!

Update with query

new_adata=tosica_obj.predicted(pre_adata=query_adata)

0
4218

Visualize the reference and mapping

We first compute the lower dimensional space of query_data, where we use omicverse's preprocessing method as well as the mde method for dimensionality reduction

To visualize the PCA's embeddings, we use the pymde package wrapper in omicverse. This is an alternative to UMAP that is GPU-accelerated.

ov.pp.scale(query_adata)
ov.pp.pca(query_adata,layer='scaled',n_pcs=50)
sc.pp.neighbors(query_adata, n_neighbors=15, n_pcs=50,
               use_rep='scaled|original|X_pca')
query_adata.obsm["X_mde"] = ov.utils.mde(query_adata.obsm["scaled|original|X_pca"])
query_adata

computing neighbors
    finished: added to `.uns['neighbors']`
    `.obsp['distances']`, distances for each pair of neighbors
    `.obsp['connectivities']`, weighted adjacency matrix (0:00:02)





AnnData object with n_obs × n_vars = 4218 × 3000
    obs: 'Celltype'
    var: 'Gene Symbol'
    uns: 'scaled|original|pca_var_ratios', 'scaled|original|cum_sum_eigenvalues', 'neighbors'
    obsm: 'scaled|original|X_pca', 'X_mde'
    varm: 'scaled|original|pca_loadings'
    layers: 'scaled', 'lognorm'
    obsp: 'distances', 'connectivities'

Since new_adata and query_adata have the same cells, their low-dimensional spaces are also the same. So we proceed directly to the assignment operation.

new_adata.obsm=query_adata[new_adata.obs.index].obsm.copy()
new_adata.obsp=query_adata[new_adata.obs.index].obsp.copy()
new_adata



AnnData object with n_obs × n_vars = 4218 × 299
    obs: 'Prediction', 'Probability', 'Celltype'
    obsm: 'scaled|original|X_pca', 'X_mde'
    obsp: 'distances', 'connectivities'

Since the predicted cell types are not exactly the same as the original cell types, the colours are not exactly the same. For the visualisation effect, we manually set the colour of the predicted cell type with the original cell type.

import numpy as np
col = np.array([
"#98DF8A","#E41A1C" ,"#377EB8", "#4DAF4A" ,"#984EA3" ,"#FF7F00" ,"#FFFF33" ,"#A65628" ,"#F781BF" ,"#999999","#1F77B4","#FF7F0E","#279E68","#FF9896"
]).astype('<U7')

celltype = ("alpha","beta","ductal","acinar","delta","PP","PSC","endothelial","epsilon","mast","macrophage","schwann",'t_cell')
new_adata.obs['Prediction'] = new_adata.obs['Prediction'].astype('category')
new_adata.obs['Prediction'] = new_adata.obs['Prediction'].cat.reorder_categories(list(celltype))
new_adata.uns['Prediction_colors'] = col[1:]

celltype = ("MHC class II","alpha","beta","ductal","acinar","delta","PP","PSC","endothelial","epsilon","mast")
new_adata.obs['Celltype'] = new_adata.obs['Celltype'].astype('category')
new_adata.obs['Celltype'] = new_adata.obs['Celltype'].cat.reorder_categories(list(celltype))
new_adata.uns['Celltype_colors'] = col[:11]

sc.pl.embedding(
    new_adata,
    basis="X_mde",
    color=['Celltype', 'Prediction'],
    frameon=False,
    #ncols=1,
    wspace=0.5,
    #palette=ov.utils.pyomic_palette()[11:],
    show=False,
)



[<AxesSubplot: title={'center': 'Celltype'}, xlabel='X_mde1', ylabel='X_mde2'>,
 <AxesSubplot: title={'center': 'Prediction'}, xlabel='X_mde1', ylabel='X_mde2'>]

enter image description here

Pathway attention

TOSICA has another special feature, which is the ability to computationally use self-attention mechanisms to find pathways associated with cell types. Here we demonstrate the approach of this downstream analysis.

We first need to filter out the predicted types of cells with cell counts less than 5.

cell_idx=new_adata.obs['Prediction'].value_counts()[new_adata.obs['Prediction'].value_counts()<5].index
new_adata=new_adata[~new_adata.obs['Prediction'].isin(cell_idx)]

We then used sc.tl.rank_genes_groups to calculate the differential pathways with the highest attention across cell types. This differential pathway is derived from the gmt genesets used for the previous calculation.

sc.tl.rank_genes_groups(new_adata, 'Prediction', method='wilcoxon')

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:00)




sc.pl.rank_genes_groups_dotplot(new_adata,
                                n_genes=3,standard_scale='var',)

enter image description here

If you want to get the cell-specific pathway, you can use sc.get.rank_genes_groups_df to get.

For example, we would like to obtain the pathway with the highest attention for the cell type PP

degs = sc.get.rank_genes_groups_df(new_adata, group='PP', key='rank_genes_groups',
                                            pval_cutoff=0.05)
degs.head()
<style scoped="">.dataframe tbody tr th:only-of-type { vertical-align: middle; } .dataframe tbody tr th { vertical-align: top; } .dataframe thead th { text-align: right; }</style>
names scores logfoldchanges pvals pvals_adj
0 GOBP_REGULATION_OF_MUSCLE_SYSTEM_PROCESS 27.348158 0.008944 1.135759e-164 3.395921e-162
1 GOBP_RESPONSE_TO_ALCOHOL 26.803308 0.008620 2.956726e-158 4.420305e-156
2 GOBP_REGULATION_OF_CELL_JUNCTION_ASSEMBLY 26.634876 0.008843 2.679320e-156 2.670389e-154
3 GOBP_REGULATION_OF_BLOOD_CIRCULATION 26.348188 0.008339 5.383820e-153 4.024406e-151
4 GOBP_MITOTIC_NUCLEAR_DIVISION 25.993446 0.007867 5.873348e-149 3.512262e-147
sc.pl.embedding(
    new_adata,
    basis="X_mde",
    color=['Prediction','GOBP_REGULATION_OF_MUSCLE_SYSTEM_PROCESS'],
    frameon=False,
    #ncols=1,
    wspace=0.5,
    #palette=ov.utils.pyomic_palette()[11:],
    show=False,
)



[<AxesSubplot: title={'center': 'Prediction'}, xlabel='X_mde1', ylabel='X_mde2'>,
 <AxesSubplot: title={'center': 'GOBP_REGULATION_OF_MUSCLE_SYSTEM_PROCESS'}, xlabel='X_mde1', ylabel='X_mde2'>]

enter image description here

If you call omciverse to complete a TOSICA analysis, don't forget to cite the following literature:

@article{pmid:36641532,
journal = {Nature communications},
doi = {10.1038/s41467-023-35923-4},
issn = {2041-1723},
number = {1},
pmid = {36641532},
pmcid = {PMC9840170},
address = {England},
title = {Transformer for one stop interpretable cell type annotation},
volume = {14},
author = {Chen, Jiawei and Xu, Hao and Tao, Wanyu and Chen, Zhaoxiong and Zhao, Yuxuan and Han, Jing-Dong J},
note = {[Online; accessed 2023-07-18]},
pages = {223},
date = {2023-01-14},
year = {2023},
month = {1},
day = {14},
}


@misc{doi:10.1101/2023.06.06.543913,
doi = {10.1101/2023.06.06.543913},
publisher = {Cold Spring Harbor Laboratory},
title = {OmicVerse: A single pipeline for exploring the entire transcriptome universe},
author = {Zeng, Zehua and Ma, Yuqing and Hu, Lei and Xiong, Yuanyan and Du, Hongwu},
note = {[Online; accessed 2023-07-18]},
date = {2023-06-08},
year = {2023},
month = {6},
day = {8},
}
scRNA-seq TOSICA • 583 views
ADD COMMENT
0
Entering edit mode

Did you not read the comments on your previous posts? Stop copy-pasting content verbatim from the Read The Docs site. Please consider this a final warning.

ADD REPLY
0
Entering edit mode

Hi, I was wondering why I can't copy it? We are part of the omicverse development team, the copyright of the tutorial is in our hands, and putting the tutorial on biostars will help more people I think.

ADD REPLY
1
Entering edit mode

Help people how? The tutorial can be reached via google. People usually ask specific questions, so simply dumping a wiki’s worth of text will most likely not do anyone a favor.

ADD REPLY
0
Entering edit mode

Why hasn't the Biostars section, labeled as Tutorials, been removed if its content is supposed to be pertinent to a specific problem? Additionally, the wiki encyclopedia serves as an introductory resource, whereas tutorials ought to provide comprehensive guidance, covering everything from data to code. Could you shed light on the reasoning behind retaining this section in Biostars?

ADD REPLY
0
Entering edit mode

Why hasn't the Biostars section, labeled as Tutorials, been removed if its content is supposed to be pertinent to a specific problem?

Because it's not centered around a problem, it's centered around an approach. Your Read The Docs approach is much better than posting a tutorial just on Biostars - which is what that kind of post is meant to do.

Additionally, the wiki encyclopedia serves as an introductory resource, whereas tutorials ought to provide comprehensive guidance, covering everything from data to code.

Your collaborator defeats this point by copy-pasting the content verbatim. The content here is not an extension of or complementary to the wiki, it's a copy of it. A copy without a link to the original. That is a few steps away from spam - add disclaimers when you're promoting your product. Unsolicited sneaky promos are spam.

Your collaborator was asked to change their approach to this copy-paste thing and they refused to listen. They have been warned again. You're welcome to create one Tool type post for omicverse and leave it at that. You're also welcome to use Biostars as (one of) the official means of supporting omicverse related problems. You can even create custom tutorials or extend existing tutorials or address special cases. This copy-paste promotional work is the problematic part.

ADD REPLY

Login before adding your answer.

Traffic: 2670 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6