使用Garnett包进行单细胞类型分类注释分析（一）：Train
image
Garnett包简介
Garnett是一个根据单细胞表达数据自动进行细胞类型分类注释的R包。Garnett使用单细胞表达数据和细胞类型定义（marker）文件，训练一个基于回归的分类器。一旦针对某一组织/样本类型训练成了一个分类器，它就可以将其应用于对其他相似组织的数据集进行分类注释。除了描述训练和分类的功能外，该网站还旨在成为一个存储以前训练出来的分类器仓库。
Garnett及其依赖包的安装
Garnett runs in the R statistical computing environment. You will need R version 3.5 or higher to install Garnett.
# Install from Github
# Garnett builds upon a package called Monocle. Before installing Garnett, first install Monocle using Bioconductor:
# First install Bioconductor and Monocle
if (!requireNamespace("BiocManager"))
install.packages("BiocManager")
# 安装monocle包
BiocManager::install(c("monocle"))
# Next install a few more dependencies
# 安装依赖包
BiocManager::install(c('DelayedArray', 'DelayedMatrixStats', 'org.Hs.eg.db', 'org.Mm.eg.db'))
# 使用devtools安装garnett包
install.packages("devtools")
devtools::install_github("cole-trapnell-lab/garnett")
# Load Garnett
library(garnett)
Garnett包的工作流程
image
Garnett包的工作流程主要包括以下两个部分：
-
Train/obtain the classifier: 我们可以下载现有的细胞分类器，或者训练自己的分类器。训练分类器时，Garnett首先解析一个标记文件（marker file），并选择一组训练细胞，然后训练一个多项分类器来区分不同的细胞类型。
-
Classify cells: 接下来，Garnett将训练好的分类器应用于一组细胞，以生成细胞类型分配。Garnett还可以选择将分类扩展到类似的细胞，以生成一组独立的分群扩展类型分配。
1a. Using a pre-trained classifier 使用一个预先训练好的分类器
Garnett包工作的第一步是构建一个细胞类型分类器，我们既可以使用预先训练好的分类器，也可以单独训练一个自己的分类器。
image
我们已经为各种生物和组织生成了一系列预训练好的分类器。如果您的数据类型存在预先训练好的分类器，我们建议您尝试一下，可以在此处找到可用分类器的列表（https://cole-trapnell-lab.github.io/garnett/classifiers/）。
image
要使用预先训练好的分类器，请首先下载分类器，然后使用以下命令将其加载到您的R会话中：
classifier <- readRDS("path/to/classifier.RDS")
1b. Train your own classifier 训练一个自己的分类器
如果我们使用的组织类型不存在分类器，或者我们的数据中不包含预期的细胞类型，则需要训练生成一个自己的分类器。
image
Loading your data
首先，我们加载单细胞的表达数据，并使用它训练一个细胞分类器。这里，我们使用10x平台测序的PBMC数据进行演示。
# load in the data
# NOTE: the 'system.file' file name is only necessary to read in included package data
library(monocle)
# 读取表达矩阵
mat <- Matrix::readMM(system.file("extdata", "exprs_sparse.mtx", package = "garnett"))
head(mat)
6 x 800 sparse Matrix of class "dgTMatrix"
[1,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
[2,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
[3,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
[4,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
[5,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
[6,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
.....suppressing 768 columns in show(); maybe adjust 'options(max.print= *, width = *)'
..............................
# 读取基因的metadata信息
fdata <- read.table(system.file("extdata", "fdata.txt", package = "garnett"))
head(fdata)
gene_short_name num_cells_expressed
MIR1302-10 MIR1302-10 0
FAM138A FAM138A 0
OR4F5 OR4F5 0
RP11-34P13.7 RP11-34P13.7 0
RP11-34P13.8 RP11-34P13.8 0
AL627309.1 AL627309.1 0
# 读取细胞的metadata信息
pdata <- read.table(system.file("extdata", "pdata.txt", package = "garnett"), sep="\t")
head(pdata)
tsne_1 tsne_2 Size_Factor FACS_type
AAGCACTGCACACA-1 3.840315 12.084191 0.5591814 B cells
GGCTCACTGGTCTA-1 9.970962 3.505393 0.5159340 B cells
AGCACTGATATCTC-1 3.459529 4.935273 0.6980284 B cells
ACACGTGATATTCC-1 1.743949 7.782671 0.8156310 B cells
ATATGCCTCTGCAA-1 5.783448 8.558898 1.1153280 B cells
TGACGAACCTATTC-1 10.792853 10.585274 0.6494699 B cells
row.names(mat) <- row.names(fdata)
colnames(mat) <- row.names(pdata)
# create a new CDS object
# 构建CDS对象
pd <- new("AnnotatedDataFrame", data = pdata)
fd <- new("AnnotatedDataFrame", data = fdata)
pbmc_cds <- newCellDataSet(as(mat, "dgCMatrix"),
phenoData = pd,
featureData = fd)
pbmc_cds
CellDataSet (storageMode: environment)
assayData: 32738 features, 800 samples
element names: exprs
protocolData: none
phenoData
sampleNames: AAGCACTGCACACA-1 GGCTCACTGGTCTA-1 ...
GTTGACGACCTTAT-1 (800 total)
varLabels: tsne_1 tsne_2 Size_Factor FACS_type
varMetadata: labelDescription
featureData
featureNames: MIR1302-10 FAM138A ... AC002321.1 (32738 total)
fvarLabels: gene_short_name num_cells_expressed
fvarMetadata: labelDescription
experimentData: use 'experimentData(object)'
Annotation:
# generate size factors for normalization later
pbmc_cds <- estimateSizeFactors(pbmc_cds)
Constructing a marker file
除了表达数据外，我们还需要另一个主要的输入文件（marker file）。标记文件中包含着细胞类型定义的列表，该文件告诉Garnett如何选择细胞来训练分类模型。每个细胞类型定义均以“>”号和细胞类型名称开头，后接一系列带有定义信息的行。定义行以关键字和“:”开头，不同条目之间用逗号分隔。
A simple valid example:
>B cells
expressed: CD19, MS4A1
>T cells
expressed: CD3D
通常，每个细胞类型定义可以包含三个主要部分，其中第一个部分是必须的。
Define the expression
细胞类型定义的第一个也是最重要的说明是其表达的marker基因。Garnett提供了几种用于指定marker基因的选项，详细信息如下。
image
以上是为细胞类型指定marker基因的默认方法。使用此规范时，Garnett会计算每个细胞的marker基因的得分，并考虑到总体的表达水平和测序深度。
image
这是指定marker基因表达的另一种方法，如果我们知道了期望基因表达占据的确切范围，则该方法很有用。但是，一般而言，我们不建议使用这些规范，因为它们不会考虑每个细胞中的测序深度和整体的表达情况。
Define the meta data
除了表达信息，我们还可以使用元数据进一步细分细胞的类型定义。我们可以在此处指定数据中期望的任何子类型。
image
-
subtype of: 允许我们指定某一细胞类型是定义文件中另一种细胞类型的子类型。 -
custom meta data: 允许我们为细胞类型定义提供任何其他的meta信息。
Provide your evidence
最后，我们强烈建议您记录如何选择的这些marker基因。为了便于跟踪，我们提供了一个附加的规范-reference：- 存储每种细胞类型的引文信息。添加一组URL或DOI，它们将包含在您的分类器中。
Add any comments
类似的，我们还提供了注释字符＃号，因此我们还可以在细胞类型标记文件中添加一些注释信息。Garnett会忽略掉带有＃号行之后的所有内容。
A more complex example:
>B cells
expressed: CD19, MS4A1
expressed above: CD79A 10
references: https://www.abcam.com/primary-antibodies/b-cells-basic-immunophenotyping, 10.3109/07420528.2013.775654
>T cells
expressed: CD3D
sample: blood # A meta data specification
>Helper T cells
expressed: CD4
subtype of: T cells
references: https://www.ncbi.nlm.nih.gov/pubmed/?term=12000723
Checking your markers
由于定义标记文件通常是该过程中最困难的部分，因此Garnett包提供了用于检查标记是否能正常工作的函数。相关的两个函数是check_markers和plot_markers。check_markers函数将生成有关标记的信息表，plot_markers函数绘制出最相关的信息。
>B cells
expressed: CD19, MS4A1, CD79A, ACTN, ACTB
references: https://www.ncbi.nlm.nih.gov/pubmed/?term=23688120,
https://www.ncbi.nlm.nih.gov/pubmed/?term=21149806
>T cells
expressed: CD3D, CD3E, CD3G, PTPRC
references: https://www.ncbi.nlm.nih.gov/pubmed/?term=1534551
>CD4 T cells
expressed: CD4, FOXP3, IL2RA, IL7R
subtype of: T cells
>CD8 T cells
expressed: CD8A, CD8B
subtype of: T cells
image
library(org.Hs.eg.db)
marker_file_path <- system.file("extdata", "pbmc_bad_markers.txt",
package = "garnett")
# 使用check_markers函数检查标记基因的信息
marker_check <- check_markers(pbmc_cds, marker_file_path,
db=org.Hs.eg.db,
cds_gene_id_type = "SYMBOL",
marker_file_gene_id_type = "SYMBOL")
head(marker_check)
marker_gene gene_id parent cell_type in_cds nominates
1 CD3D ENSG00000167286 root T cells TRUE 269
2 CD3E ENSG00000198851 root T cells TRUE 295
3 CD3G ENSG00000160654 root T cells TRUE 104
4 PTPRC ENSG00000081237 root T cells TRUE 179
5 CD4 ENSG00000010610 root T cells TRUE 25
6 FOXP3 ENSG00000049768 root T cells TRUE 8
total_nominated exclusion_dismisses inclusion_ambiguates most_overlap
1 486 14 10 B cells
2 486 24 16 B cells
3 486 1 1 B cells
4 486 97 94 B cells
5 486 0 0 <NA>
6 486 0 0 <NA>
ambiguity marker_score summary
1 0.037174721 11.7329351 Ok
2 0.054237288 9.4492763 Ok
3 0.009615385 10.9093843 Ok
4 0.525139665 0.6882554 High ambiguity?
5 0.000000000 5.1440329 Ok
6 0.000000000 1.6460905 Low nomination?
# 使用plot_markers函数对marker基因的信息进行可视化
plot_markers(marker_check)
image
这个标记图中提供了一些关于所选marker基因是否正确的关键信息。首先，红色标记的“not in db”告诉我们ACTN marker基因在org.Hs.eg.db注释中没有提供“SYMBOL”信息。接下来，x轴显示了每个marker基因的模糊度评分（ambiguity score）。在本例中，ACTB和PTPRC marker基因具有较高的模糊度得分，应将其排除。
注意：check_markers函数输出的值和plot_markers函数绘制的值是分类器选择的细胞数量的估计值。然而，它使用启发式快速找到候选细胞，并不能完全匹配标记所选择的细胞。请使用这些数字作为相对的度量，而不是训练集的绝对表示。
#pbmc_test.txt
>B cells
expressed: CD19, MS4A1, CD79A
references: https://www.ncbi.nlm.nih.gov/pubmed/?term=23688120, https://www.ncbi.nlm.nih.gov/pubmed/?term=21149806
>T cells
expressed: CD3D, CD3E, CD3G
references: https://www.ncbi.nlm.nih.gov/pubmed/?term=1534551
>CD4 T cells
expressed: CD4, FOXP3, IL2RA, IL7R
subtype of: T cells
>CD8 T cells
expressed: CD8A, CD8B
subtype of: T cells
image
Train the classifier
构建好细胞类型marker标记文件后，我们就可以使用train_cell_classifier函数训练分类器了，其参数与check_markers函数的参数类似。下面，我们更改了默认的一个num_unknown参数。 该参数告诉Garnett应该与多少个outgroup细胞进行比较，默认值为500，但是在此数据集中细胞数量较少，因此也设置的小一些。
library(org.Hs.eg.db)
set.seed(260)
marker_file_path <- system.file("extdata", "pbmc_test.txt",
package = "garnett")
# 使用train_cell_classifier函数进行训练构建分类器
pbmc_classifier <- train_cell_classifier(cds = pbmc_cds,
marker_file = marker_file_path,
db=org.Hs.eg.db,
cds_gene_id_type = "SYMBOL",
num_unknown = 50,
marker_file_gene_id_type = "SYMBOL")
There are 4 cell type definitions
training_sample
B cells T cells Unknown
200 200 50
Model training finished.
training_sample
CD4 T cells CD8 T cells Unknown
72 44 50
Model training finished.
# 查看训练的分类结果
pbmc_classifier
An object of class "garnett_classifier"
Slot "classification_tree":
IGRAPH 74f63cf DN-- 5 4 --
+ attr: name (v/c), classify_func (v/x), model (v/x)
+ edges from 74f63cf (vertex names):
[1] root ->B cells root ->T cells T cells->CD4 T cells
[4] T cells->CD8 T cells
Slot "cell_totals":
[1] 2.615075
Slot "gene_id_type":
[1] "ENSEMBL"
Slot "references":
$`B cells`
[1] "https://www.ncbi.nlm.nih.gov/pubmed/?term=23688120"
[2] "https://www.ncbi.nlm.nih.gov/pubmed/?term=21149806"
$`T cells`
[1] "https://www.ncbi.nlm.nih.gov/pubmed/?term=1534551"
head(pData(pbmc_cds))
# tsne_1 tsne_2 Size_Factor FACS_type
# AAGCACTGCACACA-1 3.840315 12.084191 0.5591814 B cells
# GGCTCACTGGTCTA-1 9.970962 3.505393 0.5159340 B cells
# AGCACTGATATCTC-1 3.459529 4.935273 0.6980284 B cells
# ACACGTGATATTCC-1 1.743949 7.782671 0.8156310 B cells
# ATATGCCTCTGCAA-1 5.783448 8.558898 1.1153280 B cells
# TGACGAACCTATTC-1 10.792853 10.585274 0.6494699 B cells
运行完train_cell_classifier函数后，输出对象中的“garnett_classifier”类型将包含对细胞进行分类所需的所有信息。
Viewing the classification genes
Garnett使用多项式弹性网络回归（multinomial elastic-net regression）训练进行细胞类型分类。这意味着选择某些基因作为区分细胞类型的相关基因。选择哪些基因可能会令人感兴趣，因此Garnett提供了访问所选基因的功能。
Garnett使用get_feature_genes函数查看所选的相关基因，参数是分类器，想查看哪个节点（node）(如果分类树是分层的)—使用“root”作为顶部节点，其他节点使用父细胞类型名称，使用db参数指定参考物种。如果设置convert_ids = TRUE，则该函数将自动将基因id转换为SYMBOL。
feature_genes <- get_feature_genes(pbmc_classifier,
node = "root",
db = org.Hs.eg.db)
head(feature_genes)
# B cells T cells Unknown
#(Intercept) -6.51829349 2.6403157858 3.877977704
#CD74 0.04660324 -0.0529185035 0.006315261
#MS4A1 2.42667982 -2.2194139450 -0.207265871
#CD19 4.33990156 -2.4210500042 -1.918851560
#CD79A 1.09667783 -0.7008461805 -0.395831654
#IGLL5 -0.00103241 -0.0001045559 0.001136966
Viewing references
上面，我们解释了如何在标记文件中包含关于如何选择标记的文档。为了获取这些信息—查看如何为已经训练好的分类器选择标记—可以使用get_classifier_references函数。除了分类器之外，还有一个额外的cell_type可选参数。如果传递细胞类型的名称，则只会打印该细胞类型的引用，否则将全部打印。
get_classifier_references(pbmc_classifier)
#$`B cells`
#[1] "https://www.ncbi.nlm.nih.gov/pubmed/?term=23688120" "https://www.ncbi.nlm.nih.gov/pubmed/?term=21149806"
#$`T cells`
#[1] "https://www.ncbi.nlm.nih.gov/pubmed/?term=1534551"
Submitting a classifier
我们鼓励您向我们提交高质量的分类器，以便我们可以将其提供给社区。为此，请打开一个特刊并在Garnett github存储库中填写表单。单击此处(https://github.com/cole-trapnell-lab/garnett/issues)，然后单击“New issue”按钮开始使用！
参考来源：https://cole-trapnell-lab.github.io/garnett/docs/#1b-train-your-own-classifier
