全网最详细单细胞保姆级分析教程(三) ---细胞注释

莘薪

于 2024-07-14 09:00:00 发布

阅读量483

点赞数 1

分类专栏： RNA-seq分析文章标签： Rstudio 单细胞 RNA-seq 细胞注释 SingerR

本文链接：https://blog.csdn.net/ZxVSaccount/article/details/140408587

版权

RNA-seq分析专栏收录该内容

8 篇文章 6 订阅

订阅专栏

前面我们学习了对单细胞流程化的分析过程,也了解了seurat对象的基本结构,从之前的学习中,我们知道细胞注释是及其重要的一部分,因此接下来让我们详细探讨一下细胞注释有哪些方法吧!!

1. 数据准备(仿照单细胞流程)

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) 
}

数据结构示意图

请添加图片描述

2. 细胞注释

2.1 查数据库

# 查看数据
library(dplyr)
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n=10,wt=avg_log2FC)
VlnPlot(pbmc,features = top10$gene[1:20],pt.size = 0)

请添加图片描述

# 通过标记基因及文献,可以人工确定各分类群的细胞类型,则可以如手动添加细胞群名称
new.cluster.ids <- c("Naive CD4 T", "Memory CD4 T", "CD14+ Mono", "B", "CD8 T", "FCGR3A+ Mono", 
                     "NK", "DC", "Platelet")
# 帮助单细胞测序进行注释的数据库：
#1.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
#2.cellmarker:http://biocc.hrbmu.edu.cn/CellMarker/index.jsp
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc,new.cluster.ids)
DimPlot(pbmc,label = T)

请添加图片描述

2.2 利用SingerR注释

这一种方法很鸡肋,主要是由于很难找到合适的数据集,但是对于免疫细胞的注释还是有可观的效果的

if(!require(SingleR))BiocManager::install(SingleR)
if(!require(matrixStats))BiocManager::install('matrixStats')
if(!require(celldex))BiocManager::install('celldex') # 有一些版本冲突，需要重新安装一些包。
library(SingleR)
library(matrixStats)
library(celldex)
cg <- ImmGenData()#选取我们要使用的免疫细胞参考数据集

# 取出样本中的表达序列
bfreaname.pbmc <- pbmc
# GetAssayData -- 去除样本中的assay数据进行分析
assay_for_SingleR <- GetAssayData(bfreaname.pbmc, slot="data")
# labels是参考数据集中细胞类型的标签
predictions <- SingleR(test=assay_for_SingleR,ref=cg, labels=cg$label.main)
table(predictions$labels) # 看看都注释到了哪些细胞
          B cells Endothelial cells       Eosinophils        Mast cells          NK cells 
             2555                19                18                25                11 
    Stromal cells 
               10 
# 可以明显看到结果并不是很理想,原因是没有合适的数据集

#得到seurat中编号与预测标签之间的关系
cellType=data.frame(seurat=bfreaname.pbmc@meta.data$seurat_clusters,
                    predict=predictions$labels)
sort(table(cellType[,1]))
table(cellType[,1:2])#访问celltyple的2~3列
seurat B cells Endothelial cells Eosinophils Mast cells NK cells Stromal cells
     0     679                 5           4          0        3             1
     1     461                 2           6         12        2             0
     2     435                 4           6          0        2             2
     3     342                 0           0          0        0             0
     4     308                 4           0          0        1             0
     5     134                 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

2.3 利用自制的SingerR数据集进行注释

library(SingleR)
library(ggplot2)
library(Seurat)
library(textshape)
library(devtools)
library(scater)
library(dplyr)
# 这里为了检测，我们将参考数据集与目标数据集用同一个数据进行测试
myref <- pbmc
myref$celltype <- Idents(myref)
> table(Idents(myref))
 Naive CD4 T Memory CD4 T   CD14+ Mono            B        CD8 T FCGR3A+ Mono           NK 
         692          483          449          342          313          159          154 
          DC     Platelet 
          32           14 

# 对 myref 对象中的 RNA 表达数据进行对数转换，并将结果存储在 Refassay 变量中
Refassay <- log1p(AverageExpression(myref, verbose = FALSE)$RNA)

ref_sce <- SingleCellExperiment::SingleCellExperiment(assays=list(counts=Refassay))
#参考数据集需要构建成一个SingleCellExperiment对象
ref_sce=scater::logNormCounts(ref_sce)

# 提取数据归一化后的counts的前四行和前四列,并将其显示出来
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.00948295   0.06667557 0.009128802 .        
# AP006222.2     .            0.01549977 0.012607874 .        
# RP11-206L10.2  0.01151648   .          .           0.0296705
# RP11-206L10.9  .            0.02064656 .           .        

colData(ref_sce)$Type=colnames(Refassay)
ref_sce#构建完成

# 提取自己的单细胞矩阵
testdata <- GetAssayData(myref, slot="data")
testdata <- GetAssayData(myref, slot="data")
pred <- SingleR(test=testdata, ref=ref_sce, 
                labels=ref_sce$Type,
                #clusters = scRNA@active.ident
)

# 整理数据
table(pred$labels)
head(pred) 
as.data.frame(table(pred$labels))
#           Var1 Freq
# 1            B  343
# 2   CD14+ Mono  556
# 3        CD8 T  306
# 4           DC   30
# 5 FCGR3A+ Mono  178
# 6 Memory CD4 T  466
# 7  Naive CD4 T  586
# 8           NK  159
# 9     Platelet   14

cellType=data.frame(seurat=pbmc@meta.data$seurat_clusters,
                    predict=pred$labels)#得到seurat中编号与预测标签之间的关系
sort(table(cellType[,1]))
table(cellType[,1:2])#访问celltyple的2~3列
     predict
seurat   B CD14+ Mono CD8 T  DC FCGR3A+ Mono Memory CD4 T Naive CD4 T  NK Platelet
     0   2        155     8   0            0            0         527   0        0
     1   0          0     0   1           20          462           0   0        0
     2   0        381    11   0            0            0          57   0        0
     3 341          1     0   0            0            0           0   0        0
     4   0         19   284   0            0            0           2   8        0
     5   0          0     0   0          158            1           0   0        0
     6   0          0     3   0            0            0           0 151        0
     7   0          0     0  29            0            3           0   0        0
     8   0          0     0   0            0            0           0   0       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的对应预测细胞类型

testname <- 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.4 映射处理

# 利用seurat内置的原先用于细胞整合的功能，将参考数据与待注释数据进行映射处理
library(Seurat)
pancreas.query <- pbmc # 待注释数据

# 这行代码使用了 Seurat 包中的 FindTransferAnchors 函数来寻找在参考数据集 pbmc 和查询数据集 pancreas.query 之间的转移锚点（transfer anchors）
# 转移锚点是在不同数据集之间建立对应关系的关键元素。通过找到这些锚点，可以将参考数据集的信息传递到查询数据集上，进行细胞类型注释、数据整合等操作。执行这行代码后，pancreas.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 对象添加元数据
pancreas.query <- AddMetaData(pancreas.query, metadata = predictions)

#把注释加回原来的数据集
pancreas.query$prediction.match <- pancreas.query$predicted.id
table(pancreas.query$prediction.match)

Idents(pancreas.query)<- 'prediction.match'

p1 <- DimPlot(pbmc,label = T)
p2 <- DimPlot(pancreas.query,label = T)#比较一下测试数据与参考数据集之间有没有偏差
p1|p2 # 完美，无差别注释，当然了，我们这个参考数据用的比较极端