相见恨晚！不同数据，不同注释方法，来对号入座！

细胞注释

大家好，我是风风。欢迎来到风风的从零开始单细胞系列。在之前的单细胞分析流程中，我们简单地根据细胞的标记基因，也就是markers，对聚类得到的不同亚群进行细胞注释，从而去识别不同的亚群所具有的生物学功能。今天我们就来一起聊一聊“细胞注释”这件事儿。

在单细胞数据分析中，最难解决、最具挑战、也是最具创新的部分可以说是对分析结果的解释。而其中最直接的结果就是获得细胞亚群。但是，在对这种结果进行解释的时候，存在一个较大的难题，就是要确定这些亚群中的每一个亚群所代表的生物学状态，这需要我们有极强的专业背景，比如说我们识别T细胞亚群，就需要对T细胞的背景知识有背景了解之后，根据背景知识来联系当前数据分析结果和既往生物背景之间的关系。

当然，由于我们经常阅读文献，往往在看到亚群分析的结果之后，脑海里就能够马上知道这大概是什么样的细胞类型。可是并不是所有做单细胞分析的人都能够达到这样的能力水平啊，而且即使有这样的能力，也没办法全知全能把所有细胞类型的特征都印刻在脑海里吧？当然，你可以说不是还有很多数据库嘛，比如想Cellmarkers数据库进行自动注释。这当然可以，但是数据库的注释比起使用代码来注释相对来说限制更大。

因此，我们可以使用各种计算方法，利用一些已知的背景数据或者一些综述收集的markers，对未定性的细胞亚群进行细胞注释。其实这个过程就是比对注释的过程，我们在常规bulk-seq分析中也经常用到这种思想，例如最常见的GO和KEGG富集分析，而相应地，GO和KEGG数据集的一些标记基因也可以成为我们进行细胞注释的依据。此外，我们还可以直接将我们的聚类结果对应的表达谱与已发表的参考数据集进行比较，从而得到相应的注释结果。

接下来，我们看看不同的数据集用不同的注释方法，可以怎么进行细胞注释：

利用R包数据注释

我们前面说了，细胞注释可以用单细胞表达谱与一些既往发表的参考数据集进行比较，然后可以根据最相似的参照结果对我们的聚类得到的每个细胞亚群分配细胞标签——markers，那什么样的结果才叫“最相似”呢？我们就需要借助一些机器学习的算法来判断。我们明确一下，是任何已经发表的、有markers的RNA-seq数据集都可以作为背景参考，不仅限于单细胞文献，也可以是RNA-seq或者综述总结的结果。

这种自动注释的方法优势是不需要分析人员有自己的背景知识（有时候分析人员可能是一些技术人员，不具有强的专业背景嘛），但是劣势也很明显，就是注释的结果可能不如根据专业背景和分析目的来进行注释那么准确。

方法的实现只需要使用特定的R包进行了，比如celldex包。celldex包是一个非常强大的注释R包，用来匹配的背景数据集非常多，甚至可以说，即使是人工注释，其结果也不一定能够比celldex包的注释结果更加准确。celldex包的数据大部分来源于RNA seq或芯片数据集。关于celldex包的具体信息以及相关数据集大家可以自行检索一下发表的文章。接下来我们将演示利用SingleR包和celldex包的注释过程。分析数据直接在后天回复关键词即可获取：

# 安装R包# BiocManager::install("SingleR")# BiocManager::install("celldex")# BiocManager::install("pheatmap")# BiocManager::install("scran")# BiocManager::install("GSEABase")# BiocManager::install("AUCell")# 加载R包library(SingleR)library(celldex)library(pheatmap)# 加载分析数据，这里使用前面用过的10x GEnomics的外周血数据进行演示load(file = "pbmc.Rda")sce.pbmc## 载入需要的程辑包：SingleCellExperiment## class: SingleCellExperiment ## dim: 33694 3985 ## metadata(1): Samples## assays(2): counts logcounts## rownames(33694): RP11-34P13.3 FAM138A ... AC213203.1 FAM231B## rowData names(2): ID Symbol## colnames(3985): AAACCTGAGAAGGCCT-1 AAACCTGAGACAGACC-1 ...## TTTGTCAGTTAAGACA-1 TTTGTCATCCCAAGAT-1## colData names(4): Sample Barcode sizeFactor label## reducedDimNames(2): PCA TSNE## mainExpName: NULL## altExpNames(0):

# 下载celldex的注释信息，直接将代码前面的#号去除即可下载，这里已经下载完毕，直接加载进来分析

# ref <- BlueprintEncodeData()# save(ref, file = "annotation.Rda")load(file = "annotation.Rda")ref## class: SummarizedExperiment ## dim: 19859 259 ## metadata(0):## assays(1): logcounts## rownames(19859): TSPAN6 TNMD ... LINC00550 GIMAP1-GIMAP5## rowData names(0):## colnames(259): mature.neutrophil## CD14.positive..CD16.negative.classical.monocyte ...## epithelial.cell.of.umbilical.artery.1## dermis.lymphatic.vessel.endothelial.cell.1## colData names(3): label.main label.fine label.ont

# 注释的代码就是下面这一句，使用SingleR()匹配数据集和注释信息，用注释信息对亚群进行注释

pred <- SingleR(test = sce.pbmc, ref = ref, labels = ref$label.main)table(pred$labels)## ## B-cells CD4+ T-cells CD8+ T-cells DC Eosinophils Erythrocytes ## 549 773 1274 1 1 5 ## HSC Monocytes NK cells ## 14 1117 251# 将结果进行可视化，这里绘制热图plotScoreHeatmap(pred)

# 我们发现，CD4+T细胞和CD8+T细胞的区分比较模糊，我们的聚类并不能有效区分CD4+T细胞和CD8+T细胞，这可能是由于T细胞亚群内存在其他的异质性因素（比如T细胞活化），这些因素对无监督方法的影响比预期的CD4+/CD8+的区别更大

tab <- table(Assigned = pred$pruned.labels, Cluster = colLabels(sce.pbmc))tab[1:5,1:5]## Cluster## Assigned 1 2 3 4 5## B-cells 0 0 529 0 0## CD4+ T-cells 0 232 0 1 1## CD8+ T-cells 5 274 8 0 340## DC 0 0 0 0 0## Eosinophils 0 0 0 0 0pheatmap(log2(tab + 1), color = colorRampPalette(c("white", "steelblue"))(100))

从基因集中获取markers进行标记

前面我们说过，除了R包自动注释，我们还可以采用的一种方法就是使用KEGG和GO这一类富集数据库，明确地确定在每个细胞中高表达的标记基因list，这种方法不需要将单个细胞与参考数据集的表达值相匹配，只需要是基因名进行比对，类似于富集分析的思想，这样做更快、更方便。

我们使用来自Zeisel等人在2015年发表的研究的神经元细胞类型标记物来实现这种方法

library(scran)library(GSEABase)library(AUCell)# 一样是首先加载数据load(file = "zeisel.Rda")sce.zeisel## class: SingleCellExperiment ## dim: 19839 2816 ## metadata(0):## assays(2): counts logcounts## rownames(19839): 0610005C13Rik 0610007N19Rik ... Zzef1 Zzz3## rowData names(2): featureType Ensembl## colnames(2816): 1772071015_C02 1772071017_G12 ... 1772063068_D01## 1772066098_A12## colData names(10): tissue group # ... level1class level2class## reducedDimNames(0):## mainExpName: NULL## altExpNames(2): ERCC repeat# 接下来我们先对数据进行wilcox.test差异分析wilcox.z <- pairwiseWilcox(sce.zeisel, sce.zeisel$level1class, lfc=1, # 设置logFC值筛选差异表达结果 direction="up") # 这里只保留上调的结果markers.z <- getTopMarkers(wilcox.z$statistics, wilcox.z$pairs, pairwise=FALSE, n=50)lengths(markers.z)## astrocytes_ependymal endothelial-mural interneurons ## 79 83 118 ## microglia oligodendrocytes pyramidal CA1 ## 69 81 125 ## pyramidal SS ## 149# 接下来我们将一个大脑的scRNA seq数据集作为测试数据进行比对注释library(scRNAseq)load(file = "tasic.Rda")sce.tasic## class: SingleCellExperiment ## dim: 24058 1809 ## metadata(0):## assays(1): counts## rownames(24058): 0610005C13Rik 0610007C21Rik ... mt_X57780 tdTomato## rowData names(0):## colnames(1809): Calb2_tdTpositive_cell_1 Calb2_tdTpositive_cell_2 ...## Rbp4_CTX_250ng_2 Trib2_CTX_250ng_1## colData names(13): sample_title mouse_line ... secondary_type## aibs_vignette_id## reducedDimNames(0):## mainExpName: NULL## altExpNames(1): ERCC

# 接下来，我们使用AUCell软件包来识别在每个细胞中高度表达的标记基因集，这个方法是根据基因在每个细胞内的表达值进行排序，并构建一个反应曲线，即每个标记集中出现的基因数量随着排序的增加而增加。然后，使用ROC分析计算每个标志物集的曲线下面积（AUC），量化这些标志物在该细胞中最高度表达的基因中的富集程度，这种方式大致类似于在标记集内和标记集外的基因之间进行Wilcoxon秩和检验，但只涉及每个细胞中按表达量排名靠前的基因

all.sets <- lapply(names(markers.z), function(x) { GeneSet(markers.z[[x]], setName=x) })all.sets <- GeneSetCollection(all.sets)rankings <- AUCell_buildRankings(counts(sce.tasic), plotStats = FALSE, verbose = FALSE)cell.aucs <- AUCell_calcAUC(all.sets, rankings)## Genes in the gene sets NOT available in the dataset: ## endothelial-mural: 7 (8% of 83)## interneurons: 1 (1% of 118)## oligodendrocytes: 2 (2% of 81)## pyramidal CA1: 4 (3% of 125)## pyramidal SS: 4 (3% of 149)results <- t(assay(cell.aucs))head(results)## gene sets## cells astrocytes_ependymal endothelial-mural interneurons## Calb2_tdTpositive_cell_1 0.1386528 0.04264310 0.5306168## Calb2_tdTpositive_cell_2 0.1366283 0.04884635 0.4538432## Calb2_tdTpositive_cell_3 0.1087323 0.07269892 0.3459148## Calb2_tdTpositive_cell_4 0.1321658 0.04993108 0.5112665## Calb2_tdTpositive_cell_5 0.1512892 0.07161420 0.4929771## Calb2_tdTpositive_cell_6 0.1342012 0.09161375 0.3378310## gene sets## cells microglia oligodendrocytes pyramidal CA1## Calb2_tdTpositive_cell_1 0.04845394 0.1317958 0.2317668## Calb2_tdTpositive_cell_2 0.02682648 0.1211293 0.2062570## Calb2_tdTpositive_cell_3 0.03583482 0.1566986 0.3218726## Calb2_tdTpositive_cell_4 0.05387632 0.1480893 0.2546714## Calb2_tdTpositive_cell_5 0.06655747 0.1385766 0.2088478## Calb2_tdTpositive_cell_6 0.03201310 0.1552946 0.4010508## gene sets## cells pyramidal SS## Calb2_tdTpositive_cell_1 0.3476778## Calb2_tdTpositive_cell_2 0.2762466## Calb2_tdTpositive_cell_3 0.5244614## Calb2_tdTpositive_cell_4 0.3505890## Calb2_tdTpositive_cell_5 0.3009582## Calb2_tdTpositive_cell_6 0.5392798

# 对测试数据集中的每个细胞分配细胞类型进行注释，将AUC分数最高的标记基因集作为该细胞的标签

new.labels <- colnames(results)[max.col(results)]tab <- table(new.labels, sce.tasic$broad_type)tab## ## new.labels Astrocyte Endothelial Cell GABA-ergic Neuron## astrocytes_ependymal 43 2 0## endothelial-mural 0 27 0## interneurons 0 0 759## microglia 0 0 0## oligodendrocytes 0 0 1## pyramidal SS 0 0 1## ## new.labels Glutamatergic Neuron Microglia Oligodendrocyte## astrocytes_ependymal 0 0 0## endothelial-mural 0 0 0## interneurons 2 0 0## microglia 0 22 0## oligodendrocytes 0 0 38## pyramidal SS 810 0 0## ## new.labels Oligodendrocyte Precursor Cell Unclassified## astrocytes_ependymal 20 4## endothelial-mural 0 2## interneurons 0 15## microglia 0 1## oligodendrocytes 2 0## pyramidal SS 0 60

# 检查每个标签的AUC在细胞间的分布情况。在异质群体中，每个标签的分布应该呈双峰出现，一个高分峰包含该细胞类型的细胞，一个低分峰包含其他类型的细胞。这两个峰值之间的差距可以用来得出一个阈值，即一个标签对一个特定的细胞是否 "活跃"。在这种情况下，我们只需取每个细胞中得分最高的一个标签，因为这些标签之间应该是相互排斥）。在预期有特定细胞类型的人群中，相应标签缺乏明显的双峰性，可能表明其基因集的信息量不够大

par(mfrow=c(3,3))AUCell_exploreThresholds(cell.aucs, plotHist = TRUE, assign = TRUE)## $astrocytes_ependymal## $astrocytes_ependymal$aucThr## $astrocytes_ependymal$aucThr$selected## R_k3 ## 0.5804594 ## ## $astrocytes_ependymal$aucThr$thresholds## threshold nCells## Global_k1 0.2156583 94## L_k2 0.2187295 94## R_k3 0.5804594 43## ## $astrocytes_ependymal$aucThr$comment## [1] ""## ## ## $astrocytes_ependymal$assignment

# 当标记物集是相互排斥的时候，AUCell结果的解释是最直接的，比如上面提到的细胞类型标记物。在任何给定的细胞中，可以使用AUC来量化这种细胞的活跃程度，以便在不同的细胞中进行比较，但是，这种结果也可能会出现假阳性，所以需要我们谨慎处理。我们前面也说了，AUCell方法的优势在于它不需要参考表达值，这对于在处理来自文献或其他定性形式的生物知识的基因集时就特别有用，比如，我们可以使用从MSigDB中定义的单细胞标签：

library(BiocFileCache)## 载入需要的程辑包：dbplyrbfc <- BiocFileCache(ask=FALSE)## Warning: DEPRECATION: As of BiocFileCache (>1.15.1), default caching location has changed.## Problematic cache: C:\Users\Welford\AppData\Local/BiocFileCache/BiocFileCache/Cache

## See https://www.bioconductor.org/packages/devel/bioc/vignettes/BiocFileCache/inst/doc/BiocFileCache.html#default-caching-location-update

scsig.path <- bfcrpath(bfc, file.path("http://software.broadinstitute.org", "gsea/msigdb/supplemental/scsig.all.v1.0.symbols.gmt"))scsigs <- getGmt(scsig.path)# 用Muraro数据集进行标注load(file = "muraro.Rda")muraro.mat <- counts(sce.muraro)rownames(muraro.mat) <- rowData(sce.muraro)$symbolmuraro.rankings <- AUCell_buildRankings(muraro.mat, plotStats=FALSE, verbose=FALSE)scsig.aucs <- AUCell_calcAUC(scsigs, muraro.rankings)scsig.results <- t(assay(scsig.aucs))full.labels <- colnames(scsig.results)[max.col(scsig.results)]tab <- table(full.labels, sce.muraro$label)fullheat <- pheatmap(log10(tab + 1), color = viridis::viridis(100), silent=TRUE)scsigs.sub <- scsigs[grep("Pancreas", names(scsigs))]sub.aucs <- AUCell_calcAUC(scsigs.sub, muraro.rankings)## Genes in the gene sets NOT available in the dataset: ## Muraro_Pancreas_Alpha_Cell: 39 (7% of 566)## Muraro_Pancreas_Beta_Cell: 50 (5% of 948)## Muraro_Pancreas_Delta_Cell: 14 (6% of 249)## Muraro_Pancreas_Pancreatic_Polypeptide_Cell: 9 (6% of 155)## Muraro_Pancreas_Epsilon_Cell: 1 (2% of 44)## Muraro_Pancreas_Ductal_Cell: 58 (5% of 1277)## Muraro_Pancreas_Acinar_Cell: 20 (3% of 731)## Muraro_Pancreas_Mesenchymal_Stromal_Cell: 25 (4% of 680)## Muraro_Pancreas_Endothelial_Cell: 18 (5% of 362)sub.results <- t(assay(sub.aucs))sub.labels <- colnames(sub.results)[max.col(sub.results)]tab <- table(sub.labels, sce.muraro$label)subheat <- pheatmap(log10(tab + 1), color = viridis::viridis(100), silent=TRUE)gridExtra::grid.arrange(fullheat[[4]], subheat[[4]])

# 当然从这里看的话，这个标记的结果并不理想，不过细胞注释的过程本来就是需要不断调整的过程

好了，本期的内容到此结束，在本期中，我介绍了单细胞数据注释的内容，也简单了解了两种不同的单细胞注释方式。当然，对于自己的数据集以及聚类结果，我们在进行单细胞注释的时候需要用多种不同的方法，甚至多种不同的标记基因进行比对注释，再从中找到最符合我们预期的结果，这些就需要大家根据自己的数据进行调整并且选择合适的方法啦！

后台回复“feng 53”获取代码和数据，我们下期见！

单细胞分析专栏传送门

大年初一，单细胞分析从nature communication开始

8+单细胞文章从解读到套路总结全部送给你

30分多的单细胞分析文章是怎么做的？

不用R，点点点就能get的单细胞公共数据下载，你“学废”了吗

单细胞都有两个“对象”，你一个都没有

8成生信审稿人必问！单细胞第一张可视化图，你漏了吗？

我从Science中，“偷”到了这个实用的聚类分析技能！真舍不得分享！（附代码）

不骗你！学了这个R包，搞定九成高大上的单细胞分析！包治百病！一包封神！小白赶紧收起来！

Cell的生信分析能有多高级？单细胞研究必备分析！不会做么？简单，我来教你！（附代码）

放大招！单细胞筛分子标志物，让你的课题逼格up! up! up! 你一定要试试！（附详细代码）