细胞注释——SingleR包
Hello~大家好，还记得上期我们介绍过的CellMarker数据库进行细胞注释的流程吗？
毕竟嘛，手动注释尽管准确性更高，也更自由，但是耗时费力，而且在面对大量的细胞亚群的时候，有时候难免会造成流程上的繁琐，这里就给大家介绍一个进行细胞注释的老牌R包——SingleR包
晨曦解读
本教程参考bioconductor官网SingleR包帮助文档以及“The SingleR book”等多渠道关于SingleR资源，再保证质量的同时，希望本教程可以为SingleR的普及助力添瓦，为scRNA-seq细胞注释的普及有稍微一丢丢的助力，那将是我作为编者的荣幸~
首次揭秘！不做实验也能发10+SCI，CNS级别空间转录组套路全解析（附超详细代码！）
过关神助！99%审稿人必问，多数据集联合分析，你注意到这点了吗？
太猛了！万字长文单细胞分析全流程讲解，看完就能发文章！建议收藏！（附代码）
秀儿！10+生信分析最大的难点在这里！30多种方法怎么选？今天帮你解决！
引言
在scRNA-seq数据呈现井喷式的今天，过饱和的数据带给我们的问题就是我们无法有一个一致公认且准确高效的方式完成细胞注释，目前的细胞注释很大程度上依赖于marker gene，通常需要我们手工完成，这种方式具有直观性，但是耗时费力，且选择细胞亚群的时候具有主观性，限制了相关亚群的细致分化，所以有必要减少细胞注释所消耗的人力和物力，把我们的精力更多的集中在生物学差异上。
SingleR原理篇介绍
首先我们需要指定参考数据集(这个在SingleR包以及celldex包中都有参考数据集,后面会介绍),即已经知道细胞类型的数据,然后SingleR会基于这个数据来对测试数据(我们需要明确细胞类型的数据)进行细胞注释,流程如下:
1. 建立Marker基因集合
2. 计算测试数据集marker基因与参考基因中marker基因的相关性系数
3. 细胞注释
4. 注释结果诊断
晨曦解读
接下来咱们进行一步步的解释,首先第一步为建立marker基因集合,这一步比较繁琐,大家只需要知道,我这一步得到的结果是参考基因组各个已知细胞类型的marker基因即可
然后第二步就是用我们测试数据集的marker基因去比对参考基因组中的marker基因,之所以用marker基因,是因为这样我们可以减少非marker基因的干扰,使我们的结果更加精确且提高运算速度
随后第三步是进行细胞注释,这一步骤其实应该分为两种可能,测试数据集我给默认为X,X与参考基因组中的细胞类型A的相关性最高,所以我可以把其X注释为A(这个相关性就是X的marker基因的表达同细胞类型A中marker基因的表达相关性很高,即X里基因的种类和表达幅度都与细胞类型A高度相似);但是还有一种可能就是X与参考数据集中的细胞类型A,细胞类型B都很相似,这时候我们就需要进行新一轮的细胞注释,也就是所谓的微调
对于微调的的流程如下:
首先SingleR会逐步缩小范围,也就是说会不断提高相关性系数筛选的阈值,直到在参考数据集中只有两个细胞类型为止,并不是一步到位,而是步步为营
然后SingleR会将剩余的细胞类型重新构建参考数据集,回到第一轮重新进行构建marker基因,然后再次进行相关性系数的匹配,选择高的相关性系数作为我们对X的注释,这样做的好处显而易见,十分的可靠~
最后一步,我们需要对我们注释的结果进行判断,首先我们需要了解下面这个图
晨曦解读
一列是一个细胞,每一行为参考数据集里面的细胞标签,每一格表示细胞在该标签中所获得的打分,颜色代表分数的高低(这个打分是我们计算的相关性打分,对于注释结果比较明确的细胞,我们希望这个细胞的在该标签的打分是要显著高于其它标签的)
晨曦解读
每一个长方形格子表示一个 细胞类型,一个点则代表一个细胞,横坐标为分配到该类型的细胞,纵坐标为该细胞的delta值中位数(delta值:注释标签的分数与所有标签之间的分数(X与细胞类型A有关的标签分数)之间的差异),所以从定义来看delta值应该越高越好,这样我们注释的标签和别的标签才可以更好的区分开
晨曦解读
每一行为一个基因,每一列为一个细胞,labels表示的为细胞注释的结果,这个图的含义是,我们如果把X注释为细胞类型A,我们就认为细胞类型A中的marker基因应该在X中高表达,,即粗略的看就是每一个labels下的基因表达应该呈现一致,不同labels下的基因表达应该呈现分割分明
总结一下就是:建立Marker基因集合→计算测试数据集marker基因与参考基因中marker基因的相关性系数→细胞注释→注释结果诊断(可视化方式有很多,选择其中一个即可)
#相关学习http://bioconductor.org/books/release/SingleRBook/introduction.html#开源许可http://bioconductor.org/packages/release/bioc/vignettes/SingleR/inst/doc/SingleR.html#标准教程#相关文献Aran, D., A. P. Looney, L. Liu, E. Wu, V. Fong, A. Hsu, S. Chak, et al. 2019. “Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage.” Nat. Immunol. 20 (2): 163–72.Mabbott, Neil A., J. K. Baillie, Helen Brown, Tom C. Freeman, and David A. Hume. 2013. “An expression atlas of human primary cells: Inference of gene function from coexpression networks.” BMC Genomics 14. https://doi.org/10.1186/1471-2164-14-632.Zheng, G. X., J. M. Terry, P. Belgrader, P. Ryvkin, Z. W. Bent, R. Wilson, S. B. Ziraldo, et al. 2017. “Massively parallel digital transcriptional profiling of single cells.” Nat Commun 8 (January): 14049.
SingleR代码篇实战
#首先进行一下准备工作#已知singleR包的内置数据集存放在了celldex包中,所以我们现在就来加载这两个R包library(celldex)library(SingleR)BiocManager::install("scRNAseq")library(scRNAseq)#load("/home/data/refdir/singleR/ImmGenData.Rdata")#上一个代码是我个人下载完celldex参考数据集后,因为下载下来直接就是Rdata,所以可以这样直接使用,也算是另一种曲线救国的办法吧~
晨曦解读
注意这两个R包非常的不好安装(貌似scRNA-seq就没有好安装的R包?),所以多尝试几次,如果实在不行,就只能放弃这个方法注释,参考上一篇推文进行注释,效果也是差不多的
福利
我把我下载好的几个数据库放在了百度云网盘（常用的），回复”晨曦05“即可获得哦~
#解读一下celldex包中的七大内置数据集(调用数据集的方式都是相似的,调用结果就是会在环境中出现一个类似于Seurat的对象)ref <- HumanPrimaryCellAtlasData()#The HPCA reference consists of publicly available microarray datasets derived from human primary cell.Most of the labels refer to blood subpopulations but cell types from other tissues are also availableref <- BlueprintEncodeData()#The Blueprint/ENCODE reference consists of bulk RNA-seq data for pure stroma and immune cells generated by Blueprint and ENCODE projects ref <- MouseRNAseqData()#This reference consists of a collection of mouse bulk RNA-seq data sets downloaded from the gene expression omnibus.A variety of cell types are available, again mostly from blood but also covering several other tissuesref <- ImmGenData()#The ImmGen reference consists of microarray profiles of pure mouse immune cells from the project of the same name.This is currently the most highly resolved immune reference - possibly overwhelmingly so, given the granularity of the fine labelsref <- DatabaseImmuneCellExpressionData()#The DICE reference consists of bulk RNA-seq samples of sorted cell populations from the project of the same name.ref <- NovershternHematopoieticData()#The Novershtern reference (previously known as Differentiation Map) consists of microarray datasets for sorted hematopoietic cell populations from GSE24759ref <- MonacoImmuneData()#The Monaco reference consists of bulk RNA-seq samples of sorted immune cell populations from GSE107011 #更详细的学习链接见:http://www.bioconductor.org/packages/release/data/experiment/vignettes/celldex/inst/doc/userguide.html
晨曦解读
这里面每个数据集里面的类型都很多,这个就告诉我们对于同一种细胞类型,如果我们选择不同的参考数据集,其鉴定的结果很可能也就不一样,这里面就衍生出两种思路,一种就是慎重的选择参考数据集,一种就是多数据集联合鉴定~
#回归正轨,接下来我们选择调用的参考数据集sceM <- MuraroPancreasData()#我们惊讶的发现这个数据集并没有在我们上面中,其实这个是运用scRNAseq这个包中的两个人类胰腺数据集,目的是使用一个预先标记的数据集来注释其他未标记的数据集,并没有使用上面的数据集
晨曦解读
这里告诉我们,scRNAseq包内的数据集也可以被我们用来当作参考数据集,后面的流程基本都是类似的
#对参考数据集进行质控sceM <- sceM[,!is.na(sceM$label) & sceM$label!="unclear"] #然后我们来看一下这个数据
晨曦解读
基本上和我们Seurat对象很类似,但是这里面的数据结构基本上就不需要我们了解了,因为后面的代码都是一些很标准的流程,当然有一些大家也是知道的,比如说assays,int_metadata等
#对参考数据集进行标准化(不要把参考数据集想的太神秘,它其实就是一个已经注释好的数据集而已)library(scater)sceM <- logNormCounts(sceM)#查看一下参考数据集里面每类细胞注释的分布情况table(sceM$label)
#建立测试数据集,也就是我们要进行注释的数据集,这里我们只是简单选择100个细胞来进行操作#SingleR本身就是针对单细胞来进行注释,当然也有人用来对亚群进行注释(貌似这个更常见),但是这样可能会缺失掉部分生物学差异sceG <- GrunPancreasData()sceG <- sceG[,colSums(counts(sceG)) > 0] # Remove libraries with no counts.sceG <- logNormCounts(sceG) sceG <- sceG[,1:100]
晨曦解读
SingleR包大家可以把它当作一个流程包,就是说这个包只是给我们提供了一个流程,仅此而已,然后这个R包支持任何归一化后的数据,也就是说如果我们用我们自己scRNA-seq的数据时候,需要经过归一化,也就是Seurat标准流程聚类后的结果~
#进行注释pred.grun <- SingleR(test=sceG, ref=sceM, labels=sceM$label, de.method="wilcox")#方法默认即可#查看注释结果table(pred.grun$labels)
晨曦解读
可以看到成功注释了100个细胞~
#上面使用了scRNAseq中已经注释成功的数据,这个其实就是类似于假设我检测了两个组织,这两个组织本身相似,然后我把其中一个组织进行了scRNA-seq并且注释成功后,我就可以用这个数据来注释另一个没有注释的数据#下面我们用celldex包中的数据集来进行注释library(celldex)hpca.se <- HumanPrimaryCellAtlasData()hpca.se
晨曦解读
这里就是简单介绍了一下这个数据集~
#选择测试数据集library(scRNAseq)hESCs <- LaMannoBrainData('human-es')hESCs <- hESCs[,1:100]#标准的注释流程library(SingleR)pred.hesc <- SingleR(test = hESCs, ref = hpca.se, assay.type.test=1, labels = hpca.se$label.main)table(pred.hesc$labels)
晨曦解读
至此我们完成了使用内置数据集来进行细胞注释的过程~
##下面我们开始拿着我们一开始的数据集来进行一下实战#首先我们需要调取我们晨曦单细胞笔记中标准的聚类流程最后的结果数据scRNA <- readRDS("~/scRNA.rds")
晨曦解读
首先因为scRNA-seq可能是数据量比较大的缘故，传统的save函数和load函数并不能完全保留工作环境内的文件，所以这里我们需要使用其专属于scRNA-seq的函数，即saveRDS函数以及readRDS函数
#然后我们接下来进行标准的注释流程#首先因为是血液样本是，所以我们调用参考数据集中关于血液的参考数据集#这里我们可以尝试多参考数据集进行注释，反正也是走标准流程#首先我们加载需要的R包library(SingleR)#加载参考数据集img.se <- load("/home/data/refdir/singleR/ImmGenData.Rdata")hpca.se <- load("/home/data/refdir/singleR/HumanPrimaryCellAtlasData.Rdata")#这里简单的赋值并不可以直接调用，只是在工作环境中做一个标签提醒作用ref@colData
晨曦解读
有的时候代码只有跑一遍，带入自己的数据，才会知道有些细小的问题，这里我们将在细致的捋一遍
首先这里我们看一下参考数据集中的colData，之所以看这个，是因为这里面就是细胞注释信息，其中的label.main和label.fine大家可以简单理解为，后者比前者的注释分辨率更高，也就是更细致~
下面是官网对参考数据集中标签的解读，其实大致翻译过来就是我上面说的意思
#然后我们导入测试数据集，也就是我们自己的数据pbmc<-readRDS("scRNA.rds")#获取测试数据集中标化后的数据（singleR的输入数据需要是RNA中的表达矩阵，如果对于多数据集联合，也是如此）#https://www.biostars.org/p/442514/（原始回答如下，感觉这种探讨分析流程的感觉棒棒的~）pbmcdata <- GetAssayData(pbmc, slot="data")#注释Anno<-SingleR(test = pbmcdata, ref =list(HP=ref,IM=ref4), labels = list(ref$label.fine,ref4$label.fine), method = "cluster", clusters = [email protected]$seurat_clusters)Anno
晨曦解读
然后我们可以看一下注释的结果，相信很多同学也和我一样，下面有些数据压根就不知道有什么用处，当然其实我们也不用知道，画箭头两个数据，一个是亚群标签，一个是注释的标签，我们只需要把这两个提取出来，即可比较直观的看到每个亚群的亚群注释~
celltype = data.frame(ClusterID=rownames(Anno), celltype=Anno$labels, stringsAsFactors = F)celltype
晨曦解读
至此细胞注释完成，接下来我们需要把注释信息添加到tSNE上
#首先给我们Seurat对象中添加细胞注释信息[email protected]$celltype = "NA"#先挖一个坑，然后用后面的循环进行补充for(i in 1:nrow(celltype)){ [email protected][which([email protected]$seurat_clusters == celltype$ClusterID[i]),'celltype'] <- celltype$celltype[i]}#可视化p1 = DimPlot(pbmc, group.by="celltype", label=T, label.size=5, reduction='tsne')p1
晨曦解读
然后这里我们其实可以发现一个问题，就是为什么我们的注释信息这么的长，然后我对其中的一个细胞亚群注释信息进行了搜索
HSC_-G-CSF：造血干细胞-粒细胞集落刺激因子（好吧，看来还是我平时的基础知识不扎实~）
#其他可视化展示select_genes <- c('LYZ','CD79A')p1 <- FeaturePlot(pbmc, features = select_genes, reduction = "tsne", label=T, ncol=2)p2 <- VlnPlot(pbmc, features = select_genes, pt.size=0, group.by="celltype", ncol=2)
晨曦解读
第二个小提琴图可能会因为横坐标太长而导致无法完全展开，目前来看VInPlot函数并没有调节坐标轴大小的参数，所以这个目前减少注释长度，比如说根据注释进行进行自主注释然后替换到celltype上，流程跟上面一样，然后进行可视化，或者是用没注释前的细胞标号来作为横坐标~
至此标准的SingleR流程到这里就讲解完毕了，这里需要有几点强调的，需要注意
晨曦解读
参考数据集有三个标签维度，这里我们只需要知道label.fine精细，label.main比前面的稍微粗糙
注释成功后会也会有三个层次的标签，分别是微调前（first.labels）、微调后（labels）以及修剪后（pruned.labels）这里默认选择labels即可
晨曦解读
SingleR也只是一种细胞注释的方法，标准的流程就是这些，但是依旧需要我们进行主观的判断和选择，包括参考数据集的选择，多种注释联合等等，还要考虑到可视化后横坐标注释内容的长度是否可以展开，毕竟如果要在VInPlot函数中调节确实是没有太好的办法的，虽然现在支持ggplot2但是调节毕竟是有限的，其实可以提取表达矩阵后然后用ggplot2来画图，这个也是一种方法
#首先我们做一下准备library(Seurat)library(SeuratObject)#然后我们运用FetchData函数来提取LYZ基因在各个细胞中的表达并且把其转换成数据框的形式exprs <- data.frame(FetchData(object = scRNA, vars = "LYZ"))
晨曦解读
FetchData函数是可以针对我们seurat对象来抓取变量的函数，我们可以选择Gene的名称、或者nFeature_RNA（基因的数量）、nCount_RNA（RNA的表达量）等等，只要你可以在Seurat对象中找到都可以通过这个函数单独抓取出来，其抓取的形式就是变量和细胞ID之间的关系，类似下面这样~
#为了绘制小提琴图等相关ggplot2图形，我们知道了表达量后还需要知道分组信息，也就是亚群信息exprs$Barcod<-rownames(exprs)
ident<-data.frame()ident<-data.frame(Barcod=names(Idents(object = scRNA)),orig.ident=Idents(object = scRNA))#获取细胞ID和聚类信息的对应关系#合并c<-merge(exprs,ident,by='Barcod')
#因为我们后续绘图是需要分组信息的，这些分组信息用因子表示比较恰当，所以我们这里把其因子话，并且规定因子的水平c$orig.ident<-factor(c$orig.ident,levels=c(sort(unique(Idents(object = scRNA)))))#并不是排列大小，而是一步新惯性操作，比如说分组信息男女无法比较大小，但是我们还希望后续我们可以直接调用男and女，而不是女and男，我们就可以设置因子水平，然后sort()即可调用男and女
晨曦解读
有的同学可能对因子这个概念不是特别明白，对我上面说的话，可能也是一知半解，这里举个例子吧
p <- c("男","女")sort(p)
p <- factor(p,levels = c("男","女"))sort(p)
晨曦解读
通过上面的例子不知道同学们有没有明白，针对可以比较大小的数来说sort是可以精准排序的，但是对于不可以比较的字符来说，设置因子的水平可以帮助我们更好的调用信息
有的同学问这个具体有什么用处？
基本上对于我们进行差异分析的时候，是实验组比对照组，也就是cancer VS normal，我们可以在分组的时候设置因子水平，然后后期比较设置谁比谁的时候直接调用这个因子水平，这样就不至于我们陷入不确定是cancer比normal，还是normal比cancer的尴尬境地之中了~
#回到正题#我们已经获得了单个基因在各个细胞中的表达量并且已经和亚群（分组）信息对应成功#下一步我们就可以开始画图了#首先准备工作library(ggplot2)#标准流程绘图ggplot(data = c,mapping = aes(x = factor(x = orig.ident),y = LYZ)) +geom_violin(scale = "width",adjust =1,trim = TRUE,mapping = aes(fill = factor(x = orig.ident)))
晨曦解读
这个通过查看网上其他相关资料，如果用seurat中的VInPlot函数绘制同类型的图会包含一个降噪的处理，也就是所谓的把一些可能不怎么高的表达干脆给抹去，展现最直接的差异结果，但是我个人感觉加不加上这一步都是可以，因为我们绘制小提琴图的目的只有一个，就是看我们关注的变量在各个细胞亚群之间的表达差异即可，当然加上降噪也不算错误
而且我上面的这个图有一个问题，肯定是不可以直接放在文章中的，不知道同学们有没有发现？
#去除NAc <- c[!(rowSums(is.na(c))),]#判断有无NAanyNA(c)#FALSE#然后我们把这个表达数据带入到上面的流程中即可得到下图，其实也算是数据质控的一环
#这里也把降噪的代码给大家，可以用来学习，代码来自Seurat官网noise <- rnorm(n = length(x = data[,c('LYZ')])) / 100000data[,c('LYZ')] <- data[, c('LYZ')] + noiseggplot(data = c,mapping = aes(x = factor(x = orig.ident),y = LYZ)) +geom_violin(scale = "width",adjust =1,trim = TRUE,mapping = aes(fill = factor(x = orig.ident))) ggsave('test.pdf')
然后后续的绘图，我们就可以基于ggplot2来进行了，具体相关的绘图参数，其实我也没有记住太多，都是现绘图现去查ggplot2的帮助文档即可~
到这里今天的singleR教程从原理到实践再到绘图就结束了，尽可能做到了很全面，也是参考了很多网上的教程以及自己的感悟，希望可以帮助到大家在细胞注释这个难点中可以走到更加顺利~
但是，对于各种scRNA-seq的提取命令，相信同学们有些还不是很清楚，那么这个就是咱们下期的内容啦~
我是晨曦，我们下次再见~
回复“晨曦05”，即可获得今天的范例数据和代码哦~
—
END
—

撰文丨晨   曦
排版丨四金兄
主编丨小雪球
欢迎大家关注解螺旋生信频道-挑圈联靠公号~
继续阅读
阅读原文