如果您的組織類型不存在分類器我們的倉(cāng)庫(kù)中，或者數(shù)據(jù)中不包含您期望的細(xì)胞類型，那么您需要生成自己的分類器。
訓(xùn)練分類器的第一步是加載單細(xì)胞數(shù)據(jù)。
因?yàn)镚arnett 建立在 Monocle上，所以Garnett 的數(shù)據(jù)保存在CellDataSet (CDS)類的對(duì)象中。這個(gè)類派生自Bioconductor ExpressionSet類，它為那些分析過(guò)生物微陣列實(shí)驗(yàn)的人提供了一個(gè)常見(jiàn)的接口。Monocle提供了關(guān)于如何生成輸入cds的詳細(xì)文檔here。
例如，Garnett包含一個(gè)來(lái)自PBMC 10x V1表達(dá)式數(shù)據(jù)的小數(shù)據(jù)集.
library(monocle3)
library(garnett)
# load in the data
# NOTE: the 'system.file' file name is only necessary to read in
# included package data
#
mat <- Matrix::readMM(system.file("extdata", "exprs_sparse.mtx", package = "garnett"))
fdata <- read.table(system.file("extdata", "fdata.txt", package = "garnett"))
pdata <- read.table(system.file("extdata", "pdata.txt", package = "garnett"),
sep="\t")
row.names(mat) <- row.names(fdata)
colnames(mat) <- row.names(pdata)
# create a new CDS object
#pd <- new("AnnotatedDataFrame", data = pdata)
#fd <- new("AnnotatedDataFrame", data = fdata)
pbmc_cds <- new_cell_data_set(as(as.matrix(mat), 'sparseMatrix'),
cell_metadata = pdata,
gene_metadata = fdata)
# generate size factors for normalization later
#pbmc_cds <- estimateSizeFactors(pbmc_cds)#
pbmc_cds = preprocess_cds(pbmc_cds, num_dim = 100)
構(gòu)造標(biāo)記(marker )文件
除了表達(dá)數(shù)據(jù)之外，您需要的第二個(gè)主要輸入是一個(gè)標(biāo)記文件(marker file)。標(biāo)記文件包含以易于閱讀的文本格式編寫的單元類型定義列表。細(xì)胞類型定義告訴Garnett如何選擇細(xì)胞來(lái)訓(xùn)練模型。每個(gè)細(xì)胞類型定義以“>”符號(hào)和細(xì)胞類型名稱開(kāi)頭，后面是一系列帶有定義信息的行。定義行以關(guān)鍵字和':'開(kāi)頭，條目用逗號(hào)分隔。
一個(gè)簡(jiǎn)單有效的例子:
>B cells
expressed: CD19, MS4A1
>T cells
expressed: CD3D
有幾種方法可以在Garnett標(biāo)記文件格式中定義細(xì)胞類型。通常，每個(gè)細(xì)胞的定義可以包含三個(gè)主要組件。只需要第一個(gè)組件。
細(xì)胞類型的第一個(gè)也是最重要的規(guī)范是它的表達(dá)式。Garnett提供了幾種指定標(biāo)記基因的選項(xiàng)，詳情如下。
Format
expressed: gene1, gene2
not expressed: gene1, gene2
這是為你的細(xì)胞類型指定標(biāo)記基因的默認(rèn)方法。在使用這個(gè)規(guī)范時(shí)，Garnett為每個(gè)細(xì)胞計(jì)算一個(gè)標(biāo)記分?jǐn)?shù)，包括潛在的泄漏（leakage）、總體表達(dá)水平和read 深度。expressed:標(biāo)記應(yīng)該特定于所定義的細(xì)胞類型。
| Format | Example |
|---|---|
| expressed above: gene1 value, gene2 value | expressed above: MYOD1 4.2, MYH3 700 |
| expressed below: gene1 value, gene2 value | expressed below: PAX6 20, PAX3 4 |
| expressed between: gene1 value1 value2, gene2 value1 value2 | expressed between: PAX6 10 20, PAX3 4 100 |
這是指定expression的另一種方法，如果您知道某個(gè)基因的確切表達(dá)范圍，這種方法會(huì)很有用。但是，通常我們不建議使用這些規(guī)范，因?yàn)樗鼈儾粫?huì)考慮每個(gè)細(xì)胞中的read深度和總體表達(dá)。數(shù)據(jù)值與輸入數(shù)據(jù)的單位相同。
定義元數(shù)據(jù)
除了表達(dá)式信息之外，您還可以使用元數(shù)據(jù)進(jìn)一步細(xì)化細(xì)胞類型定義。這也是指定數(shù)據(jù)中所期望的任何子類型的地方。
| Format | Example |
|---|---|
| subtype of: celltype | subtype of: T cells |
| custom meta data: attribute1, attribute2 | tissue: spleen, thymus |
subtype of:允許您在定義文件中指定細(xì)胞類型是另一個(gè)細(xì)胞類型的子類型。
custom meta data:規(guī)范允許您為細(xì)胞類型提供任何進(jìn)一步的元數(shù)據(jù)需求。CDS對(duì)象的pData表中的任何列都可以用作元數(shù)據(jù)規(guī)范。在上面的示例中，pData表中有一個(gè)名為“tissue”的列。
提供你的證據(jù)
最后，我們強(qiáng)烈建議您記錄如何選擇標(biāo)記定義。為了便于跟蹤，我們提供了一個(gè)額外的規(guī)范—references:—它將存儲(chǔ)每種類型的引用信息。添加一組url或DOIs，它們將包含在分類器中。有關(guān)訪問(wèn)此信息的函數(shù)，請(qǐng)參見(jiàn)此處。
添加注釋
與R代碼類似，我們已經(jīng)包含了一個(gè)注釋字符#，這樣您就可以在您的標(biāo)記文件中添加注釋/注釋。任何在同一條線上的#之后的內(nèi)容都會(huì)被忽略。
所以一個(gè)完整的例子是這樣的：
>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=1200072
目前，標(biāo)記文件不能有回車(\r)，如果您在某些Windows文本編輯器中生成標(biāo)記文件，就會(huì)自動(dòng)包括回車。相反，您必須使用換行字符(\n)。有許多方法可以在Windows上替換回車，例如使用notepad++，您可以使用replace、choose extended和replace all \r with \n。我打算盡快以更自動(dòng)化的方式解決這個(gè)問(wèn)題。
檢查標(biāo)記
因?yàn)槎x標(biāo)記文件通常是過(guò)程中最困難的部分，所以Garnett提供了一些函數(shù)來(lái)檢查標(biāo)記是否能夠正常工作。相關(guān)的兩個(gè)函數(shù)是check_marker和plot_marker。check_marker生成關(guān)于標(biāo)記的信息表，plot_marker繪制最相關(guān)的信息。
除了包含的小數(shù)據(jù)集之外，我們還在包中包含了兩個(gè)示例標(biāo)記文件:pbmc_bad_markers.txt
在生成您的marker_check data.frame 之后，您可以使用我們內(nèi)建的繪圖函數(shù)plot_marker來(lái)查看結(jié)果。
library(org.Hs.eg.db)
marker_file_path <- system.file("extdata", "pbmc_bad_markers.txt",
package = "garnett")
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")
plot_markers(marker_check)
-
db:dbis a required argument for a Bioconductor AnnotationDb-class package used for converting gene IDs. For example, for humans use org.Hs.eg.db. See available packages at the Bioconductor website. Load your chosen db usinglibrary(db). If your species does not have an AnnotationDb-class package, see here. -
cds_gene_id_type: This argument tells Garnett the format of the gene IDs in your CDS object. It should be one of the values incolumns(db). The default is "ENSEMBL". -
marker_file_gene_id_type: Similarly to above, this argument tells Garnett the format of the gene IDs in your marker file.
這個(gè)標(biāo)記圖提供了一些關(guān)于所選標(biāo)記是否正確的關(guān)鍵信息。首先，紅色標(biāo)記“not in db”讓我們知道標(biāo)記ACTN在org.Hs.eg.db注釋中沒(méi)有作為“SYMBOL”出現(xiàn)。在本例中，它是一個(gè)打印錯(cuò)誤。接下來(lái)，x軸顯示每個(gè)標(biāo)記的模糊度評(píng)分—當(dāng)包含該標(biāo)記時(shí)，測(cè)量有多少個(gè)cell接受了模糊標(biāo)簽—在本例中，ACTB和PTPRC具有很高的模糊度，應(yīng)該排除。
check_marker輸出的值和plot_marker繪制的值是分類器可以選擇的cell 數(shù)量的估計(jì)值。然而，它使用啟發(fā)式快速找到候選細(xì)胞，并不能完全匹配標(biāo)記所選擇的細(xì)胞。請(qǐng)使用這些數(shù)字作為相對(duì)的度量，而不是訓(xùn)練集的絕對(duì)表示。
關(guān)于歧義分?jǐn)?shù)的進(jìn)一步說(shuō)明:歧義分?jǐn)?shù)是當(dāng)一個(gè)標(biāo)記被包含在標(biāo)記文件中時(shí)，一個(gè)標(biāo)記被標(biāo)記為“歧義”的cell 的分?jǐn)?shù)。然而，一個(gè)高模糊度的分?jǐn)?shù)并不一定意味著一個(gè)給定的標(biāo)記是不具體的。這可能意味著一個(gè)不同的標(biāo)記是罪魁禍?zhǔn)?，但該?biāo)記也提名了許多其他未標(biāo)記的細(xì)胞(高提名率)。在決定排除哪個(gè)標(biāo)記之前，請(qǐng)仔細(xì)查看歧義度高的標(biāo)記和最模糊的cell類型marker。
在制作標(biāo)記圖之后，您可能想要修改標(biāo)記文件。在我們的示例中，我們將刪除ACTN、ACTB和PTPRC以得到最終的“pbmc_test”。txt的標(biāo)記文件。
默認(rèn)情況下，Garnett 將基因id轉(zhuǎn)換成ENSEMBL用于分類器。這使得分類器是可移植的，這樣它們就可以對(duì)未來(lái)可能具有不同基因id的數(shù)據(jù)集進(jìn)行分類。然而，對(duì)于某些生物體，ENSEMBL id要么不可用，要么不常用。在這些情況下，可以在train_cell_classifier和check_marker中使用參數(shù)classifier_gene_id_type來(lái)指定不同的ID類型。您選擇的值將與分類器一起存儲(chǔ)，因此在對(duì)未來(lái)的數(shù)據(jù)集進(jìn)行分類時(shí)不需要再次指定它。
訓(xùn)練分類器
現(xiàn)在是訓(xùn)練分類器的時(shí)候了。參數(shù)應(yīng)該與check_marker的參數(shù)非常接近。下面我將從默認(rèn)值更改的一個(gè)參數(shù)是num_unknown參數(shù)。這告訴Garnett 應(yīng)該比較多少個(gè)外群細(xì)胞。默認(rèn)值是500，但是在這個(gè)只有很少cell的數(shù)據(jù)集中，我們需要更少的cell。
library(org.Hs.eg.db)
set.seed(260)
marker_file_path <- system.file("extdata", "pbmc_test.txt",
package = "garnett")
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")
head(pData(pbmc_cds))
> head(pData(pbmc_cds))
DataFrame with 6 rows and 4 columns
tsne_1 tsne_2
<numeric> <numeric>
AAGCACTGCACACA-1 3.8403149909359 12.0841914129204
GGCTCACTGGTCTA-1 9.97096226657347 3.50539308651821
AGCACTGATATCTC-1 3.45952940410281 4.93527280576176
ACACGTGATATTCC-1 1.74394947394641 7.78267061846286
ATATGCCTCTGCAA-1 5.78344829514223 8.55889827553495
TGACGAACCTATTC-1 10.7928530485958 10.5852739146963
Size_Factor FACS_type
<numeric> <factor>
AAGCACTGCACACA-1 0.559181445161514 B cells
GGCTCACTGGTCTA-1 0.515934033527584 B cells
AGCACTGATATCTC-1 0.698028398302026 B cells
ACACGTGATATTCC-1 0.815631008885519 B cells
ATATGCCTCTGCAA-1 1.11532798424345 B cells
TGACGAACCTATTC-1 0.649469901028841 B cells
運(yùn)行train_cell_classifier之后，類型為“garnett_classifier”的輸出對(duì)象包含對(duì)細(xì)胞進(jìn)行分類所需的所有信息。
查看分類基因
Garnett 分類是使用多項(xiàng)彈性-網(wǎng)絡(luò)回歸訓(xùn)練（multinomial elastic-net regression）。這意味著選擇某些基因作為區(qū)分細(xì)胞類型的相關(guān)基因。所選擇的基因可能是有趣的，所以Garnett 包含了一個(gè)訪問(wèn)所選擇基因的功能。注意:Garnett 沒(méi)有對(duì)輸入標(biāo)記進(jìn)行正則化，所以無(wú)論如何，它們都會(huì)被包含在分類器中。
我們用來(lái)查看相關(guān)基因的函數(shù)是get_feature_genes。參數(shù)是分類器，您想查看哪個(gè)節(jié)點(diǎn)(如果您的樹(shù)是分層的)—使用“root”作為頂部節(jié)點(diǎn)，使用父細(xì)胞類型名稱作為其他節(jié)點(diǎn)，使用db作為您的物種。如果設(shè)置convert_ids = TRUE，該函數(shù)將自動(dòng)將基因id轉(zhuǎn)換為符號(hào)。
feature_genes <- get_feature_genes(pbmc_classifier,
node = "root",
db = org.Hs.eg.db)
head(feature_genes)
B cells T cells Unknown
(Intercept) -4.44471924 2.590206989 1.85451225
ENSG00000187608 -0.10903325 -0.199658082 0.30869133
ENSG00000157873 0.06411921 -0.295797377 0.23167817
ENSG00000157191 -0.07677221 -0.155893113 0.23266532
ENSG00000173436 0.02649831 0.089393639 -0.11589195
ENSG00000117318 0.07400983 -0.009415156 -0.06459467
Viewing references
我們?cè)谏厦娼忉屃巳绾卧跇?biāo)記文件中包含關(guān)于如何選擇標(biāo)記的文檔。為了獲取這些信息—查看如何為已經(jīng)訓(xùn)練好的分類器選擇標(biāo)記—使用函數(shù)get_classifier_references。除了分類器之外，還有一個(gè)額外的可選參數(shù)，稱為cell_type。如果傳遞細(xì)胞類型的名稱，則僅打印該單元類型的引用，否則將全部打印。
get_classifier_references(pbmc_classifier )
$`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"
提交一個(gè)分類器
我們鼓勵(lì)你提交你的高質(zhì)量的分類器給我們，這樣我們可以使他們對(duì)社區(qū)可用。為此，打開(kāi)一個(gè)特刊并在Garnett github存儲(chǔ)庫(kù)中填寫表單。點(diǎn)擊這里，點(diǎn)擊“New issue”按鈕開(kāi)始!
