單細胞分析實錄(14): 細胞類型注釋的另一種思路 — CellID

之前寫過一篇細胞類型注釋的推文：單細胞分析實錄(7): 差異表達分析/細胞類型注釋，里面介紹了細胞類型注釋的兩種思路。在談到自動注釋時，我們基本都會用SingleR這個軟件，類似的軟件還有很多，原理上可能五花八門，但表現不一定有SingleR好。很多軟件在做注釋時會基於某個細胞和參考細胞表達譜的相似性，今天看到的這個軟件使用的是另一種思路：將某個細胞的特征基因集與表征細胞類型的參考基因集做富集分析，當在某個細胞類型的基因集上顯著富集時，就定義為這個細胞類型（和GO富集很像）。所以這種方法的關鍵點是：找單個細胞的特征基因集，以及找細胞類型的參考基因集。

CellID發表在NBT上面：Gene signature extraction and cell identity recognition at the single-cell level with Cell-ID。解決的主要是找單個細胞的特征基因集這個問題。

軟件的官方文檔：https://bioconductor.org/packages/devel/bioc/vignettes/CelliD/inst/doc/BioconductorVignette.html
有詳細的演示，並且在選取細胞類型參考基因集時，使用的是網站PanglaoDB：https://panglaodb.se/index.html
這個網站有人和小鼠各種細胞類型的marker，這是下載鏈接：https://panglaodb.se/markers/PanglaoDB_markers_27_Mar_2020.tsv.gz
比我之前整理的詳細得多。

題外話：在單細胞分析實錄(7)發布之后，公眾號后台斷斷續續收到幾十位讀者詢問marker清單的信息，希望這個更全面的清單能幫助到你。

下面我簡單演示一下CellID的用法：

library(CellID)
library(Seurat)
library(tidyverse)

#獲取原始表達矩陣
BaronMatrix   <- readRDS(url("https://storage.googleapis.com/cellid-cbl/BaronMatrix.rds"))
#僅考慮編碼蛋白質的基因
data("HgProteinCodingGenes")
BaronMatrixProt <- BaronMatrix[rownames(BaronMatrix) %in% HgProteinCodingGenes,]

#這幾步類似Seurat的標准流程
Baron <- CreateSeuratObject(counts = BaronMatrixProt, project = "Baron", min.cells = 5)
Baron <- NormalizeData(Baron)
Baron <- ScaleData(Baron, features = rownames(Baron)) 
Baron <- RunMCA(Baron) #該軟件的主要分析函數，將細胞和基因同時降維到一個空間，離細胞近的基因被定義為細胞的特征基因——令人窒息的操作

#下載參考基因集
panglao <- read_tsv("https://panglaodb.se/markers/PanglaoDB_markers_27_Mar_2020.tsv.gz")
#根據器官過濾，示例數據就是胰腺的
panglao_pancreas <- panglao %>% filter(organ == "Pancreas")
#物種包含人
panglao_pancreas <- panglao_pancreas %>%  filter(str_detect(species,"Hs"))
#下面兩步會得到一個列表，列表的每一個元素是由基因集構成的字符向量，且每個元素被命名為對應的細胞類型的名字
panglao_pancreas <- panglao_pancreas %>%  
  group_by(`cell type`) %>%  
  summarise(geneset = list(`official gene symbol`))
pancreas_gs <- setNames(panglao_pancreas$geneset, panglao_pancreas$`cell type`)

#富集分析，用的超幾何檢驗
HGT_pancreas_gs <- RunCellHGT(Baron, pathways = pancreas_gs, dims = 1:50, n.features = 200) #每個細胞選200個特征基因

富集分析之后會返回一個矩陣，矩陣中的每一個值表示p值在校正之后，求-log10，即-log10 corrected p-value，當然這個值越大越好。

> HGT_pancreas_gs[1:4,1:4]
4 x 4 sparse Matrix of class "dgCMatrix"
             human1_lib1.final_cell_0001 human1_lib1.final_cell_0002
Acinar cells                  60.6897442                  63.2674918
Alpha cells                    0.1241303                   0.1241303
Beta cells                     .                           .        
Delta cells                    .                           .        
             human1_lib1.final_cell_0003 human1_lib1.final_cell_0004
Acinar cells                  58.1720116                  63.2674918
Alpha cells                    0.1910771                   0.1770659
Beta cells                     0.4326425                   0.1770659
Delta cells                    .                           .

pancreas_gs_prediction <- rownames(HGT_pancreas_gs)[apply(HGT_pancreas_gs, 2, which.max)]
#矩陣的每一列依次：判斷是否是最大值，得到一列布爾值，結合矩陣的行名會返回行名中的一個元素（也就是最有可能的細胞類型）。
#所有列運行完了之后會得到所有細胞最可能的注釋
pancreas_gs_prediction_signif <- ifelse(apply(HGT_pancreas_gs, 2, max)>2, yes = pancreas_gs_prediction, "unassigned")
#如果`-log10 corrected p-value`的值小於等於2，則認為不顯著，注釋更正為"unassigned"
Baron$pancreas_gs_prediction <- pancreas_gs_prediction_signif

head([e-mail],2)
                            orig.ident nCount_RNA nFeature_RNA
human1_lib1.final_cell_0001     human1      21956         3362
human1_lib1.final_cell_0002     human1      27309         4009
                            pancreas_gs_prediction
human1_lib1.final_cell_0001           Acinar cells
human1_lib1.final_cell_0002           Acinar cells

到這兒，細胞類型預測算結束了。再來用真正的細胞類型（也可以是手動注釋的結果），評估一下軟件預測出來的准不准。走Seurat標准流程：

#加載預先知道的注釋信息
BaronMetaData <- readRDS(url("https://storage.googleapis.com/cellid-cbl/BaronMetaData.rds"))
[e-mail]=merge([e-mail],BaronMetaData,by="row.names")
rownames([e-mail])=[e-mail]$Row.names
[e-mail]$Row.names=NULL
head([e-mail],2)
                            orig.ident nCount_RNA nFeature_RNA
human1_lib1.final_cell_0001     human1      21956         3362
human1_lib1.final_cell_0002     human1      27309         4009
                            pancreas_gs_prediction      donor    cell.type
human1_lib1.final_cell_0001           Acinar cells GSM2230757 Acinar cells
human1_lib1.final_cell_0002           Acinar cells GSM2230757 Acinar cells

Baron <- FindVariableFeatures(Baron, selection.method = "vst", nfeatures = 2000)
Baron <- ScaleData(Baron)
Baron <- RunPCA(Baron, npcs = 50, verbose = FALSE)
Baron <- FindNeighbors(Baron, dims = 1:30)
Baron <- FindClusters(Baron, resolution = 0.5)
Baron <- RunUMAP(Baron, dims = 1:30)
Baron <- RunTSNE(Baron, dims = 1:30)

library(cowplot)
p1 <- DimPlot(Baron, reduction = "tsne", group.by = "cell.type", pt.size=0.5, label = TRUE,repel = TRUE)
p2 <- DimPlot(Baron, reduction = "tsne", group.by = "pancreas_gs_prediction", pt.size=0.5, label = TRUE,repel = TRUE)
fig_tsne <- plot_grid(p1, p2, labels = c('cell.type','pancreas_gs_prediction'),align = "v",ncol = 2)
ggsave(filename = "tsne.pdf", plot = fig_tsne, device = 'pdf', width = 40, height = 15, units = 'cm')

---

這種細胞注釋方法非常依賴於細胞類型的參考基因集的質量，我拿自己的數據集跑這個軟件，並沒有該軟件的示例數據效果好，原因就是前面提到的PanglaoDB網站中我這個方向的參考數據質量不高。大家還知道類似的提供了細胞類型及對應的基因集網站嗎？歡迎討論，集思廣益，我會在之后的推文中列出大家推薦的網站。

另外，這種用富集來注釋的思路可以延伸一下，比如對於某個cluster不太清楚是什么類型，可以找完差異基因后，用細胞類型參考基因集做富集分析，來定義cluster類型。

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