詳細(xì)介紹新版SingleR進(jìn)行單細(xì)胞自動(dòng)化鑒定

作者：堯小飛
審稿：童蒙
編輯：amethyst

引言

隨著10xGenomics單細(xì)胞轉(zhuǎn)錄組技術(shù)的普及堵腹，現(xiàn)在的單細(xì)胞轉(zhuǎn)錄組分析越來(lái)越普及星澳，得到的單細(xì)胞數(shù)據(jù)越來(lái)越多疚顷。然而單細(xì)胞轉(zhuǎn)錄組技術(shù)發(fā)展年限相對(duì)而言較短禁偎，在分析上還存在著一定的瓶頸，其中比較頭疼的一個(gè)瓶頸就是單細(xì)胞亞群的鑒定如暖。比如說(shuō)：得到這么多細(xì)胞，我不知道細(xì)胞是什么細(xì)胞盒至，這怎么分析酗洒，這是單細(xì)胞轉(zhuǎn)錄組分析中繞不過(guò)去的一點(diǎn)枷遂。

一般來(lái)說(shuō)，需要老師根據(jù)自己專業(yè)背景知識(shí)登淘，進(jìn)行判斷單細(xì)胞亞群是什么類型細(xì)胞箫老，但不是每個(gè)老師的專業(yè)背景都這么強(qiáng)黔州。

這個(gè)問(wèn)題并不是我們才有的問(wèn)題，全世界都有這個(gè)問(wèn)題流妻。幸好Dvir等人在2019年1月14日發(fā)表在Nature Immunology （https://www.x-mol.com/paper/journal/104）( IF 23.530 )上的文章解決了這個(gè)問(wèn)題，作者開(kāi)發(fā)了一個(gè)SingleR （http://www.bioconductor.org/packages/release/bioc/vignettes/SingleR/inst/doc/SingleR.html）工具绅这，可以自動(dòng)鑒定細(xì)胞亞群的類型，這個(gè)工具支持人和小鼠兩個(gè)物種，其實(shí)工具的本質(zhì)上計(jì)算每個(gè)單細(xì)胞的表達(dá)量數(shù)據(jù)與純的bulk RNA數(shù)據(jù)相關(guān)性度苔，相關(guān)性越大，說(shuō)明是此類細(xì)胞寇窑，當(dāng)然作者對(duì)數(shù)據(jù)有其他的處理過(guò)程。SingleR也有在線版本的甩骏。主要功能有如下：

for each test cell:

We compute the Spearman correlation between its expression profile and that of each reference sample. This is done across the union of marker genes identified between all pairs of labels.

We define the per-label score as a fixed quantile (by default, 0.8) of the distribution of correlations.

We repeat this for all labels and we take the label with the highest score as the annotation for this cell.

We optionally perform a fine-tuning step:

The reference dataset is subsetted to only include labels with scores close to the maximum.

Scores are recomputed using only marker genes for the subset of labels.

This is iterated until one label remains.

下面我們來(lái)詳細(xì)看看如何使用這個(gè)軟件。

SingleR的安裝

使用bioconductor安裝
SingleR的安裝有多種方式饮笛，這是一個(gè)R包，并且發(fā)布在biocondutor福青，當(dāng)然可以直接用biocondutor方式來(lái)安裝：

if (!requireNamespace("BiocManager", quietly = TRUE))
 install.packages("BiocManager")
BiocManager::install("SingleR")

使用github來(lái)安裝
這個(gè)R包已經(jīng)放在GitHub上面，可以直接使用GitHub安裝R包的方式來(lái)安裝：

devtools::install_github('dviraran/SingleR')

不過(guò)在查找問(wèn)題的時(shí)候素跺，發(fā)現(xiàn)這個(gè)包已經(jīng)放在別的倉(cāng)庫(kù)了，因此可以通過(guò)別的倉(cāng)庫(kù)來(lái)安裝：

devtools::install_github('LTLA/SingleR')

注意事項(xiàng)
新的R包指厌，需要在 R (version "3.6")版本以上的才能安裝，其實(shí)一直很納悶踩验，為什么R一升級(jí)，就不能兼容下面的箕憾，還得重新安裝R牡借，這是R的一個(gè)大弊病袭异。

SingleR升級(jí)

一般來(lái)說(shuō)，如果工具能夠一直維護(hù)和升級(jí)御铃，當(dāng)然是好事碴里，singleR第一版出來(lái)以后上真，解決了沒(méi)有工具自動(dòng)進(jìn)行細(xì)胞亞群鑒定的問(wèn)題，但是存在一個(gè)較大的問(wèn)題睡互，就是這個(gè)運(yùn)行速度特別的慢陵像，細(xì)胞類型少，自己定制參考數(shù)據(jù)集較困難寇壳，因此亟需改進(jìn)。作者在1月份發(fā)表文章九巡，5月23日就做了一個(gè)較大的升級(jí)蹂季。這次的升級(jí)與之前的分析很多不兼容，具體升級(jí)的功能有：

1.官方提供了部分?jǐn)?shù)據(jù)集偿洁，細(xì)胞類型增多撒汉，免疫細(xì)胞分的更細(xì)了涕滋，更符合分析預(yù)期；
2.提高了運(yùn)行速度宾肺；
3.修改了部分函數(shù)，舊版本的函數(shù)與新版本的函數(shù)可能不兼容锨用；
4.提供了參考數(shù)據(jù)集的接口丰刊，雖然之前也提供了增拥，但是現(xiàn)在的接口使用更便啄巧。

SingleR參考數(shù)據(jù)集

新版的SingleR一大優(yōu)點(diǎn)就是提供了較多的數(shù)據(jù)集掌栅，不是之前的所有數(shù)據(jù)集整合在里面，我們看不見(jiàn)摸不著猾封，而且細(xì)胞類型較少，新的數(shù)據(jù)集有7種晌缘，主要是人和小鼠物種，特別的兩個(gè)物種的免疫基因枚钓，可以把細(xì)胞分成很細(xì)的T細(xì)胞铅搓。具體如下：

每個(gè)數(shù)據(jù)集類型較多搀捷，就不一一列舉了多望。以下主要列舉兩個(gè)免疫基因數(shù)據(jù)集，主要是在單細(xì)胞轉(zhuǎn)錄組研究中怀偷，免疫細(xì)胞的研究較多，而且免疫細(xì)胞亞群鑒定相對(duì)困難椎工，下面的數(shù)據(jù)集是非常好的參考數(shù)據(jù)集。

01 DatabaseImmuneCellExpressionData Labels

02 ImmGenData Labels

當(dāng)然下面的表格不全维蒙，只是列舉了少部分，從這些細(xì)胞類型來(lái)看果覆，很多都能滿足我們常見(jiàn)的免疫細(xì)胞鑒定。

03 其他數(shù)據(jù)集的細(xì)胞類型

當(dāng)然下面的表格不全局待，只是列舉了少部分，從這些細(xì)胞類型來(lái)看钳榨，很多都能滿足我們常見(jiàn)的免疫細(xì)胞鑒定，其他數(shù)據(jù)集可以看SingleR的網(wǎng)頁(yè)薛耻，點(diǎn)擊如下的紅色箭頭指定的部分就可以看到所有細(xì)胞類型。

SingleR測(cè)試

報(bào)錯(cuò)的解決

在測(cè)試SingleR的時(shí)候昭卓，首先就是想把參考數(shù)據(jù)集下載下來(lái)，以后做分析的時(shí)候候醒，可以直接用，不用聯(lián)網(wǎng)分析倒淫。結(jié)果下載的時(shí)候，有如下所示報(bào)錯(cuò)敌土，說(shuō)明格式有問(wèn)題镜硕。

然后發(fā)現(xiàn)數(shù)據(jù)下載是通過(guò)ExperimentHub下載返干，其下載源代碼如下：

 1.create_se <- function(dataset, hub = ExperimentHub(), assays="logcounts",
 2    rm.NA = c("rows","cols","both","none"), has.rowdata=FALSE, has.coldata=TRUE) 
 3{
 4    rm.NA <- match.arg(rm.NA)
 5    host <- file.path("SingleR", dataset)
 6    # hub$rdatapath[grep('SingleR',hub$rdatapath)]
 7    ## extract normalized values --------
 8    all.assays <- list()
 9    for (a in assays) {
10        nrmcnts <- hub[hub$rdatapath==file.path(host, paste0(a, ".rds"))][[1]]
11        all.assays[[a]] <- .rm_NAs(nrmcnts, rm.NA)
12    }
13
14    ## get metadata ----------------------
15    args <- list()
16    if (has.coldata) {
17        args$colData <- hub[hub$rdatapath==file.path(host, "coldata.rds")][[1]]
18    }
19    if (has.rowdata) {
20        args$rowData <- hub[hub$rdatapath==file.path(host, "rowdata.rds")][[1]]
21    }
22
23    ## make the final SE object ----------
24    do.call(SummarizedExperiment, c(list(assays=all.assays), args))
25}
26##下載函數(shù)
27ImmGenData <- function(ensembl=FALSE) {
28    version <- "1.0.0"
29    se <- .create_se(file.path("immgen", version), 
30        assays="logcounts", rm.NA = "none",
31        has.rowdata = FALSE, has.coldata = TRUE)
32
33    if (ensembl) {
34        se <- .convert_to_ensembl(se, "Mm")
35    }
36    se
37}
38 hub$rdatapath[grep('SingleR',hub$rdatapath)]
39> hub$rdatapath[grep('SingleR',hub$rdatapath)]
40 #[1] "SingleR/blueprint_encode/1.0.0/logcounts.rds"
41 #[2] "SingleR/blueprint_encode/1.0.0/coldata.rds"
42 #[3] "SingleR/hpca/1.0.0/logcounts.rds"
43 #[4] "SingleR/hpca/1.0.0/coldata.rds"
44 #[5] "SingleR/immgen/1.0.0/logcounts.rds"
45 #[6] "SingleR/immgen/1.0.0/coldata.rds"
46 #[7] "SingleR/mouse.rnaseq/1.0.0/logcounts.rds"
47 #[8] "SingleR/mouse.rnaseq/1.0.0/coldata.rds"
48# [9] "SingleR/dice/1.0.0/logcounts.rds"
49#[10] "SingleR/dice/1.0.0/coldata.rds"
50#[11] "SingleR/dmap/1.0.0/logcounts.rds"
51#[12] "SingleR/dmap/1.0.0/coldata.rds"
52#[13] "SingleR/monaco_immune/1.0.0/logcounts.rds"
53#[14] "SingleR/monaco_immune/1.0.0/coldata.rds"

由于運(yùn)行一直有報(bào)錯(cuò)，然后就一步一步地運(yùn)行财剖，在如下代碼的時(shí)候報(bào)錯(cuò)：

1args$colData <- hub[hub$rdatapath==file.path(host, "coldata.rds")][[1]]
2h=hub[hub$rdatapath==file.path(host, paste0(a, ".rds"))]
3h[[1]]

上面圖片的報(bào)錯(cuò)悠夯，可能還是網(wǎng)絡(luò)的問(wèn)題躺坟，沒(méi)法翻墻的情況下就不能用這個(gè)方法沦补，只能繼續(xù)尋找其他方法咪橙。之后在GitHub上面發(fā)現(xiàn)有如下的數(shù)據(jù)，鑒于無(wú)法翻墻的網(wǎng)絡(luò)美侦，還是下載不了产舞，不過(guò)在GitHub上面就好辦音榜，可以先把這個(gè)項(xiàng)目導(dǎo)入到gitee捧弃，然后從gitee下載赠叼。

通過(guò)下面方式導(dǎo)入违霞，這個(gè)工具可能有點(diǎn)大嘴办，800M买鸽，估計(jì)要點(diǎn)時(shí)間涧郊。

結(jié)果示范圖

不同參考數(shù)據(jù)集的影響

01 運(yùn)行速度-不同數(shù)據(jù)集的影響

為了測(cè)試不同的參考數(shù)據(jù)集對(duì)結(jié)果的影響眼五，這里用了同一個(gè)querry數(shù)據(jù)集、不同的參考數(shù)據(jù)集進(jìn)行測(cè)試看幼。這里選擇的測(cè)試數(shù)據(jù)為小鼠，querry數(shù)據(jù)集是經(jīng)過(guò)流式sort后的T細(xì)胞诵姜，有20519個(gè)細(xì)胞汽煮，這屬于常規(guī)項(xiàng)目棚唆，細(xì)胞數(shù)不少。小鼠的參考數(shù)據(jù)集就只有兩個(gè)宵凌，ImmGenData()和MouseRNAseqData()鞋囊，前一個(gè)數(shù)據(jù)集主要集中在免疫細(xì)胞相關(guān)瞎惫，后一個(gè)的細(xì)胞范圍較廣溜腐，兩個(gè)數(shù)據(jù)集都有T細(xì)胞，因此兩個(gè)數(shù)據(jù)集都可以使用逗扒。

下面表格為測(cè)試結(jié)果，從結(jié)果來(lái)看矩肩，參考數(shù)據(jù)集的細(xì)胞類型越少，其運(yùn)行時(shí)間越短黍檩，MouseRNAseqData參考數(shù)據(jù)集所耗時(shí)間只有3 min左右。但是main labels為20個(gè)的ImmGenData所費(fèi)時(shí)間要遠(yuǎn)遠(yuǎn)高于MouseRNAseqData的時(shí)間刽酱，這說(shuō)明運(yùn)行時(shí)間與參考數(shù)據(jù)集的細(xì)胞數(shù)或者說(shuō)樣品數(shù)目有關(guān)喳逛，ImmGenData有830個(gè)細(xì)胞或者樣品棵里，而MouseRNAseqData只有358個(gè)細(xì)胞或者樣品。

02 結(jié)果準(zhǔn)確性-不同數(shù)據(jù)集的影響

其實(shí)在這里說(shuō)結(jié)果準(zhǔn)確性有點(diǎn)不合適殿怜，我們也不知道具體的結(jié)果如何，只能比較一下其結(jié)果头谜。

1. MouseRNAseqData

a .main labels

> data<-immune.combined@assays$RNA@data
> pred.hesc <- SingleR(test = data, ref =ref_data, labels = ref_data$label.main)
> table(pred.hesc$labels)
Dendritic cells    Granulocytes        NK cells         T cells 
              1               1             192           20325

b. fine labels

> pred.hesc_label.fine <- SingleR(test = data, ref =ref_data, labels = ref_data$label.fine)
> table(pred.hesc_label.fine$labels)

Granulocytes     NK cells      T cells 
           1          189        20329

從上面兩種標(biāo)簽的結(jié)果來(lái)看骏掀，其結(jié)果基本上沒(méi)有太大差異柱告，不過(guò)有部分不一致，這個(gè)數(shù)據(jù)集是經(jīng)過(guò)sort以后的T細(xì)胞际度，因此大部分細(xì)胞應(yīng)該屬于T細(xì)胞葵袭，從這個(gè)來(lái)說(shuō)甲脏，上述結(jié)果都比較符合預(yù)期眶熬，比較sort的細(xì)胞也不能保證百分比的純細(xì)胞块请，其他細(xì)胞1%以下還是可以接受的。如果非要分個(gè)上學(xué)墩新，fine labels的T細(xì)胞數(shù)目更多贸弥，不過(guò)差別不大海渊，只多了4個(gè)細(xì)胞绵疲，幾乎可以忽略不計(jì)。

2. ImmGenData

ImmGenData數(shù)據(jù)的細(xì)胞數(shù)目和細(xì)胞種類要遠(yuǎn)遠(yuǎn)大于MouseRNAseqData數(shù)據(jù)集盔憨，因此運(yùn)行的速度會(huì)遠(yuǎn)遠(yuǎn)低于它，但是細(xì)胞類型多郁岩，各種免疫細(xì)胞都有婿奔，其結(jié)果將更符合研究目的问慎。
a. main labels

> pred.hesc <- SingleR(test = data, ref =ref_data, labels = ref_data$label.main)
> table(pred.hesc$labels)
Basophils       ILC  NK cells       NKT   T cells       Tgd 
        2        81        39       721     19435       241