常规的的单细胞亚群命名应该是降维聚类分群后提前每个亚群top特异性高表达基因,以及检查已知标记基因在不同单细胞亚群的表达量情况,依据生物学背景合理命名。这样的话就需要大量阅读自己领域的相关文献,参考前人的单细胞降维聚类分群后的生物学命名以及文章里面的标记基因。
如果你的单细胞数据集里面的单细胞亚群是已知的,比如血液里面的PBMC里面的肯定是淋巴系(T,B,NK细胞)和髓系(单核,树突,巨噬,粒细胞)的两大类,在单细胞转录组技术流行之前,应该是有这些单细胞亚群各自的流式细胞仪筛选方法,也算是得到了不同单细胞亚群富集后的bulk转录组测序数据,就有表达量矩阵。比如大名鼎鼎的singleR包配套的参考数据集celldex包里面的:
ref <- MouseRNAseqData()
,是geo数据库里面的28种纯化好的不同单细胞亚群的bulk转录组测序数据ref <- NovershternHematopoieticData()
,是 microarray datasets for sorted hematopoietic cell populations from GSE24759 (Novershtern et al. 2011).ref <- MonacoImmuneData()
,是新加坡团队的bulk RNA-seq samples of sorted immune cell populations from GSE107011 (Monaco et al. 2019).
而singleR包其实就是依赖于这些已知的不同单细胞亚群的bulk转录组测序数据表达量矩阵来指导单细胞分群命名,这就是为什么很多小伙伴都觉得singleR注释不是很好用,因为这些数据集都是分散的, 而我们的单细胞转录组数据集里面的单细胞亚群是丰富多样的, 单独采用任何一个独立的参考数据集都是会有偏差的。
但是,如果你有两个相近的单细胞转录组数据集,其中一个数据集有比较好的注释了,就可以把这些注释好的单细胞转录组数据集做出Pseudo-bulk的表达量矩阵后继续指导另外一个还没有注释的单细胞转录组数据集。我们来举例说明:
基于singleR的相关性
首先需要把注释好的单细胞转录组数据集做出Pseudo-bulk的表达量矩阵
理论上我们应该是针对PBMC里面的不同单细胞亚群:
# install.packages('devtools')
# devtools::install_github('satijalab/seurat-data')
library(SeuratData) #加载seurat数据集
library(Seurat)
getOption('timeout')
options(timeout=10000)
# InstallData("pbmc3k")
data("pbmc3k")
sce <- pbmc3k.final
colnames(sce@meta.data)
sce$celltype = Idents(sce)
DimPlot(sce,group.by = "celltype",label = T,repel = T)
av <-AverageExpression(sce,
group.by = "celltype",
assays = "RNA")
av=av[[1]]
head(av)
write.csv(av,file = 'AverageExpression_celltype.csv')
cg=names(tail(sort(apply(av, 1, sd)),1000))
pheatmap::pheatmap(cor(av[cg,]))
pheatmap::pheatmap(cor(av[cg,]),
file = 'AverageExpression_celltype.pdf')
dev.off()
Ref=av
ref_sce=SingleCellExperiment::SingleCellExperiment(assays=list(counts=Ref))
library(scater)
ref_sce=scater::logNormCounts(ref_sce)
library(SingleCellExperiment)
logcounts(ref_sce)[1:4,1:4]
colData(ref_sce)$Type=colnames(Ref)
table(colnames(Ref))
ref_sce
save(ref_sce,file = 'ref_sce_by_pbmc3k_celltype.Rdata')
然后就可以依据这个把这些注释好的单细胞转录组数据集做出Pseudo-bulk的表达量矩阵,对原来的pbmc3k里面的接近三千多个细胞进行亚群命名。
library(Seurat)
library(ggplot2)
library(SingleR)
testdata <- GetAssayData(sce, slot="data")
testdata[1:4,1:4]
dim(testdata)
library(SingleR)
pred <- SingleR(test=testdata, ref=ref_sce,
labels=ref_sce$Type
)
as.data.frame(table(pred$labels))
head(pred)
labels=pred$labels
table(labels)
sce$singleR=labels
DimPlot(sce, reduction = "umap",
group.by = "singleR",label = T,repel = T) +
DimPlot(sce, reduction = "umap",group.by = "celltype", label = T,repel = T)
ggsave('umap_comp_by_singleR.pdf',width = 12)
gplots::balloonplot(table(sce$celltype,
sce$singleR ))
可以看到,准确度还是蛮不错的,比较麻烦的地方就是naive和memory的CD4细胞的区分:

naive和memory的CD4细胞的区分比较差其实是因为naive和memory的CD4细胞的区分在参考数据集里面本来就是没有边界线的,所以注释也不可能完美 :

假如你同时做某个癌症的十几个单细胞转录组数据集,那么这个基于singleR的相关性就可以派上用场啦,只需要完成其中任意一个,剩余的就可以靠相关性注释。
基于projectLSI
其中LSI是 latent semantic indexing (LSI)算法,这个projectLSI包的GitHub官网是:https://github.com/sajuukLyu/projectLSI,最早出现在文献:Single-cell multiomic analysis identifies regulatory programs in mixed-phenotype acute leukemia. Nat Biotechnol 37, 1458–1465 (2019)
但是我没有看到它可以比较好的把一个单细胞转录组数据集的单细胞亚群注释结果迁移到另外的其它单细胞数据集,如果大家有更好的工具也可以留言推荐哈。