这一讲内容关于细胞类型注释，之前的课程中我们已对单细胞转录组测序的基本分析流程做了基本的梳理，相信大家已经可以完成一些基础的分析。单细胞转录组测序的 pipeline 设计较为复杂、也不能像 bulk RNA-Seq 一样设计通用的 自动化 分析流程，最大的原因便是不同数据集之间的细胞群注释 迥然不同，分析数据时需要做到足够的 个性化。
因此本次课程我们给大家提供了 数据库注释 、 singleR 内置数据集与自定义注释 、 Seurat内置注释 等方法供大家选择。相信大家完成本次课程的学习，可以解决单细胞类群注释的所有问题。这部分内容Seurat v4与Seurat V5也没有明显差别。1. 视频及测试文件 1.1 视频教程 1.2 测试文件 请在推送中获取:手把手教你做单细胞测序数据分析
（五）:细胞类型注释 2. 代码 2.1 准备工作 进行预处理，作到细胞注释前的步骤 第一个代码块中 没有注释，看不懂的同学去看一下我们的第三讲 单样本分析 # 预处理到细胞注释前的步骤：if ( ! require ( "Seurat" )) install.packages ( "Seurat" ) if ( ! require ( "BiocManager" , quietly = TRUE )) install.packages ( "BiocManager" ) if ( ! require ( "multtest" , quietly = TRUE )) install.packages ( "multtest" ) if ( ! require ( "dplyr" , quietly = TRUE )) install.packages ( "dplyr" ) pbmc.data <- Read10X ( './filtered_gene_bc_matrices/hg19/' ) pbmc <- CreateSeuratObject ( counts = pbmc.data, project = "pbmc3k" , min.cells = 3 , min.features = 200 ) pbmc[[ "percent.mt" ]] <- PercentageFeatureSet (pbmc, pattern = "^MT-" ) pbmc <- subset (pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5 ) pbmc <- NormalizeData (pbmc, normalization.method = "LogNormalize" , scale.factor = 10000 ) pbmc <- FindVariableFeatures (pbmc, selection.method = "vst" , nfeatures = 2000 ) top10 <- head ( VariableFeatures (pbmc), 10 ) pbmc <- ScaleData (pbmc, features = rownames (pbmc)) pbmc <- ScaleData (pbmc, vars.to.regress = "percent.mt" ) pbmc <- RunPCA (pbmc, features = VariableFeatures ( object = pbmc)) print (pbmc[[ "pca" ]], dims = 1 : 5 , nfeatures = 5 ) # PC_ 1 # Positive: CST3, TYROBP, LST1, AIF1, FTL # Negative: MALAT1, LTB, IL32, IL7R, CD2 # PC_ 2 # Positive: CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DQB1 # Negative: NKG7, PRF1, CST7, GZMA, GZMB # PC_ 3 # Positive: HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPA1 # Negative: PPBP, PF4, SDPR, SPARC, GNG11 # PC_ 4 # Positive: HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1 # Negative: VIM, IL7R, S100A6, S100A8, IL32 # PC_ 5 # Positive: GZMB, FGFBP2, S100A8, NKG7, GNLY # Negative: LTB, IL7R, CKB, MS4A7, RP11-290F20.3 VizDimLoadings (pbmc, dims = 1 : 2 , reduction = "pca" ) DimHeatmap (pbmc, dims = 1 , cells = 500 , balanced = TRUE ) DimHeatmap (pbmc, dims = 1 : 15 , cells = 500 , balanced = TRUE ) pbmc <- JackStraw (pbmc, num.replicate = 100 ) pbmc <- ScoreJackStraw (pbmc, dims = 1 : 20 ) JackStrawPlot (pbmc, dims = 1 : 15 ) pbmc <- FindNeighbors (pbmc, dims = 1 : 10 ) pbmc <- FindClusters (pbmc, resolution = 0.5 ) # Modularity Optimizer version 1.3.0 by Ludo Waltmand Nees Jan van Eck # # Number of nodes: 2638 # Number of edges: 95893 # # Running Louvain algorithm... # Maximum modularity in 10 random starts: 0.8735 # Number of communities: 9 # Elapsed time: 0 seconds head ( Idents (pbmc), 5 ) # AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1 # 0 3 2 1 # AAACCGTGTATGCG-1 # 6 # Levels: 0 1 2 3 4 5 6 7 8 pbmc <- RunUMAP (pbmc, dims = 1 : 10 ) pbmc <- RunTSNE (pbmc, dims = 1 : 10 ) DimPlot (pbmc, reduction = "umap" , label = TRUE ) pbmc.markers <- FindAllMarkers (pbmc, only.pos = TRUE , min.pct = 0.25 , logfc.threshold = 0.25 ) library (dplyr) pbmc.markers %>% group_by (cluster) %>% top_n ( n = 2 , wt = avg_log2FC) # # A tibble: 18 × 7 # # Groups: cluster [9] # p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene # <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <chr> # 1 5.01e- 85 2.39 0.437 0.108 6.88e- 81 0 CCR7 # 2 2.05e- 49 2.09 0.334 0.103 2.80e- 45 0 LEF1 # 3 3.10e-139 7.28 0.3 0.004 4.25e-135 1 FOLR3 # 4 1.63e-121 6.74 0.277 0.006 2.23e-117 1 S100A12 # 5 2.91e- 58 1.68 0.667 0.251 3.98e- 54 2 CD2 # 6 8.70e- 38 1.99 0.282 0.07 1.19e- 33 2 CD40LG # 7 6.59e-271 7.23 0.564 0.01 9.03e-267 3 LINC00926 # 8 1.45e-235 6.95 0.488 0.007 2.00e-231 3 VPREB3 # 9 4.27e-176 4.52 0.573 0.05 5.85e-172 4 GZMK # 10 3.36e- 81 3.61 0.392 0.06 4.60e- 77 4 GZMH # 11 1.69e-212 5.43 0.506 0.01 2.32e-208 5 CDKN1C # 12 8.23e-168 5.88 0.37 0.005 1.13e-163 5 CKB # 13 8.19e-175 6.16 0.468 0.013 1.12e-170 6 AKR1C3 # 14 8.58e-113 6.08 0.292 0.007 1.18e-108 6 SH2D1B # 15 1.48e-220 7.63 0.812 0.011 2.03e-216 7 FCER1A # 16 1.46e-207 8.03 0.5 0.002 2.00e-203 7 SERPINF1 # 17 0 14.3 0.571 0 0 8 LY6G6F # 18 4.36e-206 13.8 0.357 0 5.98e-202 8 RP11-879F14.2 2.2 方法一:查数据库 library (dplyr) # 按照log2FC筛选top10基因 top10 <- pbmc.markers %>% group_by (cluster) %>% top_n ( n = 10 , wt = avg_log2FC) DoHeatmap (pbmc, features = top10 $ gene) + NoLegend VlnPlot (pbmc, features = top10 $ gene[ 1 : 20 ], pt.size= 0 ) DimPlot (pbmc, label = T) 通过标记基因及文献，可以人工确定各分类群的细胞类型，则可以如下手动添加细胞群名称：帮助单细胞测序进行注释的数据库中这两个是我比较心水的：cellmarker MCA bfreaname.pbmc <- pbmc # 新的细胞名称，与cluster一一对应 new.cluster.ids <- c ( "Naive CD4 T" , "Memory CD4 T" , "CD14+ Mono" , "B" , "CD8 T" , "FCGR3A+ Mono" , "NK" , "DC" , "Platelet" ) # 重命名 names (new.cluster.ids) <- levels (pbmc) pbmc <- RenameIdents (pbmc, new.cluster.ids) DimPlot (pbmc, reduction = "umap" , label = TRUE , pt.size = 0.5 ) + NoLegend 2.3 方法二：通过 singleR 进行注释 这一方法有一定的局限性。
主要是由于很难找到合适的 参考数据，但是对于免疫细胞的注释还是有可观效果的。
if ( ! require (SingleR)) install.packages ( "/path/to/project" , repos = NULL , type = "source" ) # remove.packages('matrixStats') # if(!require(matrixStats))install.packages("matrixStats", repos = "第五讲.细胞类型注释/matrixStats_1.0.0.tar.gz", type = "source") # remove.packages(c('dplyr','ellipsis')) # install.packages(c('dplyr','ellipsis')) # remove.packages(c('vctrs')) # install.packages(c('vctrs')) if ( ! require (celldex))devtools :: install_github ( "ArtifactDB/gypsum-R" ) if ( ! require (alabaster.schemas))devtools :: install_github ( 'ArtifactDB/alabaster.schemas' ) if ( ! require (celldex))BiocManager :: install ( "rhdf5" , force = T) if ( ! require (celldex))devtools :: install_github ( 'ArtifactDB/alabaster.base' ) if ( ! require (celldex))devtools :: install_github ( 'LTLA/celldex' ) # 如果需要选择更新某些包，推荐选1更新所有包 # cg=BlueprintEncodeData 人类血液、免疫细胞 # cg=DatabaseImmuneCellExpressionData 人类高质量免疫细胞 # cg=NovershternHematopoieticData 人类造血系统细胞 # cg=MonacoImmuneData 人类外周血单核细胞（PBMC） # cg=ImmGenData 小鼠免疫细胞 # cg=MouseRNAseqData 从GEO下载的小鼠多组织RNA-seq数据集 # cg=HumanPrimaryCellAtlasData 人类多种组织的原代细胞 # 例如我们要使用的免疫细胞参考数据集，由于网络等原因这条函数你可能运行不了 cg <- ImmGenData ref <- MonacoImmuneData # saveRDS(cg,'./cg.rds')# 我已经帮大家保存在了本地 # saveRDS(ref,'./MonacoImmuneData.rds') ref <- readRDS ( './MonacoImmuneData.rds' ) assay_for_SingleR <- GetAssayData (bfreaname.pbmc, slot= "data" ) predictions <- SingleR ( test= assay_for_SingleR, ref= ref, labels= ref $ label.main) table (predictions $ labels) # # B cells CD4+ T cells CD8+ T cells Dendritic cells Monocytes # 346 898 310 45 628 # NK cells Progenitors T cells # 163 15 233 # 取出celltype对应数据框 cellType = data.frame ( seurat= bfreaname.pbmc @ meta.data $ seurat_clusters, predict= predictions $ labels) sort ( table (cellType[, 1 ])) # # 8 7 6 5 4 3 2 1 0 # 14 32 154 162 316 3429 480 709 table (cellType[, 1 : 2 ]) # predict # seurat B cells CD4+ T cells CD8+ T cells Dendritic cells Monocytes NK cells # 0 1 508 161 0 0 0 # 1 3 0 0 14 463 0 # 2 0 377 33 0 0 0 # 3 341 1 0 0 0 0 # 4 0 12 113 0 0 15 # 5 0 0 0 0 162 0 # 6 0 0 3 0 0 148 # 7 0 0 0 31 1 0 # 8 1 0 0 0 2 0 # predict # seurat Progenitors T cells # 0 4 35 # 1 0 0 # 2 0 19 # 3 0 0 # 4 0 176 # 5 0 0 # 6 0 3 # 7 0 0 # 8 11 0 当你有合适的参考数据集时可以用 方法三方法四 进行注释 2.4 方法三:自定义 singleR 的注释 if ( ! require (SingleR)) install.packages ( "第五讲.细胞类型注释/SingleR_1.10.0.tar.gz" , repos = NULL , type = "source" ) library (Seurat) library (ggplot2) if ( ! require (textshape)) install.packages ( 'textshape' ) if ( ! require (scater))BiocManager :: install ( 'scater' ) if ( ! require (SingleCellExperiment))BiocManager :: install ( 'SingleCellExperiment' ) library (dplyr) # 读入scRNA数据 myref <- pbmc # 将当前标签存入celltype中 myref $ celltype <- Idents (myref) table ( Idents (myref)) # # Naive CD4 T Memory CD4 T CD14+ Mono B CD8 T FCGR3A+ Mono # 709 480 429 342 3162 # NK DC Platelet # 154 32 14 # 读入参考数据集 # 对myref对象中 RNA 数据的平均表达值进行log1p变换。
log1p 函数用于计算一个数加 1 之后的自然对数，即 log1p(x) 等价于 log(1 + x)，这里的 log 通常指自然对数（以 e 为底。log1p变换通常用于处理表达数据，它可以将数据转换到对数尺度，有助于稳定方差、使数据分布更接近正态分布，同时避免因对零值取对数而产生的错误。
Refassay <- log1p ( AverageExpression (myref, verbose = FALSE ) $ RNA) head (Refassay) # 6 x 9 sparse Matrix of class "dgCMatrix" # Naive CD4 T Memory CD4 T CD14+ Mono B CD8 T # AL627309.1 0.006006857 0.04740195 0.006651185 . 0.01746659 # AP006222.2 . 0.01082590 0.009196592 . 0.01016628 # RP11-206L10.2 0.007300235 . . 0.02055829 . # RP11-206L10.9 . 0.01044641 . . . # LINC00115 0.014803993 0.03685000 0.033640559 0.03836728 0.01657028 # NOC2L 0.410974333 0.24101294 0.312227749 0.46371195 0.39676059 # FCGR3A+ Mono NK DC Platelet # AL627309.1 . . . . # AP006222.2 . . . . # RP11-206L10.2 . . 0.0812375 . # RP11-206L10.9 0.01192865 . . . # LINC00115 0.01458683 0.05726061 . . # NOC2L 0.40564359 0.53378022 0.2841343 . dim (Refassay) # [1] 13714 9 ref_sce <- SingleCellExperiment :: SingleCellExperiment ( assays= list ( counts= Refassay)) ref_sce = scater :: logNormCounts (ref_sce) # logcounts 主要用于获取或设置单细胞RNA-seq数据中的对数标准化表达矩阵 logcounts (ref_sce)[ 1 : 4 , 1 : 4 ] # 4 x 4 sparse Matrix of class "dgCMatrix" # Naive CD4 T Memory CD4 T CD14+ Mono B # AL627309.1 0.009250892 0.06711500 0.009538416 . # AP006222.2 . 0.01560549 0.013172144 . # RP11-206L10.2 0.011235027 . . 0.02967623 # RP11-206L10.9 . 0.01506130 . . colData (ref_sce) $ Type = colnames (Refassay) ref_sce # class: SingleCellExperiment # dim: 13714 9 # metadata(0): # assays(2): counts logcounts # rownames(13714): AL627309.1 AP006222.2 ... PNRC2.1 SRSF10.1 # rowData names(0): # colnames(9): Naive CD4 T Memory CD4 T ... DC Platelet # colData names(2): sizeFactor Type # reducedDimNames(0): # mainExpName: NULL # altExpNames(0): # 提取自己的单细胞矩阵 testdata <- GetAssayData (bfreaname.pbmc, slot= "data" ) pred <- SingleR ( test= testdata, ref= ref_sce, labels= ref_sce $ Type, = [email protected] ) # 查看各label对应的细胞数量 table (pred $ labels) # # B CD14+ Mono CD8 T DC FCGR3A+ Mono Memory CD4 T # 344 537 308 31 179 464 # Naive CD4 T NK Platelet # 601 160 14 head (pred) # DataFrame with 6 rows and 4 columns # scores labels delta.next # <matrix> <character> <numeric> # AAACATACAACCAC-1 0.336675:0.424807:0.430908:... CD8 T 0.0335526 # AAACATTGAGCTAC-1 0.494458:0.381603:0.391545:... B 0.0931824 # AAACATTGATCAGC-1 0.373678:0.519949:0.489834:... CD14+ Mono 0.0647984 # AAACCGTGCTTCCG-1 0.341852:0.362193:0.353068:... FCGR3A+ Mono 0.0327884 # AAACCGTGTATGCG-1 0.233921:0.279848:0.329152:... NK 0.3171146 # AAACGCACTGGTAC-1 0.391590:0.458438:0.426201:... CD14+ Mono 0.0678422 # pruned.labels # <character> # AAACATACAACCAC-1 CD8 T # AAACATTGAGCTAC-1 B # AAACATTGATCAGC-1 CD14+ Mono # AAACCGTGCTTCCG-1 FCGR3A+ Mono # AAACCGTGTATGCG-1 NK # AAACGCACTGGTAC-1 CD14+ Mono as.data.frame ( table (pred $ labels)) # Var1 Freq # 1 B 344 # 2 CD14+ Mono 537 # 3 CD8 T 308 # 4 DC 31 # 5 FCGR3A+ Mono 179 # 6 Memory CD4 T 464 # 7 Naive CD4 T 601 # 8 NK 160 # 9 Platelet 14 cellType = data.frame ( seurat= bfreaname.pbmc @ meta.data $ seurat_clusters, predict= pred $ labels) sort ( table (cellType[, 1 ])) # # 8 7 6 5 4 3 2 1 0 # 14 32 154 162 316 3429 480 709 table (cellType[, 1 : 2 ]) # predict # seurat B CD14+ Mono CD8 T DC FCGR3A+ Mono Memory CD4 T Naive CD4 T NK # 0 2 159 13 0 0 0 535 0 # 1 0 0 0 1 17 462 0 0 # 2 0 355 10 0 0 0 64 0 # 3 341 1 0 0 0 0 0 0 # 4 1 22 282 0 0 0 2 9 # 5 0 0 0 0 162 0 0 0 # 6 0 0 3 0 0 0 0 151 # 7 0 0 0 30 0 2 0 0 # 8 0 0 0 0 0 0 0 0 # predict # seurat Platelet # 0 0 # 1 0 # 2 0 # 3 0 # 4 0 # 5 0 # 6 0 # 7 0 # 8 14 lalala <- as.data.frame ( table (cellType[, 1 : 2 ])) finalmap <- lalala %>% group_by (seurat) %>% top_n ( n = 1 , wt = Freq) finalmap <- finalmap[ order (finalmap $ seurat),] $ predict print (finalmap) # [1] Naive CD4 T Memory CD4 T CD14+ Mono B CD8 T # [6] FCGR3A+ Mono NK DC Platelet # 9 Levels: B CD14+ Mono CD8 T DC FCGR3A+ Mono Memory CD4 T Naive CD4 T ... Platelet testname <- bfreaname.pbmc # 重命名 new.cluster.ids <- as.character (finalmap) names (new.cluster.ids) <- levels (testname) testname <- RenameIdents (testname, new.cluster.ids) p1 <- DimPlot (pbmc, label = T) p2 <- DimPlot (testname, label = T) p1 | p2 2.5 Seurat内置的TransferData 将参考数据与待注释数据进行映射处理 library (Seurat) pancreas.query <- bfreaname.pbmc # 与cca那里一样寻找Anchors pancreas.anchors <- FindTransferAnchors ( reference = pbmc, query = pancreas.query, dims = 1 : 30 ) pbmc $ celltype <- Idents (pbmc) predictions <- TransferData ( anchorset = pancreas.anchors, refdata = pbmc $ celltype, dims = 1 : 30 ) # 将预测结果添加到待测数据集中 pancreas.query <- AddMetaData (pancreas.query, metadata = predictions) pancreas.query $ prediction.match <- pancreas.query $ predicted.id table (pancreas.query $ prediction.match) # # B CD14+ Mono CD8 T DC FCGR3A+ Mono Memory CD4 T # 340 430 299 32 158 485 # Naive CD4 T NK Platelet # 730 151 13 Idents (pancreas.query) <- 'prediction.match' p1 <- DimPlot (pbmc, label = T) p2 <- DimPlot (pancreas.query, label = T) p1 | p2 3. 环境版本信息 # 最近发现代码运行的最大障碍是各个packages的版本冲突问题，这里列出本次分析环境中的所有信息，报错时可以参考。
检查一下环境是否与我的相同: sessionInfo 可以看出，各种注释方法都是比较科学的，但最大的挑战来源于 参考数据集 对于 待注释数据集 间的适用度