3.7 细胞类型注释

刘小泽写于2020.7.5

1 前言

scRNA数据分析具有挑战性的工作一般就是对结果的解读。获得不同亚群的细胞很简单,但如何给每个亚群赋予生物学意义就比较困难,既需要对数据掌握,有需要有生物背景,而生物背景知识更新并不是那么快。我们现在掌握的,可能和真实水平还有很大偏差。

为了减弱自身生物背景的不足,可以用一些先验知识帮助判断。最常见的就是通过去了解参与特定生物学过程的经过验证的基因,利用GO、KEGG数据库;另外还可以将我们的表达谱与之前发布的参考数据集做对比,因为参考数据集中样本或细胞的生物学状态都经过了专家的证实。这两种形式(参考数据集、基因集)接下来都会探讨。

准备数据

下面2.2部分会用到

# 准备数据
library(scRNAseq)
sce.seger <- SegerstolpePancreasData()
# 基因注释
library(AnnotationHub)
edb <- AnnotationHub()[["AH73881"]]
symbols <- rowData(sce.seger)$symbol
ens.id <- mapIds(edb, keys=symbols, keytype="SYMBOL", column="GENEID")
ens.id <- ifelse(is.na(ens.id), symbols, ens.id)
# 去除重复行
keep <- !duplicated(ens.id)
sce.seger <- sce.seger[keep,]
rownames(sce.seger) <- ens.id[keep]
# 样本注释
emtab.meta <- colData(sce.seger)[,c("cell type",
"individual", "single cell well quality")]
colnames(emtab.meta) <- c("CellType", "Donor", "Quality")
colData(sce.seger) <- emtab.meta
sce.seger$CellType <- gsub(" cell", "", sce.seger$CellType)
sce.seger$CellType <- paste0(
toupper(substr(sce.seger$CellType, 1, 1)),
substring(sce.seger$CellType, 2))
# 质控
low.qual <- sce.seger$Quality == "low quality cell"
library(scater)
stats <- perCellQCMetrics(sce.seger)
qc <- quickPerCellQC(stats, percent_subsets="altexps_ERCC_percent",
batch=sce.seger$Donor,
subset=!sce.seger$Donor %in% c("HP1504901", "HP1509101"))
sce.seger <- sce.seger[,!(qc$discard | low.qual)]
# 归一化
# 看到quickCluster和computeSumFactors就知道使用的是去卷积化方法
library(scran)
clusters <- quickCluster(sce.seger)
sce.seger <- computeSumFactors(sce.seger, clusters=clusters)
sce.seger <- logNormCounts(sce.seger)

2 使用参考数据集

非常直接的方法,将我们的表达量与参考数据集的表达量去比对,看看和参考的样本有多少类似

这个牵扯到一个分类的问题,可以辅助以标准的机器学习(比如随机森林、支持向量机),任何之前发表过的RNA-Seq数据(bulk或single cell)中对细胞的划分都可以提供帮助。

这里主要是用了2019年发表的细胞注释R包SingleR (Aran et al. 2019) 。它会计算我们数据中的细胞与参考数据集中的细胞之间Spearman相关性,找到最相关的属性赋予我们的细胞。为了减少数据的噪音,这个R包先鉴定marker基因,然后根据marker基因的表达量去计算相关性。

2.1 使用内置的参考注释数据

SingleR中就包含了一些内置数据集,大部分是bulk RNA-Seq或芯片数据中经过筛选的细胞类型。

下面使用的参考数据集就来自Blueprint和ENCODE计划:(Martens and Stunnenberg 2013; The ENCODE Project Consortium 2012).

# 我们的数据
sce.pbmc
## class: SingleCellExperiment
## dim: 33694 3922
## metadata(1): Samples
## assays(2): counts logcounts
## rownames(33694): RP11-34P13.3 FAM138A ... AC213203.1 FAM231B
## rowData names(2): ID Symbol
## colnames(3922): AAACCTGAGAAGGCCT-1 AAACCTGAGACAGACC-1 ...
## TTTGTCACAGGTCCAC-1 TTTGTCATCCCAAGAT-1
## colData names(4): Sample Barcode sizeFactor label
## reducedDimNames(3): PCA TSNE UMAP
## altExpNames(0):
# 参考数据集
library(SingleR)
ref <- BlueprintEncodeData()
# 然后把我们的细胞在参考数据集中找对应的细胞类型
# 返回的pred结果是一个数据框,每行是我们自己数据的一个细胞
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
## 525 755 1254 1 1 5
## HSC Monocytes NK cells
## 14 1116 251
# 可视化
plotScoreHeatmap(pred)

横坐标是各种细胞类型,纵坐标是我们数据中的细胞;Score也是经过归一化后的,介于0-1

理论上,我们数据中的每个细胞都会有自己对应的细胞类型,并且在所属的细胞类型中,比其他细胞相关性更高(也就是得分更高、颜色差异更明显)。看到只有B-cells、monocytes表现突出,与其他细胞类型区分明显。而像CD4+、CD8+细胞,颜色相近,区分不明显

注意

当然,会存在一些细胞归属不是很清楚的,好像一个细胞对两种细胞类型都有点沾边,但又不好判断。函数会通过一个阈值将这些划分不清楚的细胞定为”NA”

sum(is.na(pred$pruned.labels))
## [1] 76
plotScoreDistribution(pred)

这个函数会把三种情况一起画出来,并且顺序依次是:匹配上的(黄色 assigned)、根本就不匹配的(灰色 other)、不清不楚的(橙色 pruned)

可以看到大部分都是灰色,没有匹配上的

另外,还可以看cluster对label的对应关系。理想情况是,一个cluster只对应一个label,但现实是:可能好几个clusters属于一个label;还有两个label的clusters很相似(例如CD4+、CD8+)

tab <- table(Assigned=pred$pruned.labels, Cluster=colLabels(sce.pbmc))
library(pheatmap)
pheatmap(log2(tab+10), color=colorRampPalette(c("white", "blue"))(101))

可以看到,使用已知的细胞注释可以简化细胞的生物学解释过程,但同时也会将细胞们限定在已知的框架中。如果发现我们聚类分群的结果与注释后的结果相差很大,不用着急去掉某种结果,留着下游再继续探讨差异来源

2.2 使用自定义的注释数据

比如手上哟一个之前注释过的数据,现在想用其中的注释信息对另一个数据进行注释

例如,之前有一个注释好的人类胰腺数据集 sce.muraro Muraro et al. (2016)

## class: SingleCellExperiment
## dim: 16940 2299
## metadata(0):
## assays(2): counts logcounts
## rownames(16940): ENSG00000268895 ENSG00000121410 ... ENSG00000159840
## ENSG00000074755
## rowData names(2): symbol chr
## colnames(2299): D28-1_1 D28-1_2 ... D30-8_93 D30-8_94
## colData names(4): label donor plate sizeFactor
## reducedDimNames(0):
## altExpNames(1): ERCC

去掉那些不明确的细胞类型

sce.muraro <- sce.muraro[,!is.na(sce.muraro$label) &
sce.muraro$label!="unclear"]
table(sce.muraro$label)
##
## acinar alpha beta delta duct endothelial
## 217 795 442 189 239 18
## epsilon mesenchymal pp
## 3 80 96

现在有一个新的数据集sce.seger Segerstolpe et al. (2016),它本身也包括细胞注释,在数据准备代码中有提及:在sce.seger$CellType,这个注释是领域专家设定的

# 以这个新数据集作为测试数据,sce.muraro作为参考数据,看看注释结果如何
pred.seger <- SingleR(test=sce.seger, ref=sce.muraro,
labels=sce.muraro$label, de.method="wilcox")
table(pred.seger$labels)
##
## acinar alpha beta delta duct endothelial
## 188 889 279 105 385 17
## epsilon mesenchymal pp
## 5 53 169

用这个结果与测试数据本身带的细胞注释结果进行比较:

tab <- table(pred.seger$pruned.labels, sce.seger$CellType)
library(pheatmap)
pheatmap(log2(tab+10), color=colorRampPalette(c("white", "blue"))(101))

可以看到这两种注释结果比较相近

2.3 使用marker基因集

这部分可以作为2.2的延伸 它基于的想法是:高表达的marker基因与每个细胞的生物学功能相关,因此可以用marker基因作为桥梁,连接测试数据集与参考数据集

基于marker基因的方法会更快捷,同时也不受限于参考数据集是否有专家定义好的细胞类型

参考数据集 sce.zeisel Zeisel et al. (2015)

目的是将数据集的marker基因于对应的细胞类型联系起来

需要注意的是:下面的pairwiseWilcox()代码需要基于R 4.0,需要加载归一化后的的数据,前面的数据处理在之前章节都有提及 数据在:https://share.weiyun.com/fDePnIGC

load('normalized.sce.zeisel.RData')
> table(sce.zeisel$level1class)
astrocytes_ependymal endothelial-mural
179 160
interneurons microglia
290 78
oligodendrocytes pyramidal CA1
774 938
pyramidal SS
397
library(scran)
wilcox.z <- pairwiseWilcox(sce.zeisel, sce.zeisel$level1class,
lfc=1, direction="up")
markers.z <- getTopMarkers(wilcox.z$statistics, wilcox.z$pairs,
pairwise=FALSE, n=50)
> lengths(markers.z)
astrocytes_ependymal endothelial-mural
79 83
interneurons microglia
118 69
oligodendrocytes pyramidal CA1
81 125
pyramidal SS
149

测试数据集 sce.tasic Tasic et al. (2016)

library(scRNAseq)
sce.tasic <- TasicBrainData()
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):
## altExpNames(1): ERCC

接下来使用AUCell 这个包

首先拿到参考数据集中的marker基因集markers.z,并构建一个GeneSetCollection对象

然后对测试数据集中基因表达量进行排名,最后根据排名和之前的基因集进行统计,得到每个marker基因集在每个细胞中的area under the curve (AUC)指标,就能得出哪一个marker基因集在测试数据集的细胞中占比最高

> names(markers.z)
[1] "astrocytes_ependymal" "endothelial-mural"
[3] "interneurons" "microglia"
[5] "oligodendrocytes" "pyramidal CA1"
[7] "pyramidal SS"
# 对marker基因集进行操作
library(GSEABase)
all.sets <- lapply(names(markers.z), function(x) {
GeneSet(markers.z[[x]], setName=x)
})
all.sets <- GeneSetCollection(all.sets)
> all.sets
GeneSetCollection
names: astrocytes_ependymal, endothelial-mural, ..., pyramidal SS (7 total)
unique identifiers: Apoe, Clu, ..., Gabrb2 (555 total)
types in collection:
geneIdType: NullIdentifier (1 total)
collectionType: NullCollection (1 total)
# 对测试数据集进行操作
library(AUCell)
# 基因表达量排名
rankings <- AUCell_buildRankings(counts(sce.tasic),
plotStats=FALSE, verbose=FALSE)
# 计算AUC(结果见下图)
cell.aucs <- AUCell_calcAUC(all.sets, rankings)
# 转置
results <- t(assay(cell.aucs))

这样就能清楚看到,每个细胞对应最高的AUC值属于哪个marker基因集,因此可以作为这个细胞的label

结果的检查 ——方法一:与参考数据中的cluster类型作比较

new.labels <- colnames(results)[max.col(results)]
table(new.labels, sce.tasic$broad_type)

可以看到一个cluster中的大部分细胞,都被赋予了一个相同的label

结果的检查 ——方法二:计算AUC最大值减去中位数

先说判断方法:如果注释的很明确,这个差值会较大;较小的差值说明结果不理想;当然我们可以自己定义一个阈值,过滤掉这些无法注释的细胞

下面的例子中,就是假定得到的大多都是注释较好的细胞,然后得到这个阈值

library(scater)
library(DelayedArray)
deltas <- rowMaxs(results) - rowMedians(results)
discard <- isOutlier(deltas, type="lower", batch=new.labels)
table(new.labels[discard])
##
## astrocytes_ependymal endothelial-mural interneurons
## 24 1 7
## oligodendrocytes pyramidal SS
## 10 16

看看这些被过滤掉的,在整体分布中处于什么位置

par(mar=c(10,4,1,1))
boxplot(split(deltas, new.labels), las=2)
points(attr(discard, "thresholds")[1,], col="red", pch=4, cex=2)

AUCell 包的好处在于不依赖表达量

只需要获得基因列表(可以基于文献或其他背景知识获得),例如从MsigDb上下载

library(BiocFileCache)
bfc <- BiocFileCache(ask=FALSE)
scsig.path <- bfcrpath(bfc, file.path("http://software.broadinstitute.org",
"gsea/msigdb/supplemental/scsig.all.v1.0.symbols.gmt"))
scsigs <- getGmt(scsig.path)
> scsigs
GeneSetCollection
names: Zheng_Cord_Blood_C1_Putative_Megakaryocyte_Progenitor, Zheng_Cord_Blood_C2_Putative_Basophil_Eosinophil_Mast_cell_Progenitor, ..., Muraro_Pancreas_Endothelial_Cell (256 total)
unique identifiers: ABCC3, ABCC4, ..., SEC31A (18656 total)
types in collection:
geneIdType: NullIdentifier (1 total)
collectionType: NullCollection (1 total)

不基于表达量有好有坏,好处就是避免了表达量出现的噪音干扰,处理也更快;坏处是可能会丢掉发现更精细的亚群可能性。

再看一眼另一个数据,不需要运行重点看看差异即可

# 先对测试数据集基因排名
muraro.mat <- counts(sce.muraro)
rownames(muraro.mat) <- rowData(sce.muraro)$symbol
muraro.rankings <- AUCell_buildRankings(muraro.mat,
plotStats=FALSE, verbose=FALSE)
# 然后把MsigDb的scsigs应用在测试数据上
scsig.aucs <- AUCell_calcAUC(scsigs, muraro.rankings)
scsig.results <- t(assay(scsig.aucs))
# 准备作图数据1
full.labels <- colnames(scsig.results)[max.col(scsig.results)]
tab <- table(full.labels, sce.muraro$label)
fullheat <- pheatmap(log10(tab+10), color=viridis::viridis(100), silent=TRUE)
# 然后把MsigDb的一个子集(Pancreas相关)的scsigs应用在测试数据上
scsigs.sub <- scsigs[grep("Pancreas", names(scsigs))]
sub.aucs <- AUCell_calcAUC(scsigs.sub, muraro.rankings)
sub.results <- t(assay(sub.aucs))
# 准备作图数据2
sub.labels <- colnames(sub.results)[max.col(sub.results)]
tab <- table(sub.labels, sce.muraro$label)
subheat <- pheatmap(log10(tab+10), color=viridis::viridis(100), silent=TRUE)
# 组合
gridExtra::grid.arrange(fullheat[[4]], subheat[[4]])

看到一个明显的变化就是:取了MsigDb关于Pancreas的子集后,效果变好了;而使用全部的基因集会导致混淆。因此使用哪些基因还是要基于我们对自己测试数据集的认识

3 基于marker基因的富集分析

思路:还是先得到每个cluster的标志性marker基因,然后对这些代表该cluster的基因们进行富集分析,看看它们集中在哪些通路

示例数据

再以一个新的小鼠乳腺数据 Bach et al. (2017) 为例,还是经历了:基因ID注释=》质控=》归一化=》找HVGs=》降维=》聚类=》找marker基因 数据可以直接下载

#--- loading ---#
library(scRNAseq)
sce.mam <- BachMammaryData(samples="G_1")
#--- gene-annotation ---#
library(scater)
rownames(sce.mam) <- uniquifyFeatureNames(
rowData(sce.mam)$Ensembl, rowData(sce.mam)$Symbol)
library(AnnotationHub)
ens.mm.v97 <- AnnotationHub()[["AH73905"]]
rowData(sce.mam)$SEQNAME <- mapIds(ens.mm.v97, keys=rowData(sce.mam)$Ensembl,
keytype="GENEID", column="SEQNAME")
#--- quality-control ---#
is.mito <- rowData(sce.mam)$SEQNAME == "MT"
stats <- perCellQCMetrics(sce.mam, subsets=list(Mito=which(is.mito)))
qc <- quickPerCellQC(stats, percent_subsets="subsets_Mito_percent")
sce.mam <- sce.mam[,!qc$discard]
#--- normalization ---#
library(scran)
set.seed(101000110)
clusters <- quickCluster(sce.mam)
sce.mam <- computeSumFactors(sce.mam, clusters=clusters)
sce.mam <- logNormCounts(sce.mam)
#--- variance-modelling ---#
set.seed(00010101)
dec.mam <- modelGeneVarByPoisson(sce.mam)
top.mam <- getTopHVGs(dec.mam, prop=0.1)
#--- dimensionality-reduction ---#
library(BiocSingular)
set.seed(101010011)
sce.mam <- denoisePCA(sce.mam, technical=dec.mam, subset.row=top.mam)
sce.mam <- runTSNE(sce.mam, dimred="PCA")
#--- clustering ---#
snn.gr <- buildSNNGraph(sce.mam, use.dimred="PCA", k=25)
colLabels(sce.mam) <- factor(igraph::cluster_walktrap(snn.gr)$membership)
#--- find marker ---#
markers.mam <- findMarkers(sce.mam, direction="up", lfc=1)

然后我们想看cluster2是什么细胞类型

首先得到cluster2的marker基因,并且还挑差异显著的基因

chosen <- "2"
cur.markers <- markers.mam[[chosen]]
is.de <- cur.markers$FDR <= 0.05
summary(is.de)
## Mode FALSE TRUE
## logical 27819 179

进行GO富集分析

library(org.Mm.eg.db)
library(clusterProfiler)
gn=unique(rownames(cur.markers)[is.de])
trans=bitr(gn,fromType = "SYMBOL",
toType = 'ENTREZID',
OrgDb = org.Mm.eg.db)
go.out=enrichGO(gene = trans$ENTREZID,
keyType = 'ENTREZID',
OrgDb = org.Mm.eg.db,
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05,
readable = TRUE)
View(head(summary(go.out),20))

接下来就是考验自己对数据的认知了,假设我们这里对这个数据很熟悉。那么看到这里基因主要参与了lipid synthesis、 cell adhesion 、 tube formation。再结合乳腺的研究背景,推测cluster2可能包含与乳汁产生和分泌相关的管腔上皮细胞,进一步检查一下相关的基因Csn2和Csn3,也是上调的

plotExpression(sce.mam, features=c("Csn2", "Csn3"),
x="label", colour_by="label")

如果想看cluster2中某个感兴趣通路中的基因表现情况

# 比如想看与cell adhesion相关的通路
library(stringr)
adhesion <- go.out["GO:0022408",]
adhesion_gn=str_split(adhesion$geneID,"/",simplify = T)
head(cur.markers[rownames(cur.markers) %in% adhesion_gn,1:4], 10)
## DataFrame with 10 rows and 4 columns
## Top p.value FDR summary.logFC
## <integer> <numeric> <numeric> <numeric>
## Spint2 11 3.28234e-34 1.37163e-31 2.39280
## Epcam 17 8.86978e-94 7.09531e-91 2.32968
## Cebpb 21 6.76957e-16 2.03800e-13 1.80192
## Cd24a 21 3.24195e-33 1.29669e-30 1.72318
## Btn1a1 24 2.16574e-13 6.12488e-11 1.26343
## Cd9 51 1.41373e-11 3.56592e-09 2.73785
## Ceacam1 52 1.66948e-38 7.79034e-36 1.56912
## Sdc4 59 9.15001e-07 1.75467e-04 1.84014
## Anxa1 68 2.58840e-06 4.76777e-04 1.29724
## Cdh1 69 1.73658e-07 3.54897e-05 1.31265

TODO:补充12.5 Computing gene set activities