SingleR如何使用自定義的參考集

TOP生物資訊發表於2021-10-23

在我之前的帖子單細胞分析實錄(7): 差異表達分析/細胞型別註釋裡面，我已經介紹瞭如何使用SingleR給單細胞資料做註釋，當時只講了SingleR配套的參考集。這次就講講如何使用自己定義/找到的基因集做註釋，使用場景還是比較多的，比如想根據某篇論文裡面的註釋結果，給自己的資料做註釋。本文配套的視訊講解已上傳到B站，新手UP: TOP菌。gongzhong號後臺回覆20211023可獲取本文所用到的示例資料和程式碼。

本次演示用到的資料集來自2020年發表在Nature Genetics的一篇結直腸癌研究。文中用到了韓國患者(SMC)和比利時患者(KUL3)兩套資料集，兩套資料平行分析，相互印證。

我以裡面的髓系細胞為例，用KUL3資料集中的髓系細胞為參考，來註釋SMC裡面的髓系細胞

兩套資料集中髓系細胞部分的資料，已經跑完Seurat標準流程並整理成rds檔案：SMC_mye.rds和KUL3_mye.rds。程式碼存放在0.R中，這裡就不展示了，很簡單。程式碼中fread()用於快速讀取4G矩陣檔案是可以學習的地方。

mat=fread("GSE132465_GEO_processed_CRC_10X_raw_UMI_count_matrix.txt",select = c("Index",SMC.mye.anno$Index))

之後就是建立參考集，被儲存為SummarizedExperiment物件

library(Seurat)
library(tidyverse)
library(SummarizedExperiment)
library(scuttle)

KUL3_mye=readRDS("KUL3_mye.rds")

屬性資料框中主要用到cell index和cell annotation這兩個資訊，此外還需匯入表達矩陣

KUL3_mye_count=KUL3_mye[["RNA"]]@counts

pdata=KUL3_mye@meta.data[,c("Index","Cell_subtype")]
rownames(pdata)=pdata$Index
pdata$Index=NULL
colnames(pdata)="ref_label"

KUL3_mye_SE <- SummarizedExperiment(assays=list(counts=KUL3_mye_count),colData = pdata) 
#建立SummarizedExperiment物件


KUL3_mye_SE <- logNormCounts(KUL3_mye_SE) 
#Compute log-transformed normalized expression values，要有這一行，對於單細胞資料normalize之後會再算一個log
saveRDS(KUL3_mye_SE,"KUL3_mye_SE.ref.rds")

log normalize之後就多了一個logcounts的assay，singleR的官網示例都是基於logcounts這種資料

將這個參考物件儲存為rds檔案，方便以後多次呼叫

現在匯入我們的待註釋資料集和參考資料集

library(tidyverse)
library(Seurat)
library(SingleR)
library(scuttle)
library(SummarizedExperiment)

KUL3_mye_SE=readRDS("KUL3_mye_SE.ref.rds")
SMC_mye=readRDS("SMC_mye.rds")
SMC_mye_count=SMC_mye[["RNA"]]@counts

#需要取不同資料集的基因交集
common_gene <- intersect(rownames(SMC_mye_count), rownames(KUL3_mye_SE))
KUL3_mye_SE <- KUL3_mye_SE[common_gene,]
SMC_mye_count <- SMC_mye_count[common_gene,]

#也要建立SummarizedExperiment物件，以及log normalize
SMC_mye_SE <- SummarizedExperiment(assays=list(counts=SMC_mye_count))
SMC_mye_SE <- logNormCounts(SMC_mye_SE)

#註釋程式碼
singleR_res <- SingleR(test = SMC_mye_SE, ref = KUL3_mye_SE, labels = KUL3_mye_SE$ref_label)
#匯出註釋結果
anno_df <- as.data.frame(singleR_res$labels)
anno_df$Index <- rownames(singleR_res)
colnames(anno_df)[1] <- 'ref_label_from_KUL3'

#將註釋資訊新增到Seurat物件
SMC_mye@meta.data=SMC_mye@meta.data%>%inner_join(anno_df,by="Index")
rownames(SMC_mye@meta.data)=SMC_mye@meta.data$Index
DimPlot(SMC_mye, reduction = "tsne", group.by = "Cell_subtype", pt.size=2)+
DimPlot(SMC_mye, reduction = "tsne", group.by = "ref_label_from_KUL3", pt.size=2)

左圖是原文給的註釋，右圖是依據KUL3來註釋SMC資料集，可以看出還是有一些不一致的。

像這種依據某個參考集來做註釋，對參考集質量要求很高，原文KUL3只有兩千髓系細胞，考慮到10X單細胞資料的稀疏性，這個數量是不夠的。也可能是軟體的限制，在做小類註釋時，類與類之間的表達特徵其實是比較相似的，軟體不一定能精確給出合適的label，相比之下，軟體做大類註釋一般比較準。

因水平有限，有錯誤的地方，歡迎批評指正！

相關文章