单细胞分析之细胞注释-2：Garnett
摘要
Garnett是一个单细胞自动注释软件包，输入数据包括一个单细胞数据集和细胞类型定义文件。Garnett使用弹性网回归模型的机器学习算法训练出一个基于回归的分类器。随后训练好的分类器就可以用于更多数据集的细胞类型定义。
官网： https://cole-trapnell-lab.github.io/garnett/
Garnett的工作流程包括两部分：
使用
1. 安装
Garnett的运行依赖于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. 训练分类器
2.1 载入单细胞数据，处理为CDS格式
Garnett是基于Monocle3，所以它输入的数据格式是CellDataSet (CDS)。
这一部分的操作可以参考Monocle3
##加载R包
library(Seurat)
library(monocle3)
library(garnett)
library(SingleR)
library(org.Hs.eg.db)
library(tidyverse)
library(patchwork)
#载入演示数据集，依然是熟悉的pbmc3k
pbmc <- readRDS("pbmc.rds")
#创建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函数相当于seurat中NormalizeData+ScaleData+RunPCA
cds <- preprocess_cds(cds, num_dim = 50)
2.2 marker文件准备
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. 当pre-trained分类器中的genelist不能满足需求时，也可以自己准备marker基因文件。
格式：需要包含至少如下两行
1. “>”开头的细胞类型行；
2. “expressed:” 开头行，后面跟定义细胞类型的marker基因。基因之间使用,分隔。可选以下几行的附加信息：
1. “not expressed:” 开头行，后面跟定义细胞类型的负选marker基因。例如CD4+T细胞不能表达CD8，就可以写在这一行。
2. 可以使用“expressed above:”、“expressed below:”、“expressed between:”定义marker基因的表达值范围。适用于一些marker基因是根据high/low来区分细胞群的情况。比如注释Ly6ChighCCR2highCX3CR1low的单核细胞。
3. “subtype of:” 开头的字符定义细胞类型的父类，即此类细胞属于哪种细胞的亚型。
4. “references:” 开头的字符定义marker基因的选择依据。
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)
评估结果会以红色字体提示哪些marker基因在数据库中找不到对应的Ensembl名称，以及哪些基因的特异性不高（标注“High overlap with XX cells”）。我们可以根据评估结果优化marker基因，或者添加其他信息来辅助区分细胞类型。
2.4 训练分类器
# 使用marker file和cds对象训练分类器
# 这一步比较慢❗️
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")
# 查看分类器最后选择的根节点基因，注意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
# 查看分类器中分支节点的分类基因
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. 使用训练好的分类器预测自己的数据
pData(cds)$garnett_cluster <- pData(cds)$seurat_clusters
# 使用前面训练的pbmc_classifier来对自己的数据进行细胞分型
cds <- classify_cells(cds, pbmc_classifier, db = org.Hs.eg.db, cluster_extend = TRUE, cds_gene_id_type = "SYMBOL")
## 将结果返回给seurat对象# 提取分类结果
cds.meta <- subset(pData(cds), select = c("cell_type", "cluster_ext_type")) %>% as.data.frame()
pbmc <- AddMetaData(pbmc, metadata = cds.meta)
#查看结果
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)
也可以下载已有的分类器来预测
## 下载已经训练好的分类器
hsPBMC <- readRDS("hsPBMC.rds")
## 使用下载的的hsPBMC来对自己的数据进行细胞分型
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
