3.基于python的scRNA-seq细胞类型注释-自动注释

最新推荐文章于 2024-06-17 17:35:06 发布

tzc_fly

最新推荐文章于 2024-06-17 17:35:06 发布

阅读量994

点赞数 20

分类专栏：生物计算工具文章标签： python

本文链接：https://blog.csdn.net/qq_40943760/article/details/139011204

版权

生物计算工具专栏收录该内容

49 篇文章 62 订阅

订阅专栏

对于经验不够丰富的人员，手动注释是不可行的。因此现在自动注释变得更容易让人接受。自动注释分为统计模型和深度学习模型。

参考：
[1] https://github.com/Starlitnightly/single_cell_tutorial
[2] https://github.com/theislab/single-cell-best-practices

基于统计模型的SCSA

SCSA是早期的一个细胞类型自动注释工具，可以根据簇特异性marker，查询数据库中每一种细胞类型的marker，选出吻合度最高的一类细胞，如果差异性不够大，被认为是Unknown细胞。

首先，加载数据：

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

adata = ov.read('./data/s4d8_dimensionality_reduction.h5ad')

聚类，这里把resolution设置为2，用于识别更多细胞小类。聚类结果如下：
fig1
SCSA有两个数据库可以选择：CellMarker和panglaodb。在omicverse中，SCSA具有以下参数：

foldchange：每个簇相对于别的簇的差异倍数，一般设置为1.5，值越大，使用的marker越少
pvalue：与foldchange对应，在计算差异倍数的时候会进行统计学差异显著性分析，一般为0.01
celltype：该参数包括了normal和cancer两个参数，当设置为cancer，会注释出来自cancersea数据库的12种肿瘤细胞亚型
target：要使用的数据库，目前SCSA支持CellMarker和panglaodb，其实也包括cancersea
tissue：可以使用scsa.get_model_tissue()列出所有支持的组织，默认是使用全部tissue
model_path：设置数据库的本地存放路径

使用CellMarker数据库

初始化方式如下：

scsa=ov.single.pySCSA(adata=adata,
                      foldchange=1.5,
                      pvalue=0.01,
                      celltype='normal',
                      target='cellmarker',
                      tissue='All',
                      model_path='./database/pySCSA_2023_v2_plus.db'                    
)

正式的注释环节涉及三个参数：

clustertype：注释依据的簇名，存放在adata.obs
cluster：需要注释的簇，默认是all
rank_rep: 是否需要重新计算差异表达基因，如果在注释之前已经运算了sc.tl.rank_genes_groups，我们可以设置成False

我们使用SCSA注释细胞，并打印簇的类型：

adata.uns['log1p']['base']=10
res=scsa.cell_anno(clustertype='leiden_res1',
               cluster='all',rank_rep=True)

# 自动打印注释结果
scsa.cell_anno_print()

# 注释结果添加到adata
scsa.cell_auto_anno(adata,clustertype='leiden_res1',
                    key='scsa_celltype_cellmarker')

Cluster:0 Cell_type:Natural killer cell|T cell Z-score:13.328|9.099
Nice:Cluster:1 Cell_type:B cell Z-score:17.324
Cluster:2 Cell_type:Natural killer cell|T cell Z-score:12.711|7.851
Nice:Cluster:3 Cell_type:B cell Z-score:12.564
Nice:Cluster:4 Cell_type:Natural killer cell Z-score:11.75
Nice:Cluster:5 Cell_type:B cell Z-score:17.695
Nice:Cluster:6 Cell_type:Natural killer cell Z-score:14.859
Nice:Cluster:7 Cell_type:T cell Z-score:12.816
Cluster:8 Cell_type:T cell|Naive CD4+ T cell Z-score:5.584|3.975
Cluster:9 Cell_type:T cell|Natural killer cell Z-score:8.559|5.724
Nice:Cluster:10 Cell_type:Monocyte Z-score:14.734
Nice:Cluster:11 Cell_type:B cell Z-score:15.721
Cluster:12 Cell_type:Monocyte|Natural killer T (NKT) cell Z-score:14.965|12.388
Cluster:13 Cell_type:T cell|Natural killer cell Z-score:10.451|8.418
Nice:Cluster:14 Cell_type:Natural killer T (NKT) cell Z-score:16.389
Cluster:15 Cell_type:Natural killer T (NKT) cell|T cell Z-score:13.173|11.666
Cluster:16 Cell_type:Natural killer cell|T cell Z-score:7.567|7.555
Cluster:17 Cell_type:Natural killer T (NKT) cell|B cell Z-score:16.351|14.131
Nice:Cluster:18 Cell_type:Red blood cell (erythrocyte) Z-score:12.859
Nice:Cluster:19 Cell_type:Natural killer cell Z-score:9.203
Nice:Cluster:20 Cell_type:B cell Z-score:16.841
Nice:Cluster:21 Cell_type:B cell Z-score:15.203
Cluster:22 Cell_type:Monocyte|Microglial cell Z-score:12.065|7.22
Nice:Cluster:23 Cell_type:B cell Z-score:12.986
Cluster:24 Cell_type:Hematopoietic stem cell|Natural killer T (NKT) cell Z-score:10.647|10.385
Nice:Cluster:25 Cell_type:Natural killer cell Z-score:14.398
Nice:Cluster:26 Cell_type:Red blood cell (erythrocyte) Z-score:2.196
Nice:Cluster:27 Cell_type:B cell Z-score:14.338
Cluster:28 Cell_type:T cell|Naive CD8+ T cell Z-score:8.457|6.469
Nice:Cluster:29 Cell_type:B cell Z-score:15.06
Cluster:30 Cell_type:Monocyte|Macrophage Z-score:12.186|8.279
Nice:Cluster:31 Cell_type:Natural killer cell Z-score:11.614

我们可视化注释结果：

# 可视化
ov.utils.embedding(adata,
                basis='X_mde',
                color=[ "leiden_res1","scsa_celltype_cellmarker",],
                title=['clustering', 'SCSA cell type'],
                palette=ov.palette()[15:],
                show=False,frameon='small',wspace=0.6)

fig2

使用panglaodb数据库

流程与使用CellMarker的一样：

scsa=ov.single.pySCSA(adata=adata,
                          foldchange=1.5,
                          pvalue=0.01,
                          celltype='normal',
                          target='panglaodb',
                          tissue='All',
                          model_path='./database/pySCSA_2023_v2_plus.db'
                          
)
res=scsa.cell_anno(clustertype='leiden_res1',
               cluster='all',rank_rep=True)
scsa.cell_auto_anno(adata,clustertype='leiden_res1',
                    key='scsa_celltype_panglaodb')
ov.utils.embedding(adata,
                basis='X_mde',
                color=[ "scsa_celltype_panglaodb","scsa_celltype_cellmarker",],
                #title=['Cell type'],
                palette=ov.palette()[15:],
                show=False,frameon='small',wspace=0.6)

两个数据库的注释结果对比如下所示：
fig3
我们发现，panglaodb注释出了Reticulocytes，Erythroid-like两种红细胞，T cells，T memory cells两种T细胞亚群。虽然看起来更准确了，但实际上我们在大类注释时不需要将细胞类型精细化。

omicverse中还提供了get_celltype_marker用于获取指定细胞类型的marker基因，其中会根据分类结果重新计算差异：

marker_dict=ov.single.get_celltype_marker(adata,clustertype='scsa_celltype_cellmarker')

我们在进行注释之前，可以先查看数据库中支持的tissue，目前支持的物种为Human和Mouse，以Human为例，支持的tissue可以通过以下方式查看：

scsa.get_model_tissue(species='Human')

基于深度学习的CellTypist

CellTypist来自Cross-tissue immune cell analysis reveals tissue-specific features in humans, Science, 2022，这算是第一篇发表在正刊的单细胞注释方法。CellTypist可以自定义高精度和低精度，对于不认识的细胞类型，会标记为新细胞。

首先加载数据：

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

import celltypist
from celltypist import models

adata = ov.read('./data/s4d8_dimensionality_reduction.h5ad')

CellTypist要求使用的归一化值为10,000，之前预处理的归一化值为500,000，所以我们提取原始counts，重新预处理归一化：

adata=adata.raw.to_adata()
ov.utils.retrieve_layers(adata,layers='counts')
# print(adata.X)

adata=ov.pp.preprocess(adata,mode='shiftlog|pearson',n_HVGs=2000,
                target_sum=1e4)

下载官方训练好的模型，以及查看模型保存的目录：

# Enabling `force_update = True` will overwrite existing (old) models.
models.download_models(force_update = True)
# 查看模型保存的目录
print(models.models_path)

CellTypist训练了很多模型，分别用于不同的精度和tissue，可以通过以下方式查看每个模型的描述信息：

models.models_description().head()

选择要使用的模型，比如，免疫细胞亚型的所有组织组合模型：

model = models.Model.load(model = 'Immune_All_Low.pkl')
print(model)

使用celltypist.annotate将细胞类型标签从模型知识传输到query数据集：

# Not run; predict cell identities using this loaded model.
#predictions = celltypist.annotate(adata_500, model = model, majority_voting = True, mode = 'best match')
# Alternatively, just specify the model name (recommended as this ensures the model is intact every time it is loaded).
# `p_thres` defaults to 0.5.
predictions = celltypist.annotate(adata, model = 'Immune_All_Low.pkl', majority_voting = True, mode = 'prob match', p_thres = 0.5)

结果包括预测的细胞类型标签（predicted_labels）、过度聚类结果（over_clustering）和局部子聚类中多数投票后的预测标签（majority_voting）：

adata_new = predictions.to_adata()
print(adata_new)

根据PCA表示进行UMAP可视化，然后查看这些标签的分布情况，由于部分簇是很多细胞预测类型都并行排列，可视化很长，我们只保留第一个细胞类型字符串：

adata_new.obs["simple_predict"] = adata_new.obs["majority_voting"].str.split('|').str[0]
print(adata_new.obs["simple_predict"])

sc.pp.neighbors(adata_new, n_neighbors=15,
                n_pcs=30, use_rep='scaled|original|X_pca')
sc.tl.leiden(adata_new, key_added="leiden_res1", resolution=1.0)

adata_new.obsm["X_mde"] = ov.utils.mde(adata.obsm["scaled|original|X_pca"])

ov.utils.embedding(adata_new,
                basis='X_mde',
                color=[ 
                       # "predicted_labels", 
                       # "over_clustering", 
                       # "majority_voting",
                       "simple_predict"
                       ],
                #title=['Cell type'],
                palette=sc.pl.palettes.default_102,
                show=False,frameon='small',wspace=0.45)

fig4
大部分细胞被标记为Unassigned，这是由于我们设定了预测概率0.5作为阈值，低于该阈值的细胞不被注释，表明CellTypist在注释未知的数据上效果有待提升。

Pan-cancer细胞类型注释-MetaTiME

一般的深度学习模型都是在健康样本上训练，对癌症或者疾病数据集缺乏知识。MetaTiME 学习了数百个肿瘤 scRNA-seq 数据中的数百万个单细胞。来自：MetaTiME integrates single-cell gene expression to characterize the meta-components of the tumor immune microenvironment，Nature communications，2023。

MetaTiME要求预处理过的数据，首先加载数据，这里已经是预处理过的数据：

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

adata = ov.read('./data/BCC_GSE123813_aPD1_with_metainfo.h5ad')
print(adata)
print(adata.X)

然后要去除恶性肿瘤细胞（用inferCNV识别），以获得更好的肿瘤微环境注释表现，在案例中，“isTME”保存了肿瘤微环境细胞。使用UMAP可视化肿瘤微环境细胞和身份：

ov.utils.embedding(adata,
                basis='X_umap',
                color=[ "isTME", "assign_ident"],
                # title=['Cell type'],
                #palette=sc.pl.palettes.default_102,
                show=False,frameon='small', wspace=0.45)

fig5
我们在这里就不去除恶性肿瘤了。然后，加载训练好的MetaTiME：

TiME_object=ov.single.MetaTiME(adata,mode='table')

我们可以对细胞进行过度聚类，这对细粒度的细胞注释非常有用。随着分辨率的提高，聚类的数量也会增加，为了尽最大可能地增加注释精度，我们将resolution设置为8：

TiME_object.overcluster(resolution=8,clustercol = 'overcluster',)

使用 TiME_object.predictTiME()来预测 TME 中的潜在细胞类型：

细胞亚型将存储在 adata.obs['MetaTiME'] 中
细胞类型将存储在 adata.obs['Major_MetaTiME'] 中

TiME_object.predictTiME(save_obs_name='MetaTiME')
print(adata.obs['MetaTiME'])
print(adata.obs['Major_MetaTiME'])

用UMAP可视化：

ov.utils.embedding(adata,
                  basis='X_umap',
                  color=['Major_MetaTiME','MetaTiME'],
                   show=False, add_outline=False, 
                   frameon='small',legend_fontoutline=2, wspace=0.6)