本筆記來源于B站@生信技能樹-jimmy；學(xué)習(xí)視頻鏈接: 「生信技能樹」單細(xì)胞數(shù)據(jù)挖掘

3. cell annotation

# 對(duì)腫瘤細(xì)胞來說，分群后的細(xì)胞亞群注釋是不可行的
# 這里僅僅是演示 SingleR 做 cell annotation的流程
library(SingleR)
??singleR
library(Seurat)

load("../../scRNA.Rdata")

refdata <- get(load("../../rawdata/HumanPrimaryCellAtlasData.Rdata"))
assay(refdata)[1:4,1:4]
head(refdata@colData)
head(refdata)

ref <- HumanPrimaryCellAtlasData() #下載參考數(shù)據(jù)庫，等待時(shí)間較長。建議下載后儲(chǔ)存為Rdata，方便以后使用。

testdata <- GetAssayData(scRNA, slot="data")
clusters <- scRNA@meta.data$seurat_clusters
cellpred <- SingleR(test = testdata, ref = refdata, labels = refdata$label.fine, 
                    # label.fine耗時(shí)比較長一點(diǎn)
                    method = "cluster", clusters = clusters, 
                    assay.type.test = "logcounts", assay.type.ref = "logcounts")
save(cellpred,file = "../../cellpred.Rdata")
load("../../tmp/cellpred.Rdata")
rm(refdata, HumanPrimaryCellAtlasData, testdata) #珍惜內(nèi)存

table(cellpred$labels)
celltype = data.frame(ClusterID=rownames(cellpred), celltype=cellpred$labels, stringsAsFactors = F)
table(celltype$ClusterID,celltype$celltype)  

SingleCell8_1.jpg

# 結(jié)合上述結(jié)果，給scRNA增添celltype注釋信息
# 先新增列celltype，值均為NA，然后利用下一行代碼循環(huán)填充
scRNA@meta.data$celltype = "NA"

for(i in 1:nrow(celltype)){
  scRNA@meta.data[which(scRNA@meta.data$seurat_clusters == celltype$ClusterID[i]),'celltype'] <- celltype$celltype[i]}

# 繪圖
p1 <- DimPlot(scRNA, group.by="celltype", label=F , reduction='tsne')
p1
ggsave("../../out/3.3celltype_anno.pdf", plot = p1, width = 18, height = 12) 

Fig2_1.png