3.基于python的scRNA-seq细胞类型注释-自动注释（从参考到查询）

从参考数据集到查询数据集的注释又被称为迁移注释，其实属于自动注释，只不过这里的注释更具有领域特性，因为我们想在指定的参考数据集上训练深度学习模型，然后迁移到查询数据集，这种AI注释方式适合用于自动化注释特定疾病。目前TOSICA是迁移注释的代表方法。

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

加载参考数据集

参考数据需要与查询数据保持一致的预处理方法。归一化一般受到sc.pp.normalize_total中target_sum参数的影响，在scanpy中，默认值是自适应的，这是导致注释出错的一个原因，因此，我们固定设置target_sum参数为1e4。

在下面的案例中，我们使用omicverse来预处理原始counts，target_sum参数设置为50*1e4。

高质量的参考数据集是必要的，这里介绍了Cell Blast的参考数据集，Cell Blast收集了100多个高质量的单细胞带注释数据集，下面我们使用Zheng数据集来注释前面的骨髓数据集。Zheng作为参考数据集：

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

ref_adata = ov.utils.read('./data/Zheng.h5ad')
print(ref_adata)

TOSICA要求参考数据集的细胞数量不能少于30,000个，否则模型训练效果不佳。如果参考数据集的规模为100,000，随机选择30,000个细胞作为训练集也足够保持TOSICA的性能。

因此，我们在Zheng中随机选取30,000个细胞，使用random.sample进行无放回采样：

import random
cell_idx=list(random.sample(ref_adata.obs.index.tolist(),30000))
ref_adata=ref_adata[cell_idx]
print(ref_adata)

注释的细胞类型在cell_ontology_class中查看：

ref_adata.obs['cell_ontology_class'].cat.categories

检查归一化程度：

ref_adata.X.max() # 173.0

这是一个较大的整数，表明是原始counts，因此，我们对它进行预处理：

ref_adata = ov.pp.preprocess(ref_adata, mode='shiftlog|pearson', n_HVGs=3000)

加载查询数据集

这里的查询数据集使用前面的骨髓数据集，为了使预处理过程与参考数据集一致，我们需要提取counts，然后进行预处理，由于之前预处理流程一样，所以这里不再重新预处理。直接加载查询数据集：

adata = ov.utils.read('data/s4d8_dimensionality_reduction.h5ad')
# adata=adata.raw.to_adata()
# ov.utils.retrieve_layers(adata,layers='counts')
# print('raw count adata:',adata.X.max())

# adata=ov.pp.preprocess(adata,mode='shiftlog|pearson',n_HVGs=3000)

我们发现，参考数据集和查询数据集分别提取了3000个和2000个HVG，但是它们是不一致的，我们需要提取重叠的高变基因：

# 获取参考数据集和查询数据集重叠的HVG
print(ref_adata)
ref_adata = ref_adata[:, ref_adata.var.highly_variable_features]
ref_adata.var_names_make_unique()
adata.var_names_make_unique()
ret_gene=list(set(adata.var_names) & set(ref_adata.var_names))
print(len(ret_gene))

# 根据共同HVG切片数据集
print(ret_gene)
print(adata)
adata=adata[:,ret_gene]
ref_adata=ref_adata[:,ret_gene]

print(f"The max of ref_adata is {ref_adata.X.max()}, query_data is {adata.X.max()}",)
# The max of ref_adata is 11.3140869140625, query_data is 10.989398002624512

此时看ref_adata与adata的最大值，发现其相接近，表明二者具有相似的归一化值。

TOSICA训练以及预测

TOSICA需要知识库，omicverse可以自动下载：

ov.utils.download_tosica_gmt()

然后初始化TOSICA，这里涉及到以下参数：

• adata：参考 adata 对象
• gmt_path：默认的预先准备的掩码，或者 .gmt 文件的路径。可以使用 ov.utils.download_tosica_gmt() 来获取基因集。
• depth：Transformer 模型的深度，设置为 2 时可能会发生内存泄漏。
• label_name：在 adata.obs 中表示细胞类型的参考键。
• project_path：TOSICA 模型的保存路径。
• batch_size：表示在单次传递中传递给程序进行训练的细胞数量。

初始化如下：

tosica_obj=ov.single.pyTOSICA(adata=ref_adata,
                              gmt_path='genesets/GO_bp.gmt', depth=1,
                              label_name='cell_ontology_class',
                              project_path='hGOBP_demo',
                              batch_size=8)

然后训练TOSICA：

# pre_weights： 预训练权重的路径
# lr： 学习率
# epochs：epochs 的个数
# lrf： 最后一层的学习率
tosica_obj.train(epochs=5)

tosica_obj.save()
tosica_obj.load()

在查询数据集上预测：

new_adata=tosica_obj.predicted(pre_adata=adata)

new_adata.obsm=adata[new_adata.obs.index].obsm.copy()
new_adata.obsp=adata[new_adata.obs.index].obsp.copy()
print(new_adata)
"""
AnnData object with n_obs × n_vars = 14814 × 299
    obs: 'Prediction', 'Probability', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_20_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'total_counts_hb', 'log1p_total_counts_hb', 'pct_counts_hb', 'outlier', 'mt_outlier', 'scDblFinder_score', 'scDblFinder_class', 'log2_nUMIs', 'log2_nGenes'
    obsm: 'X_tsne', 'X_tsne_counts', 'X_tsne_scaled', 'X_umap', 'X_umap_counts', 'X_umap_scaled', 'counts|original|X_pca', 'scaled|original|X_pca'
    obsp: 'connectivities', 'distances'
"""

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

TOSICA的可解释性

TOSICA根据注意力来解释预测，首先我们将预测细胞类型数量小于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)]

然后，用scanpy的差异计算，得到各细胞类型中注意力最高的差异通路：

sc.tl.rank_genes_groups(new_adata, 'Prediction', method='wilcoxon')
sc.pl.rank_genes_groups_dotplot(new_adata,
                                n_genes=3,standard_scale='var',)

如果想获取特定细胞的通路，可以使用 sc.get.rank_genes_groups_df，比如，想获取B细胞特定的通路：

degs = sc.get.rank_genes_groups_df(new_adata, group='B cell', key='rank_genes_groups',
                                            pval_cutoff=0.05)
degs.head()

可视化通路的表达情况：

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