【原】肺腺癌單細(xì)胞數(shù)據(jù)集GSE189357復(fù)現(xiàn)(二)：細(xì)胞注釋

健明 2024-10-18 發(fā)布于廣東

展開全文

前言

Hello小伙伴們大家好，我是生信技能樹的小學(xué)徒”我才不吃蛋黃“。連著三天手術(shù)日，拖更了一天。據(jù)說生信夜班俠在實(shí)驗(yàn)室被背刺（血漿和腫瘤組織的多組學(xué)分析揭示了三陰性乳腺癌抗PD-L1免疫治療的核心蛋白），見了見世面，深感同情。夜班俠是剛從臨床到實(shí)驗(yàn)室，我是剛從實(shí)驗(yàn)室來臨床。我比較幸運(yùn)，來到了一個(gè)和諧的科室和一個(gè)氣氛超級(jí)好的組，免去了很多人際上的是非。但是在臨床上，你還會(huì)面臨其他各種問題。各種辛酸，慢慢道來，加油吧，東北老鐵?。?！加油吧，各位實(shí)驗(yàn)室和臨床的同（niu）行（ma）?。?！

好了，私貨夾帶完畢，開始進(jìn)入正題。

今天是肺腺癌單細(xì)胞數(shù)據(jù)集GSE189357復(fù)現(xiàn)系列第二期。第一期我們走了Seurat V5標(biāo)準(zhǔn)流程，利用harmony整合去批次后，按標(biāo)準(zhǔn)流程進(jìn)行了降維聚類分群。本期，我們將在第一期基礎(chǔ)上選擇合適的分辨率，對(duì)細(xì)胞亞群進(jìn)行注釋。

1.背景介紹

細(xì)胞注釋是單細(xì)胞分析的魂,是一切分析的基石，我們后續(xù)的所有分析都是圍繞細(xì)胞亞群展開的。但是無論是對(duì)于新手還是老手來說，細(xì)胞注釋都充滿了艱難與挑戰(zhàn)。

為什么說細(xì)胞注釋難，個(gè)人認(rèn)為有以下幾點(diǎn)：

一是注釋的方法多：我們可以使用SingleR、Cellassign等工具，這些工具可以通過比較未知樣本與已知細(xì)胞類型的參考數(shù)據(jù)集來自動(dòng)注釋細(xì)胞類型。也可以進(jìn)行手動(dòng)注釋，基于Marker基因的注釋：通過識(shí)別特定細(xì)胞類型的標(biāo)記基因（Marker genes）來手動(dòng)注釋細(xì)胞類型；基于數(shù)據(jù)庫(kù)的注釋：利用如CellMarker、PanglaoDB、CancerSEA等數(shù)據(jù)庫(kù)，這些數(shù)據(jù)庫(kù)提供了不同細(xì)胞類型的標(biāo)記基因集合。

二是準(zhǔn)確度不好把控：?jiǎn)渭?xì)胞細(xì)胞注釋的準(zhǔn)確度受多種因素影響，包括數(shù)據(jù)質(zhì)量、分析方法、使用的標(biāo)記基因（marker genes）數(shù)據(jù)庫(kù)、以及研究人員的專業(yè)知識(shí)等。影響因素變量多，準(zhǔn)確度就難以把控。

三是無法命名的細(xì)胞的處理問題：在單細(xì)胞測(cè)序數(shù)據(jù)分析中，有時(shí)候會(huì)遇到無法直接命名的細(xì)胞群體，這可能是因?yàn)檫@些細(xì)胞群體是新的或未知的細(xì)胞類型，或者是由于數(shù)據(jù)的復(fù)雜性導(dǎo)致的。處理無法命名的細(xì)胞群體是一個(gè)復(fù)雜的過程，通常需要多方面的分析和驗(yàn)證。在某些情況下，這些細(xì)胞群體可能代表了新的細(xì)胞類型，其發(fā)現(xiàn)可能對(duì)科學(xué)領(lǐng)域有重要的貢獻(xiàn)，而新手（甚至是老手）往往可能會(huì)直接把這一部分細(xì)胞忽略或者過濾掉。

四是亞群注釋命名的若干問題：亞群注釋命名方式有很多，目前缺乏統(tǒng)一的規(guī)范。細(xì)胞具有多樣性，同時(shí)，由于我們對(duì)細(xì)胞類型的認(rèn)識(shí)仍相對(duì)處于相對(duì)初級(jí)階段，因此準(zhǔn)確的亞群注釋仍道阻且長(zhǎng)。

五是沒有參考數(shù)據(jù)庫(kù)如何注釋：例如罕見腫瘤、新測(cè)組織類型、新物種等等。

六是比例較低的細(xì)胞類型的注釋問題：質(zhì)控嚴(yán)格與否？

2.細(xì)胞注釋

常見的分群細(xì)胞及其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)

通常我們第一層次降維聚類分群，可以將細(xì)胞分為三大類，包括：

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

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

言歸正傳，開始今天的代碼運(yùn)行：

2.1 首先加載R資源和單細(xì)胞數(shù)據(jù)，并設(shè)置分辨率：

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 創(chuàng)建新的文件夾，設(shè)置需要可視化的Maker，繪制分群氣泡圖：

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

genes_to_check = c('EPCAM','KRT19','CLDN4','SCGB1A1',  #上皮
                   'PECAM1' , 'CLO1A2', 'VWF',  #基質(zhì)
                   'CDH5', 'PECAM1', 'VWF','CLDN5',  #內(nèi)皮
                   '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', #中性粒細(xì)胞
                   '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分布情況

####構(gòu)建左下角坐標(biāo)軸
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

##圖中的標(biāo)簽和方框都是可以自定義的，例如下面這副刪掉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 人工肉眼注釋亞群

#####細(xì)胞生物學(xué)命名
celltype=data.frame(ClusterID=0:22,
                    celltype= 0:22) 

# 這里強(qiáng)烈依賴于生物學(xué)背景，看dotplot的基因表達(dá)量情況來人工審查單細(xì)胞亞群名字
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('../')

結(jié)語(yǔ)

本期，我們選擇resolution=0.5對(duì)肺腺癌單細(xì)胞數(shù)據(jù)集GSE189357進(jìn)行了手工細(xì)胞分群注釋。細(xì)胞注釋的方式有很多，詳情請(qǐng)參照單細(xì)胞天地前期推文：單細(xì)胞測(cè)序最好的教程（六）：細(xì)胞類型注釋。下一期，我們將在本期基礎(chǔ)上對(duì)細(xì)胞類型及基因進(jìn)行可視化，使用多種可視化方法，繪制DimPlot，F(xiàn)eaturePlot，ggplot，DoHeatmap圖，歡迎大家持續(xù)追更，謝謝！