3-3. 自动注释(一)¶
1. 背景¶
我们在3-2的教程中,详细介绍了如何根据marker对细胞进行手动注释的方法。但是随着时代的发展,越来越多的组织,器官被人们注释完。很多初学者对marker的不熟悉,使得其手动注释的准确率较低。这时候,我们就在想,是否可以基于现有的注释数据,构建一个全面的注释图谱,使得我们只需要输入单细胞数据,就能自动输出注释结果。
在这里,我将自动注释算法分为需要GPU的深度学习模型以及不需要GPU的统计学模型两种,我们在(一)中先介绍不需要GPU的统计学模型自动注释。我们在此介绍的算法是SCSA,SCSA可以根据簇特异性marker,查询数据库中每一种细胞类型的marker,选出吻合度最高的一类细胞,如果差异性不够大,那么则选出两类细胞,这时可以被认为是Unknown细胞。
import omicverse as ov
print(f'omicverse version: {ov.__version__}')
import scanpy as sc
print(f'scanpy version: {sc.__version__}')
ov.ov_plot_set()
omicverse version: 1.5.0 scanpy version: 1.7.2
2. 加载数据¶
在这里,我们使用手动注释好的所用到的人类骨髓数据进行注释,使用注释好的数据可以方便我们比较手动注释与自动注释的准确率。
adata = ov.read('s4d8_manual_annotation.h5ad')
adata
AnnData object with n_obs × n_vars = 14814 × 2000
obs: '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', 'leiden_res1', 'major_celltype', 'manual_celltype', 'minor_celltype'
var: 'gene_ids', 'feature_types', 'genome', 'mt', 'ribo', 'hb', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'n_cells', 'percent_cells', 'robust', 'mean', 'var', 'residual_variances', 'highly_variable_rank', 'highly_variable_features'
uns: 'hvg', 'layers_counts', 'leiden', 'leiden_res1_colors', 'log1p', 'major_celltype_colors', 'manual_celltype_colors', 'minor_celltype_colors', 'neighbors', 'scaled|original|cum_sum_eigenvalues', 'scaled|original|pca_var_ratios'
obsm: 'X_mde', 'scaled|original|X_pca'
varm: 'scaled|original|pca_loadings'
layers: 'counts', 'lognorm', 'scaled', 'soupX_counts'
obsp: 'connectivities', 'distances'
3. 聚类¶
与手动注释类似,我们使用SCSA模型也需要对单细胞测序数据先进行聚类,在这里,我们的resolution=2可以设置地更大一些。
sc.pp.neighbors(adata, n_neighbors=15,
n_pcs=30,use_rep='scaled|original|X_pca')
sc.tl.leiden(adata, key_added="leiden_res1", resolution=2.0)
computing neighbors
finished (0:00:00)
running Leiden clustering
finished (0:00:01)
我们观察不同类别在UMAP图上的分布情况
adata.obsm["X_mde"] = ov.utils.mde(adata.obsm["scaled|original|X_pca"])
ov.utils.embedding(adata,
basis='X_mde',
color=[ "leiden_res1"],
title=['Clusters'],
palette=ov.palette()[:],
show=False,frameon='small',)
WARNING: Length of palette colors is smaller than the number of categories (palette length: 28, categories length: 32. Some categories will have the same color.
<AxesSubplot: title={'center': 'Clusters'}, xlabel='X_mde1', ylabel='X_mde2'>
4. SCSA自动注释¶
4.1 Cellmarker注释¶
很多paper里,将SCSA作为一个benchmark,认为其注释准确率较低,但实际这跟SCSA所依赖的数据库有关,作者一直将CellMarker数据库停留在1.0版本,直到今年omicverse将SCSA整合,并更新数据库后,这个19年发表的算法也更新了CellMarker数据库。而在omicverse中,我还添加了panglaodb作为第二个参考数据库供读者选择。
我们在这里详细介绍函数ov.single.pySCSA的参数:
- foldchange: 这是我们每个簇相对于别的簇的差异倍数,一般设置1.5即可,设置地越高用到的marker越少
- pvalue: 与foldchange对应,我们在计算差异倍数的时候会进行统计学差异显著性分析,通常由于细胞数的原因,pvalue会出现膨胀,所以p值都非常小,我们设置一个常见阈值0.01即可
- celltype: 这个参数包括了
normal和cancer两个参数,当我们设置成cancer的时候,我们可以注释出来自cancersea数据库的12种肿瘤细胞亚型 - target: 我们要使用的目的数据库,目前pySCSA中存放了
cellmarker,cancersea,panglaodb三个数据库 - tissue: 我们可以使用
scsa.get_model_tissue()列出所有支持的组织,默认是使用全部tissue - model_path: 这里可以设定我们使用的数据库,默认为'',将自动从figshare下载omicverse专用数据库,当然你也可以手动下载
注意事项:
celltype='cancer',target='cancersea'需要同时设定
数据库:
pySCSA_2023_v2_plus.db
scsa=ov.single.pySCSA(adata=adata,
foldchange=1.5,
pvalue=0.01,
celltype='normal',
target='cellmarker',
tissue='All',
model_path='temp/pySCSA_2023_v2_plus.db'
)
正式开始注释的环节,我们这里有三个参数
- clustertype: 注释依据的簇名,存放在adata.obs
- cluster: 需要注释的cluster,可以设置为list: ['1','2']表示只注释簇1和簇2
- rank_rep: 是否需要重新计算差异表达基因,由于当细胞数量上去后,差异表达基因的计算会较慢,如果我们已经运算了
sc.tl.rank_genes_groups,我们可以设置成False停止重复计算
adata.uns['log1p']['base']=10
res=scsa.cell_anno(clustertype='leiden_res1',
cluster='all',rank_rep=True)
res.head()
我们可以使用cell_anno_print来打印自动注释的结果
scsa.cell_anno_print()
Cluster:0 Cell_type:Natural killer cell|T cell Z-score:13.273|9.396 Nice:Cluster:1 Cell_type:B cell Z-score:17.404 Cluster:2 Cell_type:Natural killer cell|T cell Z-score:11.593|6.514 Cluster:3 Cell_type:Natural killer cell|T cell Z-score:12.76|8.216 Nice:Cluster:4 Cell_type:B cell Z-score:12.598 Cluster:5 Cell_type:Natural killer cell|T cell Z-score:9.54|8.235 Nice:Cluster:6 Cell_type:B cell Z-score:17.663 Nice:Cluster:7 Cell_type:Natural killer cell Z-score:14.884 Cluster:8 Cell_type:T cell|CD4+ T cell Z-score:11.135|5.887 Nice:Cluster:9 Cell_type:Natural killer cell Z-score:11.839 Cluster:10 Cell_type:T cell|Natural killer cell Z-score:8.536|5.758 Nice:Cluster:11 Cell_type:Monocyte Z-score:14.702 Nice:Cluster:12 Cell_type:Natural killer T (NKT) cell Z-score:16.376 Cluster:13 Cell_type:Monocyte|Natural killer T (NKT) cell Z-score:14.982|12.323 Cluster:14 Cell_type:T cell|Naive CD4+ T cell Z-score:4.286|3.828 Nice:Cluster:15 Cell_type:B cell Z-score:15.455 Cluster:16 Cell_type:Natural killer T (NKT) cell|B cell Z-score:16.348|14.152 Cluster:17 Cell_type:Natural killer T (NKT) cell|T cell Z-score:12.837|11.628 Cluster:18 Cell_type:T cell|Natural killer cell Z-score:12.177|9.216 Nice:Cluster:19 Cell_type:Red blood cell (erythrocyte) Z-score:12.92 Nice:Cluster:20 Cell_type:B cell Z-score:15.318 Nice:Cluster:21 Cell_type:B cell Z-score:13.077 Cluster:22 Cell_type:Monocyte|Natural killer T (NKT) cell Z-score:12.582|6.915 Nice:Cluster:23 Cell_type:B cell Z-score:16.847 Cluster:24 Cell_type:Hematopoietic stem cell|Natural killer T (NKT) cell Z-score:10.647|10.385 Nice:Cluster:25 Cell_type:Red blood cell (erythrocyte) Z-score:2.196 Nice:Cluster:26 Cell_type:Natural killer cell Z-score:14.092 Nice:Cluster:27 Cell_type:B cell Z-score:14.312 Nice:Cluster:28 Cell_type:B cell Z-score:15.943 Cluster:29 Cell_type:T cell|Naive CD8+ T cell Z-score:7.808|6.713 Cluster:30 Cell_type:Monocyte|Macrophage Z-score:12.19|8.303 Nice:Cluster:31 Cell_type:Natural killer cell Z-score:11.614
我们使用cell_auto_anno将分数最高的细胞类型注释上,但这个函数不会区分非Nice的细胞注释结果,还存在一些误差。
scsa.cell_auto_anno(adata,clustertype='leiden_res1',
key='scsa_celltype_cellmarker')
...cell type added to scsa_celltype_cellmarker on obs of anndata
我们发现自动注释中多注释出了Hematopoietic stem cell,这类细胞可能是注释偏差,其对应的簇为24,我们发现其z-score没有高于第二类两倍,10.647|10.385。此外,我们自动注释没有注释出neuron,neuron对应的是胶质细胞,这是自动注释的局限性
ov.utils.embedding(adata,
basis='X_mde',
color=[ "major_celltype","scsa_celltype_cellmarker",],
title=['Cell type'],
palette=ov.palette()[15:],
show=False,frameon='small',wspace=0.45)
WARNING: The title list is shorter than the number of panels. Using 'color' value instead for some plots.
[<AxesSubplot: title={'center': 'Cell type'}, xlabel='X_mde1', ylabel='X_mde2'>,
<AxesSubplot: title={'center': 'scsa_celltype_cellmarker'}, xlabel='X_mde1', ylabel='X_mde2'>]
4.2 panglaodb注释¶
我们还可以尝试panglaodb数据库来注释
scsa=ov.single.pySCSA(adata=adata,
foldchange=1.5,
pvalue=0.01,
celltype='normal',
target='panglaodb',
tissue='All',
model_path='temp/pySCSA_2023_v2_plus.db'
)
res=scsa.cell_anno(clustertype='leiden_res1',
cluster='all',rank_rep=True)
res.head()
scsa.cell_auto_anno(adata,clustertype='leiden_res1',
key='scsa_celltype_panglaodb')
...cell type added to scsa_celltype_panglaodb on obs of anndata
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.45)
[<AxesSubplot: title={'center': 'scsa_celltype_panglaodb'}, xlabel='X_mde1', ylabel='X_mde2'>,
<AxesSubplot: title={'center': 'scsa_celltype_cellmarker'}, xlabel='X_mde1', ylabel='X_mde2'>]
有趣的是,panglaodb注释出了Reticulocytes,Erythroid-like两种红细胞,T cells,T memory cells两种T细胞亚群。虽然看起来更准确了,但实际上我们在大类注释时不需要将细胞类型精细化。我们还提供了一个函数get_celltype_marker获取指定细胞类型的marker基因
marker_dict=ov.single.get_celltype_marker(adata,clustertype='scsa_celltype_cellmarker')
marker_dict.keys()
...get cell type marker
ranking genes
finished (0:00:15)
dict_keys(['B cell', 'Hematopoietic stem cell', 'Monocyte', 'Natural killer T (NKT) cell', 'Natural killer cell', 'Red blood cell (erythrocyte)', 'T cell'])
marker_dict['B cell']
array(['AFF3', 'BACH2', 'CD74', 'EBF1', 'BANK1', 'RALGPS2', 'IGHM', 'LYN',
'FCHSD2', 'ZCCHC7', 'PAX5', 'MEF2C', 'UBE2E2'], dtype=object)
我们可以查询所有支持的tissue,目前支持两种物种:Human与Mouse。我们在后续的更新中会陆续加入斑马鱼,植物等数据库
scsa.get_model_tissue(species='Human')
Version V2.1 [2023/06/27] DB load: GO_items:47347,Human_GO:3,Mouse_GO:3, CellMarkers:82887,CancerSEA:1574,PanglaoDB:24223 Ensembl_HGNC:61541,Ensembl_Mouse:55414 ######################################################################################################################## ------------------------------------------------------------------------------------------------------------------------ Species:Human Num:298 ------------------------------------------------------------------------------------------------------------------------ 1: Abdomen 2: Abdominal adipose tissue 3: Abdominal fat pad 4: Acinus 5: Adipose tissue 6: Adrenal gland 7: Adventitia 8: Airway 9: Airway epithelium 10: Allocortex 11: Alveolus 12: Amniotic fluid 13: Amniotic membrane 14: Ampullary 15: Anogenital tract 16: Antecubital vein 17: Anterior cruciate ligament 18: Anterior presomitic mesoderm 19: Aorta 20: Aortic valve 21: Artery 22: Arthrosis 23: Articular Cartilage 24: Ascites 25: Ascitic fluid 26: Atrium 27: Basal airway 28: Basilar membrane 29: Beige Fat 30: Bile duct 31: Biliary tract 32: Bladder 33: Blood 34: Blood vessel 35: Bone 36: Bone marrow 37: Brain 38: Breast 39: Bronchial vessel 40: Bronchiole 41: Bronchoalveolar lavage 42: Bronchoalveolar system 43: Bronchus 44: Brown adipose tissue 45: Calvaria 46: Capillary 47: Cardiac atrium 48: Cardiovascular system 49: Carotid artery 50: Carotid plaque 51: Cartilage 52: Caudal cortex 53: Caudal forebrain 54: Caudal ganglionic eminence 55: Cavernosum 56: Central amygdala 57: Central nervous system 58: Cerebellum 59: Cerebral organoid 60: Cerebrospinal fluid 61: Cervix 62: Choriocapillaris 63: Chorionic villi 64: Chorionic villus 65: Choroid 66: Choroid plexus 67: Colon 68: Colon epithelium 69: Colorectum 70: Cornea 71: Corneal endothelium 72: Corneal epithelium 73: Coronary artery 74: Corpus callosum 75: Corpus luteum 76: Cortex 77: Cortical layer 78: Cortical thymus 79: Decidua 80: Deciduous tooth 81: Dental pulp 82: Dermis 83: Diencephalon 84: Distal airway 85: Dorsal forebrain 86: Dorsal root ganglion 87: Dorsolateral prefrontal cortex 88: Ductal tissue 89: Duodenum 90: Ectocervix 91: Ectoderm 92: Embryo 93: Embryoid body 94: Embryonic Kidney 95: Embryonic brain 96: Embryonic heart 97: Embryonic prefrontal cortex 98: Embryonic stem cell 99: Endocardium 100: Endocrine 101: Endoderm 102: Endometrium 103: Endometrium stroma 104: Entorhinal cortex 105: Epidermis 106: Epithelium 107: Esophageal 108: Esophagus 109: Eye 110: Fat pad 111: Fetal brain 112: Fetal gonad 113: Fetal heart 114: Fetal ileums 115: Fetal kidney 116: Fetal liver 117: Fetal lung 118: Fetal thymus 119: Fetal umbilical cord 120: Fetus 121: Foreskin 122: Frontal cortex 123: Fundic gland 124: Gall bladder 125: Gastric corpus 126: Gastric epithelium 127: Gastric gland 128: Gastrointestinal tract 129: Germ 130: Germinal center 131: Gingiva 132: Gonad 133: Gut 134: Hair follicle 135: Head 136: Head and neck 137: Heart 138: Heart muscle 139: Hippocampus 140: Ileum 141: Iliac crest 142: Inferior colliculus 143: Intervertebral disc 144: Intestinal crypt 145: Intestine 146: Intrahepatic cholangio 147: Jejunum 148: Kidney 149: Lacrimal gland 150: Large Intestine 151: Large intestine 152: Larynx 153: Lateral ganglionic eminence 154: Left lobe 155: Ligament 156: Limb bud 157: Limbal epithelium 158: Liver 159: Lumbar vertebra 160: Lung 161: Lymph 162: Lymph node 163: Lymphatic vessel 164: Lymphoid tissue 165: Malignant pleural effusion 166: Mammary epithelium 167: Mammary gland 168: Medial ganglionic eminence 169: Medullary thymus 170: Meniscus 171: Mesenchyme 172: Mesoblast 173: Mesoderm 174: Microvascular endothelium 175: Microvessel 176: Midbrain 177: Middle temporal gyrus 178: Milk 179: Molar 180: Muscle 181: Myenteric plexus 182: Myocardium 183: Myometrium 184: Nasal concha 185: Nasal epithelium 186: Nasal mucosa 187: Nasal polyp 188: Nasopharyngeal mucosa 189: Nasopharynx 190: Neck 191: Neocortex 192: Nerve 193: Nose 194: Nucleus pulposus 195: Olfactory neuroepithelium 196: Omentum 197: Optic nerve 198: Oral cavity 199: Oral mucosa 200: Osteoarthritic cartilage 201: Ovarian cortex 202: Ovarian follicle 203: Ovary 204: Oviduct 205: Palatine tonsil 206: Pancreas 207: Pancreatic acinar tissue 208: Pancreatic duct 209: Pancreatic islet 210: Periodontal ligament 211: Periodontium 212: Periosteum 213: Peripheral blood 214: Peritoneal fluid 215: Peritoneum 216: Pituitary 217: Pituitary gland 218: Placenta 219: Plasma 220: Pleura 221: Pluripotent stem cell 222: Polyp 223: Posterior fossa 224: Posterior presomitic mesoderm 225: Prefrontal cortex 226: Premolar 227: Presomitic mesoderm 228: Primitive streak 229: Prostate 230: Pulmonary arteriy 231: Pyloric gland 232: Rectum 233: Renal glomerulus 234: Respiratory tract 235: Retina 236: Retinal organoid 237: Retinal pigment epithelium 238: Right ventricle 239: Saliva 240: Salivary gland 241: Scalp 242: Sclerocorneal tissue 243: Seminal plasma 244: Septum transversum 245: Serum 246: Sinonasal mucosa 247: Sinus tissue 248: Skeletal muscle 249: Skin 250: Small intestine 251: Soft tissue 252: Sperm 253: Spinal cord 254: Spleen 255: Sputum 256: Stomach 257: Subcutaneous adipose tissue 258: Submandibular gland 259: Subpallium 260: Subplate 261: Subventricular zone 262: Superior frontal gyrus 263: Sympathetic ganglion 264: Synovial fluid 265: Synovium 266: Taste bud 267: Tendon 268: Testis 269: Thalamus 270: Thymus 271: Thyroid 272: Tongue 273: Tonsil 274: Tooth 275: Trachea 276: Transformed artery 277: Trophoblast 278: Umbilical cord 279: Umbilical cord blood 280: Umbilical vein 281: Undefined 282: Urine 283: Urothelium 284: Uterine cervix 285: Uterus 286: Vagina 287: Vein 288: Venous blood 289: Ventral thalamus 290: Ventricular and atrial 291: Ventricular zone 292: Vessel 293: Visceral adipose tissue 294: Vocal cord 295: Vocal fold 296: White adipose tissue 297: White matter ########################################################################################################################