單細(xì)胞分析之細(xì)胞注釋-1：Azimuth

摘要

Garnett是一個(gè)單細(xì)胞自動注釋軟件包浙于，輸入數(shù)據(jù)包括一個(gè)單細(xì)胞數(shù)據(jù)集和細(xì)胞類型定義文件棍好。Garnett使用彈性網(wǎng)回歸模型的機(jī)器學(xué)習(xí)算法訓(xùn)練出一個(gè)基于回歸的分類器贬养。隨后訓(xùn)練好的分類器就可以用于更多數(shù)據(jù)集的細(xì)胞類型定義芭届。
官網(wǎng)： https://cole-trapnell-lab.github.io/garnett/

Garnett的工作流程包括兩部分：

使用

1. 安裝

Garnett的運(yùn)行依賴于monole3，因此在安裝garnett前需要先安裝monole3和其他依賴包。

# First install Bioconductor and Monocle
if (!requireNamespace("BiocManager"))
    install.packages("BiocManager")

BiocManager::install()
BiocManager::install(c("monocle"))

# Next install a few more dependencies
BiocManager::install(c('DelayedArray', 'DelayedMatrixStats', 'org.Hs.eg.db', 'org.Mm.eg.db'))

install.packages("devtools")
devtools::install_github("cole-trapnell-lab/garnett")

2. 訓(xùn)練分類器

2.1 載入單細(xì)胞數(shù)據(jù)窜醉，處理為CDS格式

Garnett是基于Monocle3荤胁，所以它輸入的數(shù)據(jù)格式是CellDataSet (CDS)铜邮。
這一部分的操作可以參考Monocle3

##加載R包
library(Seurat)
library(monocle3)
library(garnett)
library(SingleR)
library(org.Hs.eg.db)
library(tidyverse)
library(patchwork)
#載入演示數(shù)據(jù)集，依然是熟悉的pbmc3k
pbmc <- readRDS("pbmc.rds")

#創(chuàng)建CDS對象
data <- GetAssayData(pbmc, assay = 'RNA', slot = 'counts')
cell_metadata <- pbmc@meta.data
gene_annotation <- data.frame(gene_short_name = rownames(data))
rownames(gene_annotation) <- rownames(data)
cds <- new_cell_data_set(data,
                         cell_metadata = cell_metadata,
                         gene_metadata = gene_annotation)
#preprocess_cds函數(shù)相當(dāng)于seurat中NormalizeData+ScaleData+RunPCA
cds <- preprocess_cds(cds, num_dim = 50)

2.2 marker文件準(zhǔn)備

a. 直接下載
pre-trained分類器鏈接：https://cole-trapnell-lab.github.io/garnett/classifiers/寨蹋，可以直接下載這些marker基因列表（第二列）使用松蒜。

#如下載上面hsPBMC_markers.txt文件
download.file(url="https://cole-trapnell-lab.github.io/garnett/marker_files/hsPBMC_markers.txt",
              destfile = "hsPBMC_markers.txt")

b. 當(dāng)pre-trained分類器中的genelist不能滿足需求時(shí)，也可以自己準(zhǔn)備marker基因文件已旧。

格式：需要包含至少如下兩行
1. “>”開頭的細(xì)胞類型行秸苗；
2. “expressed:” 開頭行，后面跟定義細(xì)胞類型的marker基因运褪【ィ基因之間使用,分隔玖瘸。

可選以下幾行的附加信息：
1. “not expressed:” 開頭行，后面跟定義細(xì)胞類型的負(fù)選marker基因檀咙。例如CD4⁺T細(xì)胞不能表達(dá)CD8雅倒，就可以寫在這一行。
2. 可以使用“expressed above:”弧可、“expressed below:”蔑匣、“expressed between:”定義marker基因的表達(dá)值范圍。適用于一些marker基因是根據(jù)high/low來區(qū)分細(xì)胞群的情況棕诵。比如注釋Ly6C^highCCR2^highCX3CR1^low的單核細(xì)胞裁良。
3. “subtype of:” 開頭的字符定義細(xì)胞類型的父類，即此類細(xì)胞屬于哪種細(xì)胞的亞型校套。
4. “references:” 開頭的字符定義marker基因的選擇依據(jù)价脾。
5. #后可添加注釋

2.3 marker基因評估

# 對marker file中marker基因評分
marker_check <- check_markers(cds, "hsPBMC_markers.txt",
                              db=org.Hs.eg.db,
                              cds_gene_id_type = "SYMBOL",
                              marker_file_gene_id_type = "SYMBOL")
plot_markers(marker_check)

評估結(jié)果會以紅色字體提示哪些marker基因在數(shù)據(jù)庫中找不到對應(yīng)的Ensembl名稱，以及哪些基因的特異性不高（標(biāo)注“High overlap with XX cells”）笛匙。我們可以根據(jù)評估結(jié)果優(yōu)化marker基因侨把，或者添加其他信息來輔助區(qū)分細(xì)胞類型。

2.4 訓(xùn)練分類器

# 使用marker file和cds對象訓(xùn)練分類器
# 這一步比較慢??
pbmc_classifier <- train_cell_classifier(cds = cds,
                                         marker_file = "hsPBMC_markers.txt",
                                         db=org.Hs.eg.db,
                                         cds_gene_id_type = "SYMBOL",
                                         num_unknown = 50,
                                         marker_file_gene_id_type = "SYMBOL")
saveRDS(pbmc_classifier, "pbmc_classifier.rds")

# 查看分類器最后選擇的根節(jié)點(diǎn)基因妹孙，注意markerfile的基因都會在其中
feature_genes_root <- get_feature_genes(pbmc_classifier, node = "root",  
                                        db = org.Hs.eg.db, convert_ids = TRUE)
head(feature_genes_root)
#                 B cells Dendritic cells   Monocytes     NK cells    T cells     Unknown
# (Intercept) -0.58587469       2.5432315  0.39392524 -1.273097026  1.3280862 -2.40627129
# MS4A1        4.22744788      -1.7427981 -2.92309279 -0.481741931 -2.6971375  3.61732250
# CD3E         0.06193712      -0.4754740  0.32149610 -0.401413523  0.8284622 -0.33500786
# CD3D        -0.08966444      -0.9919422  0.08437099 -0.002444114  0.7886462  0.21103359
# CD3G         0.31612388      -0.6283401  0.23584751 -0.749523853  0.8628022 -0.03690959
# CD19         4.96368329      -0.5244165 -5.40305396 -0.120560423 -4.2104908  5.29483834
# 查看分類器中分支節(jié)點(diǎn)的分類基因
feature_genes_branch <- get_feature_genes(pbmc_classifier, node = "T cells",  
                                        db = org.Hs.eg.db, convert_ids = TRUE)
head(feature_genes_branch)  
#              CD4 T cells   CD8 T cells      Unknown
# (Intercept) -2.288375871 -3.213676e-01  2.609743456
# IL7R         3.927857793 -4.670505e+00  0.742647421
# PPP6C       -0.000225455  1.163007e-06  0.000224292
# IL2RA        4.370269532 -2.698651e+00 -1.671618461
# CD4          3.844007505 -3.589102e+00 -0.254905668
# CD8A        -4.915576863  5.761798e+00 -0.846220932

3. 使用訓(xùn)練好的分類器預(yù)測自己的數(shù)據(jù)

pData(cds)$garnett_cluster <- pData(cds)$seurat_clusters
# 使用前面訓(xùn)練的pbmc_classifier來對自己的數(shù)據(jù)進(jìn)行細(xì)胞分型
cds <- classify_cells(cds, pbmc_classifier, db = org.Hs.eg.db, cluster_extend = TRUE, cds_gene_id_type = "SYMBOL")
## 將結(jié)果返回給seurat對象# 提取分類結(jié)果
cds.meta <- subset(pData(cds), select = c("cell_type", "cluster_ext_type")) %>% as.data.frame()
pbmc <- AddMetaData(pbmc, metadata = cds.meta)
#查看結(jié)果
p <- DimPlot(pbmc, group.by = "cluster_ext_type", label = T, 
              label.size = 3,repel = T) +  ggtitle("Classified by Garnett")
ggsave("Garnett.pdf", p, width = 8, height = 6)

也可以下載已有的分類器來預(yù)測

## 下載已經(jīng)訓(xùn)練好的分類器
hsPBMC <- readRDS("hsPBMC.rds")
## 使用下載的的hsPBMC來對自己的數(shù)據(jù)進(jìn)行細(xì)胞分型
pData(cds)$garnett_cluster <- pData(cds)$seurat_clusters
cds <- classify_cells(cds, hsPBMC, db = org.Hs.eg.db, cluster_extend = TRUE, cds_gene_id_type = "SYMBOL")
cds.meta <- subset(pData(cds), select = c("cell_type", "cluster_ext_type")) %>% as.data.frame()
pbmc <- AddMetaData(pbmc, metadata = cds.meta)

參考：
https://cole-trapnell-lab.github.io/garnett/docs_m3/
https://mp.weixin.qq.com/s/tYNW86UsjM9quytLTxbmMA