肺腺癌单细胞数据集GSE189357复现(二)：细胞注释

前言

Hello小伙伴们大家好，我是生信技能树的小学徒”我才不吃蛋黄“。连着三天手术日，拖更了一天。据说生信夜班侠在实验室被背刺（血浆和肿瘤组织的多组学分析揭示了三阴性乳腺癌抗PD-L1免疫治疗的核心蛋白），见了见世面，深感同情。夜班侠是刚从临床到实验室，我是刚从实验室来临床。我比较幸运，来到了一个和谐的科室和一个气氛超级好的组，免去了很多人际上的是非。但是在临床上，你还会面临其他各种问题。各种辛酸，慢慢道来，加油吧，东北老铁！！！加油吧，各位实验室和临床的同（niu）行（ma）！！！

好了，私货夹带完毕，开始进入正题。

今天是肺腺癌单细胞数据集GSE189357复现系列第二期。第一期我们走了Seurat V5标准流程，利用harmony整合去批次后，按标准流程进行了降维聚类分群。本期，我们将在第一期基础上选择合适的分辨率，对细胞亚群进行注释。

1.背景介绍

细胞注释是单细胞分析的魂,是一切分析的基石，我们后续的所有分析都是围绕细胞亚群展开的。但是无论是对于新手还是老手来说，细胞注释都充满了艰难与挑战。

为什么说细胞注释难，个人认为有以下几点：

一是注释的方法多：我们可以使用SingleR、Cellassign等工具，这些工具可以通过比较未知样本与已知细胞类型的参考数据集来自动注释细胞类型。也可以进行手动注释，基于Marker基因的注释：通过识别特定细胞类型的标记基因（Marker genes）来手动注释细胞类型；基于数据库的注释：利用如CellMarker、PanglaoDB、CancerSEA等数据库，这些数据库提供了不同细胞类型的标记基因集合。

二是准确度不好把控：单细胞细胞注释的准确度受多种因素影响，包括数据质量、分析方法、使用的标记基因（marker genes）数据库、以及研究人员的专业知识等。影响因素变量多，准确度就难以把控。

三是无法命名的细胞的处理问题：在单细胞测序数据分析中，有时候会遇到无法直接命名的细胞群体，这可能是因为这些细胞群体是新的或未知的细胞类型，或者是由于数据的复杂性导致的。处理无法命名的细胞群体是一个复杂的过程，通常需要多方面的分析和验证。在某些情况下，这些细胞群体可能代表了新的细胞类型，其发现可能对科学领域有重要的贡献，而新手（甚至是老手）往往可能会直接把这一部分细胞忽略或者过滤掉。

四是亚群注释命名的若干问题：亚群注释命名方式有很多，目前缺乏统一的规范。细胞具有多样性，同时，由于我们对细胞类型的认识仍相对处于相对初级阶段，因此准确的亚群注释仍道阻且长。

五是没有参考数据库如何注释：例如罕见肿瘤、新测组织类型、新物种等等。

六是比例较低的细胞类型的注释问题：质控严格与否？

2.细胞注释

常见的分群细胞及其Maker：

# T Cells (CD3D, CD3E, CD8A), 
# B cells (CD19, CD79A, MS4A1 [CD20]), 
# Plasma cells (IGHG1, MZB1, SDC1, CD79A), 
# macrophages (CD68, CD163),
# 'CCL3L1' ,   #M2
# 'FABP4',  #M1
# Monocytes  (CD14),
# NK Cells (FGFBP2, FCG3RA, CX3CR1),  
# Photoreceptor cells (RCVRN), 
# Fibroblasts (FGF7, MME), 
# Neutrophil ('G0S2', 'S100A9','S100A8','CXCL8')
# Endothelial cells (PECAM1, VWF). 
# epi or tumor (EPCAM, KRT19, PROM1, ALDH1A1, CD24).
# immune (CD45+,PTPRC), epithelial/cancer (EpCAM+,EPCAM), 
# stromal (CD10+,MME,fibo or CD31+,PECAM1,endo) 

通常我们第一层次降维聚类分群，可以将细胞分为三大类，包括：

• immune (CD45+,PTPRC),
• epithelial/cancer (EpCAM+,EPCAM),
• stromal (CD10+,MME,fibro or CD31+,PECAM1,endo)

原文中提供的都是常见的细胞群：上皮细胞（EPCAM、KRT19、CLDN4）、基质（PECAM1、CLO1A2、VWF）、增殖性（MKI67、STMN1、PCNA）、T（CD3D、CD3E、CD2）、B（CD79A，IGHG1，MS4A1），NK（KLRD1、GNLY、KLRF1）和髓系（CSF1R、CSF3R、CD68）细胞。

言归正传，开始今天的代码运行：

2.1 首先加载R资源和单细胞数据，并设置分辨率：

rm(list=ls())
source('scRNA_scripts/lib.R')
source('scRNA_scripts/mycolors.R')
sce.all.int = readRDS('2-harmony/sce.all_int.rds') 
sel.clust = "RNA_snn_res.0.5"
sce.all.int <- SetIdent(sce.all.int, value = sel.clust)
table(sce.all.int@active.ident) 
colnames(sce.all.int@meta.data) 

2.2 创建新的文件夹，设置需要可视化的Maker，绘制分群气泡图：

dir.create("./3-Celltype")
setwd("./3-Celltype")
scRNA=sce.all.int

genes_to_check = c('EPCAM','KRT19','CLDN4','SCGB1A1',  #上皮
                   'PECAM1' , 'CLO1A2', 'VWF',  #基质
                   'CDH5', 'PECAM1', 'VWF','CLDN5',  #内皮
                   'LUM' , 'FGF7', 'MME',  #成纤维
                   'CD3D', 'CD3E', 'CD8A', 'CD4','CD2', #T
                   'AIF1', 'C1QC','C1QB','LYZ',  #巨噬
                   'MKI67', 'STMN1', 'PCNA',  #增殖
                   'CPA3' ,'CST3', 'KIT', 'TPSAB1','TPSB2',#肥大
                   'GOS2', 'S100A9','S100A8','CXCL8', #中性粒细胞
                   'KLRD1', 'GNLY', 'KLRF1','AREG', 'XCL2','HSPA6', #NK
                   'MS4A1','CD19', 'CD79A','IGHG1','MZB1', 'SDC1',  #B
                   'IGHD',  #MALT B
                   'CSF1R', 'CSF3R', 'CD68') #髓系

library(stringr)  
genes_to_check=str_to_upper(genes_to_check)
genes_to_check

p = DotPlot(scRNA, features = unique(genes_to_check),
                assay='RNA'  )  + coord_flip()
p
#ggsave('check_last_markers.pdf',height = 11,width = 11)

2.3 看看选择的resolution的umap分布情况

####构建左下角坐标轴
source('../scRNA_scripts/Bottom_left_axis.R')
result <- left_axes(scRNA)
axes <- result$axes
label <- result$label
umap =DimPlot(scRNA, reduction = "umap",cols = my36colors,pt.size = 0.8,
                  group.by = "RNA_snn_res.0.5",label = T,label.box = T) +
  NoAxes() + 
  theme(aspect.ratio = 1) +
  geom_line(data = axes,
            aes(x = x,y = y,group = group),
            arrow = arrow(length = unit(0.1, "inches"),
                          ends="last", type="closed")) +
  geom_text(data = label,
            aes(x = x,y = y,angle = angle,label = lab),fontface = 'italic')+
  theme(plot.title = element_blank())
umap
ggsave('RNA_snn_res.0.5_umap.pdf',width = 9,height = 7)

2.4 看看样本分组的umap

##图中的标签和方框都是可以自定义的，例如下面这副删掉label和label.box
sample_umap =DimPlot(scRNA, reduction = "umap",cols = my36colors,pt.size = 0.8,
                     group.by = "sample") +
  NoAxes() + 
  theme(aspect.ratio = 1) +
  geom_line(data = axes,
            aes(x = x,y = y,group = group),
            arrow = arrow(length = unit(0.1, "inches"),
                          ends="last", type="closed")) +
  geom_text(data = label,
            aes(x = x,y = y,angle = angle,label = lab),fontface = 'italic')+
  theme(plot.title = element_blank())
sample_umap
ggsave('sample_umap.pdf',width = 9,height = 7)

umap+sample_umap

2.5 人工肉眼注释亚群

#####细胞生物学命名
celltype=data.frame(ClusterID=0:22,
                    celltype= 0:22) 

# 这里强烈依赖于生物学背景，看dotplot的基因表达量情况来人工审查单细胞亚群名字
celltype[celltype$ClusterID %in% c(3,6,9,13,14,15,16,20,22 ),2]='Myeloid'
celltype[celltype$ClusterID %in% c( 5,11,12,17,19 ),2]='Epithelial'
celltype[celltype$ClusterID %in% c(0,1,21),2]='T&NK'
celltype[celltype$ClusterID %in% c( 8 ),2]='Fibro'
celltype[celltype$ClusterID %in% c( 18 ),2]='Proliferative'
celltype[celltype$ClusterID %in% c( 10 ),2]='Plasma'
celltype[celltype$ClusterID %in% c( 2),2]='B'
celltype[celltype$ClusterID %in% c( 7 ),2]='Endothelial'
celltype[celltype$ClusterID %in% c( 4),2]='Mast'

table(scRNA@meta.data$RNA_snn_res.0.5)
table(celltype$celltype)

scRNA@meta.data$celltype = "NA"
for(i in 1:nrow(celltype)){
  scRNA@meta.data[which(scRNA@meta.data$RNA_snn_res.0.5 == celltype$ClusterID[i]),'celltype'] <- celltype$celltype[i]}
table(scRNA@meta.data$celltype)


th=theme(axis.text.x = element_text(angle = 45, 
                                    vjust = 0.5, hjust=0.5)) 
library(patchwork)
celltype_umap =DimPlot(scRNA, reduction = "umap",cols = my36colors,pt.size = 0.8,
                       group.by = "celltype",label = T) +
  NoAxes() + 
  theme(aspect.ratio = 1) +
  geom_line(data = axes,
            aes(x = x,y = y,group = group),
            arrow = arrow(length = unit(0.1, "inches"),
                          ends="last", type="closed")) +
  geom_text(data = label,
            aes(x = x,y = y,angle = angle,label = lab),fontface = 'italic')+
  theme(plot.title = element_blank())
celltype_umap
ggsave('umap_by_celltype.pdf',width = 9,height = 7)

saveRDS(scRNA, "sce_celltype.rds")

setwd('../')

结语

本期，我们选择resolution=0.5对肺腺癌单细胞数据集GSE189357进行了手工细胞分群注释。细胞注释的方式有很多，详情请参照单细胞天地前期推文：单细胞测序最好的教程（六）：细胞类型注释。下一期，我们将在本期基础上对细胞类型及基因进行可视化，使用多种可视化方法，绘制DimPlot，FeaturePlot，ggplot，DoHeatmap图，欢迎大家持续追更，谢谢！