单细胞/时空文章代码复现——细胞类型注释及UMAP图优化

在上一篇推文中我们对所有分组样本进行了预处理及组内整合，并进行了聚类和差异基因筛选，本篇我们将着重介绍如何对单细胞转录组数据进行细胞类型注释以及UMAP图的优化技巧。

1. 单细胞转录组数据自动注释

在对人或小鼠进行单细胞细胞类型注释时，我们可以借用已发表的软件工具来基于相应的数据库进行自动注释。文章的方法部分尽管没有提及此步骤，但在作者上传的脚本中，我们可以看到作者使用SingleR软件尝试了自动注释。
SingleR是一个十分常用的细胞自动注释软件，它内置了7个人/小鼠的参考数据集，但是这7个数据集的下载非常不方便，建议大家从网络上搜索相应的资源保存至本地方便使用。我是从“生信菜鸟团”的网站中找到了一个长期的网盘下载链接。
言归正传，我们使用SingleR对损伤愈合组的数据进行自动注释。

library(Seurat)
library(SingleR)

###导入上一篇整合的多样本Seurat对象
obj <- readRDS('sc_inte.rds')
###导入人的参考数据集
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为参考数据集的变量名，label.main储存了细胞类型名称，我们指定针对seurat_clusters进行注释，也可以对每个细胞进行注释。

这样我们就通过SingleR获得了对每个cluster注释的结果。我们将该结果保存为表格，并使用UMAP图进行可视化。

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注释结果的UMAP图

自动注释软件的注释结果是非常粗糙的，比如，在SingleR注释结果中出现了软骨细胞（Chondrocytes）和神经元细胞（Neurons），这显然不符合常理。因此，如果需要获得一个精确的细胞注释结果，还是需要加入人工的校正或调整。

2. 单细胞转录组数据手动注释

依据文章正文部分，首先依据已有的细胞类型marker的表达来将cluster注释为细胞大类。因此，我们首先画出给出的marker在cluster间表达的小提琴图，并根据文中给出的标准进行大类注释。

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的小提琴图

以角质形成细胞（keratinocytes）为例，文中将KTR5或KRT10高表达的细胞定义为角质形成细胞，从小提琴图中可以看出，cluster1、2、4、6、7、11、14、15、17、20、21、22都可以被定义为角质形成细胞。我们将所有的marker表达的小提琴图都画出后，发现除了cluster19之外的所有cluster都可以根据已知marker的表达来定义细胞大类，但是我们并没能成功定义施旺细胞。
之后，文中对每个大类进行了亚类细分，将细胞大类进一步定义为更细致的细胞类型。根据上述31个cluster找到的差异基因，依照文章附表给出的信息将除cluster19之外的cluster都定义为更精细的细胞类型，这里我们依然未成功定义施旺细胞，以及FB-II、Bas-II类型。
我们将定义好的类型形成一个列表，并加入SingleR的结果作为对比。

我们可以看到，对于大部分cluster的细胞大类来说，SingleR的注释是准确的，只不过如果我们需要对亚类进行精确划分时，我们还是需要借助cluster间的差异表达基因。因此，我们可以在对自己的数据进行细胞注释时使用SingleR等自动注释软件进行细胞大类的划分，然后再根据marker基因进行更精细的注释。
然后将上述列表中的细胞类型注释结果添加至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图优化

在完成细胞类型注释后，我们来画UMAP图并进行修饰。其中，UMAP图中cluster的轮廓是由ggunchull包实现的，每个细胞大类的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)

###为每个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))
###计算每个label所处的位置 (细胞类型标签的定位，可以在生成后根据实际图像进行调整)
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为曲线外扩的距离，lty指定了线条类型为虚线。
因为我们是根据细胞亚类进行着色，所以细胞大类的label是使用geom_text_repel添加的，geom_text_repel中data即为我们希望label添加到的位置。
箭头状的坐标轴是通过scale_x_continuous(breaks = NULL)，scale_y_continuous(breaks = NULL)和theme中的axis.line实现的，只是不知道为何我的坐标轴的长短没有变化，只能后期通过AI修饰了。
经过上述的优化，来看下最终我们的注释结果的UMAP图。

复现图

文章原图

可以看到我们复现出的细胞类型注释结果和UMAP降维的大体分布与原文基本是一致的。
以上就是我们复现出的一个基本的注释范式。接下来，我们将对文中提到的所有单细胞数据进行整合分析。

往期内容：
单细胞/时空文章代码复现——数据预处理及整合