热门标签
热门文章
- 1NLP自然语言处理（二）—— 语料及词性标注 & 分词 & TFIDF_语料标注
- 2分布式系统的跟踪系统：Dapper设计给我们的启示_谷歌哪篇论文提到了 dapper
- 3Go 首次冲进前八！曾两次夺得年度编程语言，也曾跌至百名开外 | TIOBE 2 月榜单发布_golang市场份额
- 4百度AI攻略：图像对比度增强_ai contrast
- 5java 服务证书问题-生成jssecacerts CA证书
- 6腾讯云4核8G服务器性能怎么样？4核8G12M服务器能干嘛？_轻量应用服务器4核8g有gpu
- 7Pascal VOC 数据集介绍_pascal voc2007/2012数据集
- 8什么是AI大模型?大模型应用以及未来发展趋势
- 9预测房价：人工智能回归问题_“橙现智能”房价预测模型说明
- 10【AI大模型应用开发】【LangChain系列】5. 实战LangChain的智能体Agents模块_hwchase17/structured-chat-agent
当前位置: article > 正文
手把手教你做单细胞测序数据分析（五）——细胞类型注释_单细胞测序生信分析代码
作者：Monodyee | 2024-03-31 16:21:57
赞
踩
单细胞测序生信分析代码
前面的课程中基本已带领大家将单细胞测序预处理部分打通了，这里给大家介绍的是单细胞测序中最让人头疼的细胞类型注释。第一个代码块中没有注释，看不懂的同学去看一下我们的第三讲 单样本分析，测试数据也与那里的相同。
视频教程：
手把手教你做单细胞测序数据分析（五）——细胞类型注释
（B站同步播出，先看一遍视频再跟着代码一起操作，建议每个视频至少看三遍）
准备工作
先进行预处理，作到细胞注释前的步骤 if(T){rm(list = ls()) 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") download.file('https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz','pbmc3k_filtered_gene_bc_matrices.tar.gz') library(R.utils) gunzip('pbmc3k_filtered_gene_bc_matrices.tar.gz') untar('pbmc3k_filtered_gene_bc_matrices.tar') 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) 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) head(Idents(pbmc), 5) 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) } ## Loading required package: Seurat ## Warning: package 'Seurat' was built under R version 4.0.5 ## Attaching SeuratObject ## Warning: package 'BiocManager' was built under R version 4.0.5 ## Bioconductor version '3.12' is out-of-date; the current release version '3.13' ## is available with R version '4.1'; see https://bioconductor.org/install ## Warning: package 'BiocGenerics' was built under R version 4.0.5 ## ## Attaching package: 'BiocGenerics' ## The following objects are masked from 'package:parallel': ## ## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ, ## clusterExport, clusterMap, parApply, parCapply, parLapply, ## parLapplyLB, parRapply, parSapply, parSapplyLB ## The following objects are masked from 'package:stats': ## ## IQR, mad, sd, var, xtabs ## The following objects are masked from 'package:base': ## ## anyDuplicated, append, as.data.frame, basename, cbind, colnames, ## dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep, ## grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget, ## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank, ## rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply, ## union, unique, unsplit, which.max, which.min ## Welcome to Bioconductor ## ## Vignettes contain introductory material; view with ## 'browseVignettes()'. To cite Bioconductor, see ## 'citation("Biobase")', and for packages 'citation("pkgname")'. ## Warning: package 'dplyr' was built under R version 4.0.5 ## ## Attaching package: 'dplyr' ## The following object is masked from 'package:Biobase': ## ## combine ## The following objects are masked from 'package:BiocGenerics': ## ## combine, intersect, setdiff, union ## The following objects are masked from 'package:stats': ## ## filter, lag ## The following objects are masked from 'package:base': ## ## intersect, setdiff, setequal, union ## Loading required package: R.oo ## Loading required package: R.methodsS3 ## R.methodsS3 v1.8.1 (2020-08-26 16:20:06 UTC) successfully loaded. See ?R.methodsS3 for help. ## R.oo v1.24.0 (2020-08-26 16:11:58 UTC) successfully loaded. See ?R.oo for help. ## ## Attaching package: 'R.oo' ## The following object is masked from 'package:R.methodsS3': ## ## throw ## The following objects are masked from 'package:methods': ## ## getClasses, getMethods ## The following objects are masked from 'package:base': ## ## attach, detach, load, save ## R.utils v2.10.1 (2020-08-26 22:50:31 UTC) successfully loaded. See ?R.utils for help. ## ## Attaching package: 'R.utils' ## The following object is masked from 'package:utils': ## ## timestamp ## The following objects are masked from 'package:base': ## ## cat, commandArgs, getOption, inherits, isOpen, nullfile, parse, ## warnings ## Warning: Feature names cannot have underscores ('_'), replacing with dashes ## ('-') ## Centering and scaling data matrix ## Regressing out percent.mt ## Centering and scaling data matrix ## PC_ 1 ## Positive: CST3, TYROBP, LST1, AIF1, FTL, FTH1, LYZ, FCN1, S100A9, TYMP ## FCER1G, CFD, LGALS1, LGALS2, SERPINA1, S100A8, CTSS, IFITM3, SPI1, CFP ## PSAP, IFI30, COTL1, SAT1, S100A11, NPC2, GRN, LGALS3, GSTP1, PYCARD ## Negative: MALAT1, LTB, IL32, IL7R, CD2, B2M, ACAP1, CTSW, STK17A, CD27 ## CD247, CCL5, GIMAP5, GZMA, AQP3, CST7, TRAF3IP3, SELL, GZMK, HOPX ## MAL, MYC, ITM2A, ETS1, LYAR, GIMAP7, KLRG1, NKG7, ZAP70, BEX2 ## PC_ 2 ## Positive: CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DQB1, HLA-DRA, LINC00926, CD79B, HLA-DRB1, CD74 ## HLA-DMA, HLA-DPB1, HLA-DQA2, CD37, HLA-DRB5, HLA-DMB, HLA-DPA1, FCRLA, HVCN1, LTB ## BLNK, P2RX5, IGLL5, IRF8, SWAP70, ARHGAP24, FCGR2B, SMIM14, PPP1R14A, C16orf74 ## Negative: NKG7, PRF1, CST7, GZMA, GZMB, FGFBP2, CTSW, GNLY, B2M, SPON2 ## CCL4, GZMH, FCGR3A, CCL5, CD247, XCL2, CLIC3, AKR1C3, SRGN, HOPX ## TTC38, CTSC, APMAP, S100A4, IGFBP7, ANXA1, ID2, IL32, XCL1, RHOC ## PC_ 3 ## Positive: HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPA1, HLA-DPB1, CD74, MS4A1, HLA-DRB1, HLA-DRA ## HLA-DRB5, HLA-DQA2, TCL1A, LINC00926, HLA-DMB, HLA-DMA, CD37, HVCN1, FCRLA, IRF8 ## PLAC8, BLNK, MALAT1, SMIM14, PLD4, IGLL5, SWAP70, P2RX5, LAT2, FCGR3A ## Negative: PPBP, PF4, SDPR, SPARC, GNG11, NRGN, GP9, RGS18, TUBB1, CLU ## HIST1H2AC, AP001189.4, ITGA2B, CD9, TMEM40, PTCRA, CA2, ACRBP, MMD, TREML1 ## NGFRAP1, F13A1, SEPT5, RUFY1, TSC22D1, CMTM5, MPP1, MYL9, RP11-367G6.3, GP1BA ## PC_ 4 ## Positive: HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1, CD74, HLA-DPB1, HIST1H2AC, HLA-DPA1, HLA-DRB1 ## TCL1A, PF4, HLA-DQA2, SDPR, HLA-DRA, LINC00926, PPBP, GNG11, HLA-DRB5, SPARC ## GP9, PTCRA, CA2, AP001189.4, CD9, NRGN, RGS18, GZMB, CLU, TUBB1 ## Negative: VIM, IL7R, S100A6, S100A8, IL32, S100A4, GIMAP7, S100A10, S100A9, MAL ## AQP3, CD14, CD2, LGALS2, FYB, GIMAP4, ANXA1, RBP7, CD27, FCN1 ## LYZ, S100A12, MS4A6A, GIMAP5, S100A11, FOLR3, TRABD2A, AIF1, IL8, TMSB4X ## PC_ 5 ## Positive: GZMB, FGFBP2, S100A8, NKG7, GNLY, CCL4, PRF1, CST7, SPON2, GZMA ## GZMH, LGALS2, S100A9, CCL3, XCL2, CD14, CLIC3, CTSW, MS4A6A, GSTP1 ## S100A12, RBP7, IGFBP7, FOLR3, AKR1C3, TYROBP, CCL5, TTC38, XCL1, APMAP ## Negative: LTB, IL7R, CKB, MS4A7, RP11-290F20.3, AQP3, SIGLEC10, VIM, CYTIP, HMOX1 ## LILRB2, PTGES3, HN1, CD2, FAM110A, CD27, ANXA5, CTD-2006K23.1, MAL, VMO1 ## CORO1B, TUBA1B, LILRA3, GDI2, TRADD, ATP1A1, IL32, ABRACL, CCDC109B, PPA1 ## 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
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 79
- 80
- 81
- 82
- 83
- 84
- 85
- 86
- 87
- 88
- 89
- 90
- 91
- 92
- 93
- 94
- 95
- 96
- 97
- 98
- 99
- 100
- 101
- 102
- 103
- 104
- 105
- 106
- 107
- 108
- 109
- 110
- 111
- 112
- 113
- 114
- 115
- 116
- 117
- 118
- 119
- 120
- 121
- 122
- 123
- 124
- 125
- 126
- 127
- 128
- 129
- 130
- 131
- 132
- 133
- 134
- 135
- 136
- 137
- 138
- 139
- 140
- 141
- 142
- 143
- 144
- 145
- 146
- 147
- 148
- 149
- 150
- 151
- 152
- 153
- 154
- 155
- 156
- 157
- 158
- 159
- 160
- 161
- 162
- 163
## Computing nearest neighbor graph ## Computing SNN ## Modularity Optimizer version 1.3.0 by Ludo Waltman and 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 ## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric ## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation' ## This message will be shown once per session ## 09:15:44 UMAP embedding parameters a = 0.9922 b = 1.112 ## 09:15:44 Read 2638 rows and found 10 numeric columns ## 09:15:44 Using Annoy for neighbor search, n_neighbors = 30 ## 09:15:44 Building Annoy index with metric = cosine, n_trees = 50 ## 0% 10 20 30 40 50 60 70 80 90 100% ## [----|----|----|----|----|----|----|----|----|----| ## **************************************************| ## 09:15:44 Writing NN index file to temp file C:\Users\53513\AppData\Local\Temp\Rtmpo9Xcsm\file56685416913 ## 09:15:44 Searching Annoy index using 1 thread, search_k = 3000 ## 09:15:45 Annoy recall = 100% ## 09:15:45 Commencing smooth kNN distance calibration using 1 thread ## 09:15:46 Initializing from normalized Laplacian + noise ## 09:15:46 Commencing optimization for 500 epochs, with 106338 positive edges ## 09:15:52 Optimization finished ## Calculating cluster 0 ## Calculating cluster 1 ## Calculating cluster 2 ## Calculating cluster 3 ## Calculating cluster 4 ## Calculating cluster 5 ## Calculating cluster 6 ## Calculating cluster 7 ## Calculating cluster 8
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
## Registered S3 method overwritten by 'cli': ## method from ## print.boxx spatstat.geom ## # A tibble: 18 x 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 1.88e-117 1.08 0.913 0.588 2.57e-113 0 LDHB ## 2 5.01e- 85 1.34 0.437 0.108 6.88e- 81 0 CCR7 ## 3 0 5.57 0.996 0.215 0 1 S100A9 ## 4 0 5.48 0.975 0.121 0 1 S100A8 ## 5 1.93e- 80 1.27 0.981 0.65 2.65e- 76 2 LTB ## 6 2.91e- 58 1.27 0.667 0.251 3.98e- 54 2 CD2 ## 7 0 4.31 0.939 0.042 0 3 CD79A ## 8 1.06e-269 3.59 0.623 0.022 1.45e-265 3 TCL1A ## 9 3.60e-221 3.21 0.984 0.226 4.93e-217 4 CCL5 ## 10 4.27e-176 3.01 0.573 0.05 5.85e-172 4 GZMK ## 11 3.51e-184 3.31 0.975 0.134 4.82e-180 5 FCGR3A ## 12 2.03e-125 3.09 1 0.315 2.78e-121 5 LST1 ## 13 3.17e-267 4.83 0.961 0.068 4.35e-263 6 GZMB ## 14 1.04e-189 5.28 0.961 0.132 1.43e-185 6 GNLY ## 15 1.48e-220 3.87 0.812 0.011 2.03e-216 7 FCER1A ## 16 1.67e- 21 2.87 1 0.513 2.28e- 17 7 HLA-DPB1 ## 17 7.73e-200 7.24 1 0.01 1.06e-195 8 PF4 ## 18 3.68e-110 8.58 1 0.024 5.05e-106 8 PPBP
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
方法一：查数据库
library(dplyr)
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()#展示前10个标记基因的热图
## Warning in DoHeatmap(pbmc, features = top10$gene): The following features were
## omitted as they were not found in the scale.data slot for the RNA assay: CD8A,
## VPREB3, CD40LG, PIK3IP1, PRKCQ-AS1, NOSIP, LEF1, CD3E, CD3D, CCR7, LDHB, RPS3A
- 1
- 2
- 3
- 4
- 5
- 6
VlnPlot(pbmc, features = top10$gene[1:20],pt.size=0)
- 1
DimPlot(pbmc,label = T)
- 1
#通过标记基因及文献，可以人工确定各分类群的细胞类型，则可以如下手动添加细胞群名称
bfreaname.pbmc <- pbmc
new.cluster.ids <- c("Naive CD4 T", "Memory CD4 T", "CD14+ Mono", "B", "CD8 T", "FCGR3A+ Mono",
"NK", "DC", "Platelet")
#帮助单细胞测序进行注释的数据库：
#http://mp.weixin.qq.com/s?__biz=MzI5MTcwNjA4NQ==&mid=2247502903&idx=2&sn=fd21e6e111f57a4a2b6c987e391068fd&chksm=ec0e09bddb7980abf038f62d03d3beea6249753c8fba69b69f399de9854fc208ca863ca5bc23&mpshare=1&scene=24&srcid=1110SJhxDL8hmNB5BThrgOS9&sharer_sharetime=1604979334616&sharer_shareid=853c5fb0f1636baa0a65973e8b5db684#rd
#cellmarker: http://biocc.hrbmu.edu.cn/CellMarker/index.jsp
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
方法二：通过singleR
这一方法很鸡肋，主要是由于很难找到合适的参考数据，但是对于免疫细胞的注释还是有可观效果的
if(!require(SingleR))BiocManager::install(SingleR) ## Loading required package: SingleR ## Loading required package: SummarizedExperiment ## Loading required package: MatrixGenerics ## Loading required package: matrixStats ## Warning: package 'matrixStats' was built under R version 4.0.5 ## ## Attaching package: 'matrixStats' ## The following object is masked from 'package:dplyr': ## ## count ## The following objects are masked from 'package:Biobase': ## ## anyMissing, rowMedians ## ## Attaching package: 'MatrixGenerics' ## The following objects are masked from 'package:matrixStats': ## ## colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse, ## colCounts, colCummaxs, colCummins, colCumprods, colCumsums, ## colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs, ## colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats, ## colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds, ## colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads, ## colWeightedMeans, colWeightedMedians, colWeightedSds, ## colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet, ## rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods, ## rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps, ## rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins, ## rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks, ## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars, ## rowWeightedMads, rowWeightedMeans, rowWeightedMedians, ## rowWeightedSds, rowWeightedVars ## The following object is masked from 'package:Biobase': ## ## rowMedians ## Loading required package: GenomicRanges ## Loading required package: stats4 ## Loading required package: S4Vectors ## ## Attaching package: 'S4Vectors' ## The following objects are masked from 'package:dplyr': ## ## first, rename ## The following object is masked from 'package:base': ## ## expand.grid ## Loading required package: IRanges ## ## Attaching package: 'IRanges' ## The following object is masked from 'package:R.oo': ## ## trim ## The following objects are masked from 'package:dplyr': ## ## collapse, desc, slice ## The following object is masked from 'package:grDevices': ## ## windows ## Loading required package: GenomeInfoDb ## Warning: package 'GenomeInfoDb' was built under R version 4.0.5 ## ## Attaching package: 'SummarizedExperiment' ## The following object is masked from 'package:SeuratObject': ## ## Assays ## The following object is masked from 'package:Seurat': ## ## Assays # remove.packages('matrixStats')#移除后可能需要重新打开 if(!require(matrixStats))BiocManager::install('matrixStats') # remove.packages(c('dplyr','ellipsis')) # install.packages(c('dplyr','ellipsis')) # remove.packages(c('vctrs')) # install.packages(c('vctrs')) if(!require(celldex))BiocManager::install('celldex')#有一些版本冲突，需要重新安装一些包。 ## Loading required package: celldex ## ## Attaching package: 'celldex' ## The following objects are masked from 'package:SingleR': ## ## BlueprintEncodeData, DatabaseImmuneCellExpressionData, ## HumanPrimaryCellAtlasData, ImmGenData, MonacoImmuneData, ## MouseRNAseqData, NovershternHematopoieticData ###内置数据集： # cg=BlueprintEncodeData() # cg=DatabaseImmuneCellExpressionData() # cg=NovershternHematopoieticData() # cg=MonacoImmuneData() # cg=ImmGenData() # cg=MouseRNAseqData() # cg=HumanPrimaryCellAtlasData() cg<-ImmGenData()#选取我们要使用的免疫细胞参考数据集 ## snapshotDate(): 2020-10-27 ## see ?celldex and browseVignettes('celldex') for documentation ## Could not check id: EH3494 for updates. ## Using previously cached version. ## loading from cache ## Could not check id: EH3494 for updates. ## Using previously cached version. ## see ?celldex and browseVignettes('celldex') for documentation ## Could not check id: EH3495 for updates. ## Using previously cached version. ## loading from cache ## Could not check id: EH3495 for updates. ## Using previously cached version. assay_for_SingleR <- GetAssayData(bfreaname.pbmc, slot="data")#取出样本中的表达序列 predictions <- SingleR(test=assay_for_SingleR, ref=cg, labels=cg$label.main) #以kidney中提取的阵列为输入数据，以小鼠的阵列作为参考，predict细胞类型 table(predictions$labels)#看看都注释到了哪些细胞 ## ## B cells Endothelial cells Eosinophils Mast cells ## 2555 19 18 25 ## NK cells Stromal cells ## 11 10 cellType=data.frame(seurat=bfreaname.pbmc@meta.data$seurat_clusters, predict=predictions$labels)#得到seurat中编号与预测标签之间的关系 sort(table(cellType[,1])) ## ## 8 7 6 5 4 3 2 1 0 ## 14 32 154 162 316 342 429 480 709 table(cellType[,1:2])#访问celltyple的2~3列 ## predict ## seurat B cells Endothelial cells Eosinophils Mast cells NK cells Stromal cells ## 0 696 5 4 0 3 1 ## 1 458 2 6 12 2 0 ## 2 415 4 6 0 2 2 ## 3 342 0 0 0 0 0 ## 4 311 4 0 0 1 0 ## 5 137 3 1 13 1 7 ## 6 151 1 1 0 1 0 ## 7 31 0 0 0 1 0 ## 8 14 0 0 0 0 0 #可以看出 singleR如果没有合适的数据集，得到的结果有多拉跨了吧 #这里就没必要再往下重命名了
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 79
- 80
- 81
- 82
- 83
- 84
- 85
- 86
- 87
- 88
- 89
- 90
- 91
- 92
- 93
- 94
- 95
- 96
- 97
- 98
- 99
- 100
- 101
- 102
- 103
- 104
- 105
- 106
- 107
- 108
- 109
- 110
- 111
- 112
- 113
- 114
- 115
- 116
- 117
- 118
- 119
- 120
- 121
- 122
- 123
- 124
- 125
- 126
- 127
- 128
- 129
- 130
- 131
- 132
- 133
- 134
- 135
- 136
- 137
你有合适的参考数据集时可以用方法三、方法四进行注释
方法三：自定义singleR注释
#我们走个极端，拿它自己作为自己的参考数据集，看看注释的准不准 ##########利用singleR构建自己的数据作为参考数据集######## library(SingleR) library(Seurat) library(ggplot2) ## Warning: package 'ggplot2' was built under R version 4.0.5 if(!require(textshape))install.packages('textshape') ## Loading required package: textshape ## Warning: package 'textshape' was built under R version 4.0.5 ## ## Attaching package: 'textshape' ## The following object is masked from 'package:dplyr': ## ## combine ## The following object is masked from 'package:Biobase': ## ## combine ## The following object is masked from 'package:BiocGenerics': ## ## combine if(!require(scater))BiocManager::install('scater') ## Loading required package: scater ## Warning: package 'scater' was built under R version 4.0.4 ## Loading required package: SingleCellExperiment if(!require(SingleCellExperiment))BiocManager::install('SingleCellExperiment') library(dplyr) # 读入scRNA数据 ------- myref <- pbmc##这里为了检测，我们将参考数据集与目标数据集用同一个数据进行测试 myref$celltype <- Idents(myref) table(Idents(myref)) ## ## Naive CD4 T Memory CD4 T CD14+ Mono B CD8 T FCGR3A+ Mono ## 709 480 429 342 316 162 ## NK DC Platelet ## 154 32 14 # 读入参考数据集 ------- Refassay <- log1p(AverageExpression(myref, verbose = FALSE)$RNA)#求 #Ref <- textshape::column_to_rownames(Ref, loc = 1)#另一种得到参考矩阵的办法 head(Refassay)#看看表达矩阵长啥样 ## Naive CD4 T Memory CD4 T CD14+ Mono B CD8 T ## AL627309.1 0.006006857 0.04740195 0.006651185 0.00000000 0.01746659 ## AP006222.2 0.000000000 0.01082590 0.009196592 0.00000000 0.01016628 ## RP11-206L10.2 0.007300235 0.00000000 0.000000000 0.02055829 0.00000000 ## RP11-206L10.9 0.000000000 0.01044641 0.000000000 0.00000000 0.00000000 ## 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 0.00000000 0.00000000 0.0000000 0 ## AP006222.2 0.00000000 0.00000000 0.0000000 0 ## RP11-206L10.2 0.00000000 0.00000000 0.0812375 0 ## RP11-206L10.9 0.01192865 0.00000000 0.0000000 0 ## LINC00115 0.01458683 0.05726061 0.0000000 0 ## NOC2L 0.40564359 0.53378022 0.2841343 0 ref_sce <- SingleCellExperiment::SingleCellExperiment(assays=list(counts=Refassay)) #参考数据集需要构建成一个SingleCellExperiment对象 ref_sce=scater::logNormCounts(ref_sce) logcounts(ref_sce)[1:4,1:4] ## Naive CD4 T Memory CD4 T CD14+ Mono B ## AL627309.1 0.009250892 0.06711500 0.009538416 0.00000000 ## AP006222.2 0.000000000 0.01560549 0.013172144 0.00000000 ## RP11-206L10.2 0.011235027 0.00000000 0.000000000 0.02967623 ## RP11-206L10.9 0.000000000 0.01506130 0.000000000 0.00000000 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): ## altExpNames(0): ###提取自己的单细胞矩阵########## testdata <- GetAssayData(bfreaname.pbmc, slot="data") pred <- SingleR(test=testdata, ref=ref_sce, labels=ref_sce$Type, #clusters = scRNA@active.ident ) 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 5 columns ## scores first.labels tuning.scores ## <matrix> <character> <DataFrame> ## AAACATACAACCAC-1 0.336675:0.424807:0.430908:... CD8 T 0.321456:0.2879031 ## AAACATTGAGCTAC-1 0.494458:0.381603:0.391545:... B 0.494458:0.4012752 ## AAACATTGATCAGC-1 0.373678:0.519949:0.489834:... CD14+ Mono 0.411459:0.3466608 ## AAACCGTGCTTCCG-1 0.341852:0.362193:0.353068:... Memory CD4 T 0.422782:0.3899934 ## AAACCGTGTATGCG-1 0.233921:0.279848:0.329152:... NK 0.384988:0.0678734 ## AAACGCACTGGTAC-1 0.391590:0.458438:0.426201:... CD14+ Mono 0.335794:0.2679515 ## labels pruned.labels ## <character> <character> ## AAACATACAACCAC-1 CD8 T CD8 T ## AAACATTGAGCTAC-1 B B ## AAACATTGATCAGC-1 CD14+ Mono CD14+ Mono ## AAACCGTGCTTCCG-1 FCGR3A+ Mono FCGR3A+ Mono ## AAACCGTGTATGCG-1 NK NK ## AAACGCACTGGTAC-1 CD14+ Mono 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 #pred@listData[["scores"]] #预测评分，想看看结构的可以自己看看 #同上，我们找一下seurat中类群与注释结果直接的关系 cellType=data.frame(seurat=bfreaname.pbmc@meta.data$seurat_clusters, predict=pred$labels)#得到seurat中编号与预测标签之间的关系 sort(table(cellType[,1])) ## ## 8 7 6 5 4 3 2 1 0 ## 14 32 154 162 316 342 429 480 709 table(cellType[,1:2])#访问celltyple的2~3列 ## 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)#找出每种seurat_cluster注释比例最高的对应类型 finalmap <-finalmap[order(finalmap$seurat),]$predict#找到seurat中0：8的对应预测细胞类型 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#完美，无差别注释，当然了，我们这个参考数据用的比较极端
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 79
- 80
- 81
- 82
- 83
- 84
- 85
- 86
- 87
- 88
- 89
- 90
- 91
- 92
- 93
- 94
- 95
- 96
- 97
- 98
- 99
- 100
- 101
- 102
- 103
- 104
- 105
- 106
- 107
- 108
- 109
- 110
- 111
- 112
- 113
- 114
- 115
- 116
- 117
- 118
- 119
- 120
- 121
- 122
- 123
- 124
- 125
- 126
- 127
- 128
- 129
- 130
- 131
- 132
- 133
- 134
- 135
- 136
- 137
- 138
- 139
- 140
- 141
- 142
- 143
- 144
- 145
- 146
- 147
- 148
- 149
- 150
- 151
- 152
- 153
- 154
- 155
- 156
- 157
- 158
- 159
- 160
- 161
- 162
- 163
方法四、Seurat内置的投影
######利用seurat内置的原先用于细胞整合的功能，将参考数据与待注释数据进行映射处理 library(Seurat) pancreas.query <- bfreaname.pbmc#待注释数据 pancreas.anchors <- FindTransferAnchors(reference = pbmc, query = pancreas.query, dims = 1:30) ## Performing PCA on the provided reference using 2000 features as input. ## Projecting cell embeddings ## Finding neighborhoods ## Finding anchors ## Found 7206 anchors ## Filtering anchors ## Retained 5268 anchors pbmc$celltype <- Idents(pbmc) predictions <- TransferData(anchorset = pancreas.anchors, refdata = pbmc$celltype, dims = 1:30) ## Finding integration vectors ## Finding integration vector weights ## Predicting cell labels 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 431 294 32 158 485 ## Naive CD4 T NK Platelet ## 732 153 13 Idents(pancreas.query)<- 'prediction.match' p1 <- DimPlot(pbmc,label = T) p2 <- DimPlot(pancreas.query,label = T)#比较一下测试数据与参考数据集之间有没有偏差 p1|p2#完美，无差别注释，当然了，我们这个参考数据用的比较极端
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
版本问题
#### 最近发现代码运行的最大障碍是各个packages的版本冲突问题，这里列出本次分析环境中的所有信息，报错时可以参考 sessionInfo() ## R version 4.0.3 (2020-10-10) ## Platform: x86_64-w64-mingw32/x64 (64-bit) ## Running under: Windows 10 x64 (build 19042) ## ## Matrix products: default ## ## locale: ## [1] LC_COLLATE=Chinese (Simplified)_China.936 ## [2] LC_CTYPE=Chinese (Simplified)_China.936 ## [3] LC_MONETARY=Chinese (Simplified)_China.936 ## [4] LC_NUMERIC=C ## [5] LC_TIME=Chinese (Simplified)_China.936 ## ## attached base packages: ## [1] stats4 parallel stats graphics grDevices utils datasets ## [8] methods base ## ## other attached packages: ## [1] scater_1.18.6 SingleCellExperiment_1.12.0 ## [3] textshape_1.7.3 ggplot2_3.3.5 ## [5] celldex_1.0.0 SingleR_1.4.1 ## [7] SummarizedExperiment_1.20.0 GenomicRanges_1.42.0 ## [9] GenomeInfoDb_1.26.7 IRanges_2.24.1 ## [11] S4Vectors_0.28.1 MatrixGenerics_1.2.1 ## [13] matrixStats_0.60.1 R.utils_2.10.1 ## [15] R.oo_1.24.0 R.methodsS3_1.8.1 ## [17] dplyr_1.0.7 multtest_2.46.0 ## [19] Biobase_2.50.0 BiocGenerics_0.36.1 ## [21] BiocManager_1.30.16 SeuratObject_4.0.2 ## [23] Seurat_4.0.2 ## ## loaded via a namespace (and not attached): ## [1] utf8_1.2.2 reticulate_1.20 ## [3] tidyselect_1.1.1 RSQLite_2.2.8 ## [5] AnnotationDbi_1.52.0 htmlwidgets_1.5.4 ## [7] grid_4.0.3 BiocParallel_1.24.1 ## [9] Rtsne_0.15 munsell_0.5.0 ## [11] codetools_0.2-16 ica_1.0-2 ## [13] future_1.22.1 miniUI_0.1.1.1 ## [15] withr_2.4.2 colorspace_2.0-2 ## [17] highr_0.9 knitr_1.34 ## [19] rstudioapi_0.13 ROCR_1.0-11 ## [21] tensor_1.5 listenv_0.8.0 ## [23] labeling_0.4.2 GenomeInfoDbData_1.2.4 ## [25] polyclip_1.10-0 bit64_4.0.5 ## [27] farver_2.1.0 parallelly_1.28.1 ## [29] vctrs_0.3.8 generics_0.1.0 ## [31] xfun_0.25 BiocFileCache_1.14.0 ## [33] R6_2.5.1 ggbeeswarm_0.6.0 ## [35] rsvd_1.0.5 bitops_1.0-7 ## [37] spatstat.utils_2.2-0 cachem_1.0.6 ## [39] DelayedArray_0.16.3 assertthat_0.2.1 ## [41] promises_1.2.0.1 scales_1.1.1 ## [43] beeswarm_0.4.0 gtable_0.3.0 ## [45] beachmat_2.6.4 globals_0.14.0 ## [47] goftest_1.2-2 rlang_0.4.11 ## [49] splines_4.0.3 lazyeval_0.2.2 ## [51] spatstat.geom_2.2-2 yaml_2.2.1 ## [53] reshape2_1.4.4 abind_1.4-5 ## [55] httpuv_1.6.2 tools_4.0.3 ## [57] ellipsis_0.3.2 spatstat.core_2.3-0 ## [59] jquerylib_0.1.4 RColorBrewer_1.1-2 ## [61] ggridges_0.5.3 Rcpp_1.0.7 ## [63] plyr_1.8.6 sparseMatrixStats_1.2.1 ## [65] zlibbioc_1.36.0 purrr_0.3.4 ## [67] RCurl_1.98-1.4 rpart_4.1-15 ## [69] deldir_0.2-10 viridis_0.6.1 ## [71] pbapply_1.4-3 cowplot_1.1.1 ## [73] zoo_1.8-9 ggrepel_0.9.1 ## [75] cluster_2.1.0 magrittr_2.0.1 ## [77] data.table_1.14.0 RSpectra_0.16-0 ## [79] scattermore_0.7 lmtest_0.9-38 ## [81] RANN_2.6.1 fitdistrplus_1.1-5 ## [83] patchwork_1.1.1 mime_0.11 ## [85] evaluate_0.14 xtable_1.8-4 ## [87] gridExtra_2.3 compiler_4.0.3 ## [89] tibble_3.1.4 KernSmooth_2.23-17 ## [91] crayon_1.4.1 htmltools_0.5.2 ## [93] mgcv_1.8-33 later_1.3.0 ## [95] tidyr_1.1.3 DBI_1.1.1 ## [97] ExperimentHub_1.16.1 dbplyr_2.1.1 ## [99] MASS_7.3-53 rappdirs_0.3.3 ## [101] Matrix_1.3-4 cli_3.0.1 ## [103] igraph_1.2.6 pkgconfig_2.0.3 ## [105] scuttle_1.0.4 plotly_4.9.4.1 ## [107] spatstat.sparse_2.0-0 vipor_0.4.5 ## [109] bslib_0.3.0 XVector_0.30.0 ## [111] stringr_1.4.0 digest_0.6.27 ## [113] sctransform_0.3.2 RcppAnnoy_0.0.19 ## [115] spatstat.data_2.1-0 rmarkdown_2.10 ## [117] leiden_0.3.9 uwot_0.1.10 ## [119] DelayedMatrixStats_1.12.3 curl_4.3.2 ## [121] shiny_1.6.0 lifecycle_1.0.0 ## [123] nlme_3.1-149 jsonlite_1.7.2 ## [125] BiocNeighbors_1.8.2 viridisLite_0.4.0 ## [127] limma_3.46.0 fansi_0.5.0 ## [129] pillar_1.6.2 lattice_0.20-41 ## [131] fastmap_1.1.0 httr_1.4.2 ## [133] survival_3.2-7 interactiveDisplayBase_1.28.0 ## [135] glue_1.4.2 png_0.1-7 ## [137] BiocVersion_3.12.0 bit_4.0.4 ## [139] stringi_1.7.4 sass_0.4.0 ## [141] blob_1.2.2 BiocSingular_1.6.0 ## [143] AnnotationHub_2.22.1 memoise_2.0.0 ## [145] irlba_2.3.3 future.apply_1.8.1
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 79
- 80
- 81
- 82
- 83
- 84
- 85
- 86
- 87
- 88
- 89
- 90
- 91
- 92
- 93
- 94
- 95
- 96
- 97
- 98
- 99
- 100
- 101
- 102
- 103
- 104
- 105
- 106
- 107
可以看出，各种注释方法都是比较科学的，但最大的挑战来源于参考数据集与待注释数据集间的适用度。
OK 吃席！
本系列其他课程
手把手教你做单细胞测序数据分析（六）——组间差异分析及可视化
欢迎关注同名公众号~
声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/Monodyee/article/detail/345210
推荐阅读
相关标签
