在上一篇推文中我們對(duì)所有分組樣本進(jìn)行了預(yù)處理及組內(nèi)整合，并進(jìn)行了聚類和差異基因篩選，本篇我們將著重介紹如何對(duì)單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行細(xì)胞類型注釋以及UMAP圖的優(yōu)化技巧。

1. 單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)自動(dòng)注釋

在對(duì)人或小鼠進(jìn)行單細(xì)胞細(xì)胞類型注釋時(shí)，我們可以借用已發(fā)表的軟件工具來(lái)基于相應(yīng)的數(shù)據(jù)庫(kù)進(jìn)行自動(dòng)注釋。文章的方法部分盡管沒(méi)有提及此步驟，但在作者上傳的腳本中，我們可以看到作者使用SingleR軟件嘗試了自動(dòng)注釋。
SingleR是一個(gè)十分常用的細(xì)胞自動(dòng)注釋軟件，它內(nèi)置了7個(gè)人/小鼠的參考數(shù)據(jù)集，但是這7個(gè)數(shù)據(jù)集的下載非常不方便，建議大家從網(wǎng)絡(luò)上搜索相應(yīng)的資源保存至本地方便使用。我是從“生信菜鳥(niǎo)團(tuán)”的網(wǎng)站中找到了一個(gè)長(zhǎng)期的網(wǎng)盤(pán)下載鏈接。
言歸正傳，我們使用SingleR對(duì)損傷愈合組的數(shù)據(jù)進(jìn)行自動(dòng)注釋。

library(Seurat)
library(SingleR)

###導(dǎo)入上一篇整合的多樣本Seurat對(duì)象
obj <- readRDS('sc_inte.rds')
###導(dǎo)入人的參考數(shù)據(jù)集
load('ref_Human_all.RData')

###SingleR注釋
obj_sr <- GetAssayData(obj, slot="data")
meta_cell <- SingleR(test = obj_sr, ref = ref_Human_all, labels = ref_Human_all$label.main, clusters=obj$seurat_clusters)
### ref_Human_all為參考數(shù)據(jù)集的變量名，label.main儲(chǔ)存了細(xì)胞類型名稱，我們指定針對(duì)seurat_clusters進(jìn)行注釋，也可以對(duì)每個(gè)細(xì)胞進(jìn)行注釋。

這樣我們就通過(guò)SingleR獲得了對(duì)每個(gè)cluster注釋的結(jié)果。我們將該結(jié)果保存為表格，并使用UMAP圖進(jìn)行可視化。

cell_type <- DataFrame(seurat_clusters=rownames(meta_cell), SingleR_label=meta_cell$labels)
write.csv(cell_type, 'SingleR_celltype.csv', quote=F)

cell_name <- cell_type$SingleR_label
names(cell_name) <- levels(obj)
obj <- RenameIdents(obj, cell_name)
obj$singleR_label <- Idents(obj)

DimPlot(obj, group.by='singleR_label', label=T, pt.size=0.2, label.size=3)

SingleR注釋結(jié)果的UMAP圖

自動(dòng)注釋軟件的注釋結(jié)果是非常粗糙的，比如，在SingleR注釋結(jié)果中出現(xiàn)了軟骨細(xì)胞（Chondrocytes）和神經(jīng)元細(xì)胞（Neurons），這顯然不符合常理。因此，如果需要獲得一個(gè)精確的細(xì)胞注釋結(jié)果，還是需要加入人工的校正或調(diào)整。

2. 單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)手動(dòng)注釋

依據(jù)文章正文部分，首先依據(jù)已有的細(xì)胞類型marker的表達(dá)來(lái)將cluster注釋為細(xì)胞大類。因此，我們首先畫(huà)出給出的marker在cluster間表達(dá)的小提琴圖，并根據(jù)文中給出的標(biāo)準(zhǔn)進(jìn)行大類注釋。

known_marker <- c('KRT5', 'KRT10', 
                  'COL1A1', 
                  'LYZ', 'HLA-DRA', 
                  'CD3D', 'NKG7', 
                  'PECAM1', 
                  'TPSAB1', 'TPSB2', 
                  'ACTA2', 'MYH11', 
                  'TYRP1', 'PMEL', 
                  'SOX10', 'SOX2')
VlnPlot(obj, features = c("KRT5", "KRT10"), ncol=1)

Marker gene的小提琴圖

以角質(zhì)形成細(xì)胞（keratinocytes）為例，文中將KTR5或KRT10高表達(dá)的細(xì)胞定義為角質(zhì)形成細(xì)胞，從小提琴圖中可以看出，cluster1、2、4、6、7、11、14、15、17、20、21、22都可以被定義為角質(zhì)形成細(xì)胞。我們將所有的marker表達(dá)的小提琴圖都畫(huà)出后，發(fā)現(xiàn)除了cluster19之外的所有cluster都可以根據(jù)已知marker的表達(dá)來(lái)定義細(xì)胞大類，但是我們并沒(méi)能成功定義施旺細(xì)胞。
之后，文中對(duì)每個(gè)大類進(jìn)行了亞類細(xì)分，將細(xì)胞大類進(jìn)一步定義為更細(xì)致的細(xì)胞類型。根據(jù)上述31個(gè)cluster找到的差異基因，依照文章附表給出的信息將除cluster19之外的cluster都定義為更精細(xì)的細(xì)胞類型，這里我們依然未成功定義施旺細(xì)胞，以及FB-II、Bas-II類型。
我們將定義好的類型形成一個(gè)列表，并加入SingleR的結(jié)果作為對(duì)比。

我們可以看到，對(duì)于大部分cluster的細(xì)胞大類來(lái)說(shuō)，SingleR的注釋是準(zhǔn)確的，只不過(guò)如果我們需要對(duì)亞類進(jìn)行精確劃分時(shí)，我們還是需要借助cluster間的差異表達(dá)基因。因此，我們可以在對(duì)自己的數(shù)據(jù)進(jìn)行細(xì)胞注釋時(shí)使用SingleR等自動(dòng)注釋軟件進(jìn)行細(xì)胞大類的劃分，然后再根據(jù)marker基因進(jìn)行更精細(xì)的注釋。
然后將上述列表中的細(xì)胞類型注釋結(jié)果添加至SeuratObject中。

obj <- subset(obj, seurat_clusters !='19')
obj$seurat_clusters <- droplevels(obj$seurat_clusters)

celist <- read.table('cell_type.list', header=T, sep='\t')
cell_name <- celist$cell_type
names(cell_name) <- levels(obj)
obj <- RenameIdents(obj, cell_name)
obj$cell_type <- Idents(obj)

Idents(obj) <- obj$seurat_clusters
cell_name <- celist$main_type
names(cell_name) <- levels(obj)
obj <- RenameIdents(obj, cell_name)
obj$main_type <- Idents(obj)

saveRDS(obj, 'cell_type.rds')

3. UMAP圖優(yōu)化

在完成細(xì)胞類型注釋后，我們來(lái)畫(huà)UMAP圖并進(jìn)行修飾。其中，UMAP圖中cluster的輪廓是由ggunchull包實(shí)現(xiàn)的，每個(gè)細(xì)胞大類的label是由ggrepel包添加的。

library(ggunchull)
library(ggrepel)
library(ggplot2)
library(dplyr)
###定義色板
cols <- c(colorRampPalette(c("#1f78b4", "#a6cee3"))(3),
          colorRampPalette(c("#33a02c", "#b2df8a"))(3), "#df65b0",  
          colorRampPalette(c("#fdbf6f", "#ff7f00"))(3), 
          "#696969", "#d8bfd8", "#008b00", "#fb9a99", "#8b5a2b", "#da70d6",
          colorRampPalette(c("#66cdaa", "#5f9ea0"))(2), 
          colorRampPalette(c("#cab2d6", "#6a3d9a"))(6),
          "#1c6597", "#cc6805", "#c6abd4", "#569395", "#dc65d8", "#f99190", "#057605","#575757"
          )
###為cell type重新排序
obj$cell_type <- factor(obj$cell_type, 
                        levels=c("FB-I","FB-III","FB-prolif",
                                 "Spi-I","Spi-II","Spi-mig",
                                 "Gra",
                                 "Bas-I","Bas-mig","Bas-prolif",
                                 "MEL","HF","PC-vSMC","LE","VE","Mast-cell",
                                 "Th","NK-cell",
                                 "Mono-Mac","cDC1","cDC2","DC3","Plasma_Bcell","LC"), 
                        ordered=T)

meta <- cbind(obj@meta.data, obj@reductions$umap@cell.embeddings)

###為每個(gè)cell type指定顏色
names(cols) <- levels(factor(c(levels(meta$cell_type), levels(meta$main_type)), 
                      levels=c(levels(meta$cell_type), levels(meta$main_type)), 
                      ordered=T))
###計(jì)算每個(gè)label所處的位置 (細(xì)胞類型標(biāo)簽的定位，可以在生成后根據(jù)實(shí)際圖像進(jìn)行調(diào)整)
main_type_med <- meta %>% 
                 group_by(main_type) %>% 
                 summarise(UMAP_1 = median(UMAP_1)-1, UMAP_2 = median(UMAP_2)-1)

ggplot(meta, aes(x = UMAP_1, y = UMAP_2)) +
       stat_unchull(aes(fill = main_type, color = main_type), 
                    alpha = 0.05, size = 1, lty = 2, delta=0.25, show.legend=F) +
       geom_point(aes(color = cell_type), size = 0.2, show.legend = FALSE) +
       scale_x_continuous(breaks = NULL) + 
       scale_y_continuous(breaks = NULL) +
       scale_color_manual(values = cols) +
       geom_text_repel(aes(label=main_type, color=main_type), 
                       fontface="bold", data=main_type_med, show.legend=F, size=6) +
       theme(aspect.ratio = 1,
             panel.background = element_blank(),
             panel.grid = element_blank(),
             axis.line = element_line(arrow = arrow(type = "closed")),
             axis.title = element_text(hjust = 0.05, size=12))

stat_unchull給cluster加上外面的圈，delta為曲線外擴(kuò)的距離，lty指定了線條類型為虛線。
因?yàn)槲覀兪歉鶕?jù)細(xì)胞亞類進(jìn)行著色，所以細(xì)胞大類的label是使用geom_text_repel添加的，geom_text_repel中data即為我們希望label添加到的位置。
箭頭狀的坐標(biāo)軸是通過(guò)scale_x_continuous(breaks = NULL)，scale_y_continuous(breaks = NULL)和theme中的axis.line實(shí)現(xiàn)的，只是不知道為何我的坐標(biāo)軸的長(zhǎng)短沒(méi)有變化，只能后期通過(guò)AI修飾了。
經(jīng)過(guò)上述的優(yōu)化，來(lái)看下最終我們的注釋結(jié)果的UMAP圖。

復(fù)現(xiàn)圖

文章原圖

可以看到我們復(fù)現(xiàn)出的細(xì)胞類型注釋結(jié)果和UMAP降維的大體分布與原文基本是一致的。
以上就是我們復(fù)現(xiàn)出的一個(gè)基本的注釋范式。接下來(lái)，我們將對(duì)文中提到的所有單細(xì)胞數(shù)據(jù)進(jìn)行整合分析。

往期內(nèi)容：
單細(xì)胞/時(shí)空文章代碼復(fù)現(xiàn)——數(shù)據(jù)預(yù)處理及整合