# 3.7 细胞类型注释 ## 1 前言 scRNA数据分析具有挑战性的工作一般就是对结果的解读。获得不同亚群的细胞很简单，但如何给每个亚群赋予生物学意义就比较困难，既需要对数据掌握，有需要有生物背景，而生物背景知识更新并不是那么快。我们现在掌握的，可能和真实水平还有很大偏差。 为了减弱自身生物背景的不足，可以用一些先验知识帮助判断。最常见的就是通过去了解参与特定生物学过程的经过验证的基因，利用GO、KEGG数据库；另外还可以将我们的表达谱与之前发布的参考数据集做对比，因为参考数据集中样本或细胞的生物学状态都经过了专家的证实。这两种形式（参考数据集、基因集）接下来都会探讨。 ### **准备数据** 下面2.2部分会用到 ```r # 准备数据 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`](https://bioconductor.org/packages/release/bioc/vignettes/SingleR/inst/doc/SingleR.html) ([Aran et al. 2019](https://pubmed.ncbi.nlm.nih.gov/30643263/)) 。它会计算我们数据中的细胞与参考数据集中的细胞之间Spearman相关性，找到最相关的属性赋予我们的细胞。为了减少数据的噪音，这个R包先鉴定marker基因，然后根据marker基因的表达量去计算相关性。 ### **2.1 使用内置的参考注释数据** SingleR中就包含了一些内置数据集，大部分是bulk RNA-Seq或芯片数据中经过筛选的细胞类型。 下面使用的参考数据集就来自Blueprint和ENCODE计划：([Martens and Stunnenberg 2013](https://pubmed.ncbi.nlm.nih.gov/24091925/); [The ENCODE Project Consortium 2012](https://pubmed.ncbi.nlm.nih.gov/22955616/)). ```r # 我们的数据 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) ``` ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-05-083229.png) 横坐标是各种细胞类型，纵坐标是我们数据中的细胞；Score也是经过归一化后的，介于0-1 理论上，我们数据中的每个细胞都会有自己对应的细胞类型，并且在所属的细胞类型中，比其他细胞相关性更高（也就是得分更高、颜色差异更明显）。看到只有B-cells、monocytes表现突出，与其他细胞类型区分明显。而像CD4+、CD8+细胞，颜色相近，区分不明显 **注意** 当然，会存在一些细胞归属不是很清楚的，好像一个细胞对两种细胞类型都有点沾边，但又不好判断。函数会通过一个阈值将这些划分不清楚的细胞定为”NA” ```r sum(is.na(pred$pruned.labels)) ## [1] 76 plotScoreDistribution(pred) ``` 这个函数会把三种情况一起画出来，并且顺序依次是：匹配上的（黄色 assigned）、根本就不匹配的（灰色 other）、不清不楚的（橙色 pruned） 可以看到大部分都是灰色，没有匹配上的 ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-06-070758.png) 另外，还可以看cluster对label的对应关系。理想情况是，一个cluster只对应一个label，但现实是：可能好几个clusters属于一个label；还有两个label的clusters很相似（例如CD4+、CD8+） ```r tab <- table(Assigned=pred$pruned.labels, Cluster=colLabels(sce.pbmc)) library(pheatmap) pheatmap(log2(tab+10), color=colorRampPalette(c("white", "blue"))(101)) ``` ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-06-072204.png) 可以看到，使用已知的细胞注释可以简化细胞的生物学解释过程，但同时也会将细胞们限定在已知的框架中。如果发现我们聚类分群的结果与注释后的结果相差很大，不用着急去掉某种结果，留着下游再继续探讨差异来源 ### **2.2 使用自定义的注释数据** 比如手上哟一个之前注释过的数据，现在想用其中的注释信息对另一个数据进行注释 例如，之前有一个注释好的人类胰腺数据集 `sce.muraro` [Muraro et al. (2016)](https://pubmed.ncbi.nlm.nih.gov/27693023/) ```r ## 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 ``` 去掉那些不明确的细胞类型 ```r 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)](https://pubmed.ncbi.nlm.nih.gov/27667667/)，它本身也包括细胞注释，在数据准备代码中有提及：在`sce.seger$CellType`，这个注释是领域专家设定的 ```r # 以这个新数据集作为测试数据，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 ``` 用这个结果与测试数据本身带的细胞注释结果进行比较： ```r tab <- table(pred.seger$pruned.labels, sce.seger$CellType) library(pheatmap) pheatmap(log2(tab+10), color=colorRampPalette(c("white", "blue"))(101)) ``` 可以看到这两种注释结果比较相近 ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-07-013106.png) ### **2.3 使用marker基因集** > 这部分可以作为2.2的延伸 它基于的想法是：高表达的marker基因与每个细胞的生物学功能相关，因此可以用marker基因作为桥梁，连接测试数据集与参考数据集 基于marker基因的方法会更快捷，同时也不受限于参考数据集是否有专家定义好的细胞类型 **参考数据集 sce.zeisel** [Zeisel et al. (2015)](https://pubmed.ncbi.nlm.nih.gov/25700174/) 目的是将数据集的marker基因于对应的细胞类型联系起来 需要注意的是：下面的`pairwiseWilcox()`代码**需要基于R 4.0**，需要加载归一化后的的数据，前面的数据处理在之前章节都有提及 数据在： ```r 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)](https://pubmed.ncbi.nlm.nih.gov/26727548/) ```r 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**](https://bioconductor.org/packages/3.11/bioc/vignettes/AUCell/inst/doc/AUCell.html) **这个包** 首先拿到参考数据集中的marker基因集`markers.z`，并构建一个`GeneSetCollection`对象 然后对测试数据集中基因表达量进行排名，最后根据排名和之前的基因集进行统计，得到每个marker基因集在每个细胞中的area under the curve (AUC)指标，就能得出哪一个marker基因集在测试数据集的细胞中占比最高 ```r > 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)) ``` ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-07-035814.png) 这样就能清楚看到，**每个细胞对应最高的AUC值属于哪个marker基因集**，因此可以作为这个细胞的label **结果的检查 ——方法一：与参考数据中的cluster类型作比较** ```r new.labels <- colnames(results)[max.col(results)] table(new.labels, sce.tasic$broad_type) ``` ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-07-040813.png) 可以看到一个cluster中的大部分细胞，都被赋予了一个相同的label **结果的检查 ——方法二：计算AUC最大值减去中位数** 先说判断方法：如果注释的很明确，这个差值会较大；较小的差值说明结果不理想；当然我们可以自己定义一个阈值，过滤掉这些无法注释的细胞 下面的例子中，就是假定得到的大多都是注释较好的细胞，然后得到这个阈值 ```r 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 ``` 看看这些被过滤掉的，在整体分布中处于什么位置 ```r 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) ``` ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-07-060634.png) **AUCell 包的好处在于不依赖表达量** 只需要获得基因列表（可以基于文献或其他背景知识获得），例如从MsigDb上下载 ```r 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) ``` 不基于表达量有好有坏，好处就是避免了表达量出现的噪音干扰，处理也更快；坏处是可能会丢掉发现更精细的亚群可能性。 再看一眼另一个数据，不需要运行重点看看差异即可 ```r # 先对测试数据集基因排名 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的子集后，效果变好了；而使用全部的基因集会导致混淆。因此使用哪些基因还是要基于我们对自己测试数据集的认识 ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-07-092916.png) ## 3 基于marker基因的富集分析 > 思路：还是先得到每个cluster的标志性marker基因，然后对这些代表该cluster的基因们进行富集分析，看看它们集中在哪些通路 ### **示例数据** 再以一个新的小鼠乳腺数据 [Bach et al. (2017) ](https://pubmed.ncbi.nlm.nih.gov/29225342/)为例，还是经历了：基因ID注释=》质控=》归一化=》找HVGs=》降维=》聚类=》找marker基因\ [数据可以直接下载](https://share.weiyun.com/HXK6DM5U) ```r #--- 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基因，并且还挑差异显著的基因 ```r chosen <- "2" cur.markers <- markers.mam[[chosen]] is.de <- cur.markers$FDR <= 0.05 summary(is.de) ## Mode FALSE TRUE ## logical 27819 179 ``` 进行GO富集分析 ```r 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)) ``` ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-07-101643.png) 接下来就是考验自己对数据的认知了，假设我们这里对这个数据很熟悉。那么看到这里基因主要参与了lipid synthesis、 cell adhesion 、 tube formation。再结合乳腺的研究背景，推测cluster2可能包含与乳汁产生和分泌相关的管腔上皮细胞，进一步检查一下相关的基因Csn2和Csn3，也是上调的 ```r plotExpression(sce.mam, features=c("Csn2", "Csn3"), x="label", colour_by="label") ``` ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-07-102239.png) 如果想看cluster2中某个感兴趣通路中的基因表现情况 ```r # 比如想看与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 ## ## 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