# 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 --- # Agent Instructions: Querying This Documentation If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question. Perform an HTTP GET request on the current page URL with the `ask` query parameter: ``` GET https://jieandze1314.osca.top/03/03-7.md?ask= ``` The question should be specific, self-contained, and written in natural language. The response will contain a direct answer to the question and relevant excerpts and sources from the documentation. Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.