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