作者：堯小飛
審稿：童蒙
編輯：amethyst
引言
隨著10xGenomics單細(xì)胞轉(zhuǎn)錄組技術(shù)的普及，現(xiàn)在的單細(xì)胞轉(zhuǎn)錄組分析越來越普及，得到的單細(xì)胞數(shù)據(jù)越來越多。然而單細(xì)胞轉(zhuǎn)錄組技術(shù)發(fā)展年限相對(duì)而言較短，在分析上還存在著一定的瓶頸，其中比較頭疼的一個(gè)瓶頸就是單細(xì)胞亞群的鑒定。比如說：得到這么多細(xì)胞，我不知道細(xì)胞是什么細(xì)胞，這怎么分析，這是單細(xì)胞轉(zhuǎn)錄組分析中繞不過去的一點(diǎn)。
一般來說，需要老師根據(jù)自己專業(yè)背景知識(shí)，進(jìn)行判斷單細(xì)胞亞群是什么類型細(xì)胞，但不是每個(gè)老師的專業(yè)背景都這么強(qiáng)。
這個(gè)問題并不是我們才有的問題，全世界都有這個(gè)問題。幸好Dvir等人在2019年1月14日發(fā)表在Nature Immunology （https://www.x-mol.com/paper/journal/104）( IF 23.530 )上的文章解決了這個(gè)問題，作者開發(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)性越大，說明是此類細(xì)胞，當(dāng)然作者對(duì)數(shù)據(jù)有其他的處理過程。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.
下面我們來詳細(xì)看看如何使用這個(gè)軟件。
SingleR的安裝
- 使用bioconductor安裝
SingleR的安裝有多種方式，這是一個(gè)R包，并且發(fā)布在biocondutor，當(dāng)然可以直接用biocondutor方式來安裝：
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("SingleR")
- 使用github來安裝
這個(gè)R包已經(jīng)放在GitHub上面，可以直接使用GitHub安裝R包的方式來安裝：
devtools::install_github('dviraran/SingleR')
不過在查找問題的時(shí)候，發(fā)現(xiàn)這個(gè)包已經(jīng)放在別的倉(cāng)庫(kù)了，因此可以通過別的倉(cāng)庫(kù)來安裝：
devtools::install_github('LTLA/SingleR')
- 注意事項(xiàng)
新的R包，需要在 R (version "3.6")版本以上的才能安裝，其實(shí)一直很納悶，為什么R一升級(jí)，就不能兼容下面的，還得重新安裝R，這是R的一個(gè)大弊病。
SingleR升級(jí)
一般來說，如果工具能夠一直維護(hù)和升級(jí)，當(dāng)然是好事，singleR第一版出來以后，解決了沒有工具自動(dòng)進(jìn)行細(xì)胞亞群鑒定的問題，但是存在一個(gè)較大的問題，就是這個(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ù)集整合在里面，我們看不見摸不著，而且細(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ì)胞類型來看，很多都能滿足我們常見的免疫細(xì)胞鑒定。
03 其他數(shù)據(jù)集的細(xì)胞類型
當(dāng)然下面的表格不全，只是列舉了少部分，從這些細(xì)胞類型來看，很多都能滿足我們常見的免疫細(xì)胞鑒定，其他數(shù)據(jù)集可以看SingleR的網(wǎng)頁(yè)，點(diǎn)擊如下的紅色箭頭指定的部分就可以看到所有細(xì)胞類型。
SingleR測(cè)試
報(bào)錯(cuò)的解決
在測(cè)試SingleR的時(shí)候，首先就是想把參考數(shù)據(jù)集下載下來，以后做分析的時(shí)候，可以直接用，不用聯(lián)網(wǎng)分析。結(jié)果下載的時(shí)候，有如下所示報(bào)錯(cuò)，說明格式有問題。
然后發(fā)現(xiàn)數(shù)據(jù)下載是通過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ò)的問題，沒法翻墻的情況下就不能用這個(gè)方法，只能繼續(xù)尋找其他方法。之后在GitHub上面發(fā)現(xiàn)有如下的數(shù)據(jù)，鑒于無法翻墻的網(wǎng)絡(luò)，還是下載不了，不過在GitHub上面就好辦，可以先把這個(gè)項(xiàng)目導(dǎo)入到gitee，然后從gitee下載。
通過下面方式導(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)過流式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é)果來看，參考數(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í)間，這說明運(yùn)行時(shí)間與參考數(shù)據(jù)集的細(xì)胞數(shù)或者說樣品數(shù)目有關(guān)，ImmGenData有830個(gè)細(xì)胞或者樣品，而MouseRNAseqData只有358個(gè)細(xì)胞或者樣品。
02 結(jié)果準(zhǔn)確性-不同數(shù)據(jù)集的影響
其實(shí)在這里說結(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é)果來看，其結(jié)果基本上沒有太大差異，不過有部分不一致，這個(gè)數(shù)據(jù)集是經(jīng)過sort以后的T細(xì)胞，因此大部分細(xì)胞應(yīng)該屬于T細(xì)胞，從這個(gè)來說，上述結(jié)果都比較符合預(yù)期，比較sort的細(xì)胞也不能保證百分比的純細(xì)胞，其他細(xì)胞1%以下還是可以接受的。如果非要分個(gè)上學(xué)，fine labels的T細(xì)胞數(shù)目更多，不過差別不大，只多了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
b. fine labels
1> pred.hesc_label.fine <- SingleR(test = data, ref =ref_data, labels = ref_data$label.fine)
2> as.matrix(table(pred.hesc_label.fine$labels))
3T cells (T.4SP69+) 1
4T cells (T.8EFF.OT1.12HR.LISOVA) 114
5T cells (T.8EFF.OT1.24HR.LISOVA) 135
6T cells (T.8EFF.OT1.48HR.LISOVA) 66
7T cells (T.8Mem) 1515
8T cells (T.8MEM) 61
9T cells (T.8MEMKLRG1-CD127+.D8.LISOVA) 445
10T cells (T.8NVE.OT1) 386
11T cells (T.8Nve) 574
12T cells (T.8NVE) 19
13T cells (T.CD4.1H) 1707
14T cells (T.CD4.24H) 6
15T cells (T.CD4.48H) 1
16T cells (T.CD4.5H) 40
17T cells (T.CD4.CTR) 22
18T cells (T.CD4+TESTDB) 1
統(tǒng)計(jì)了一下，T cells 一共有19761個(gè)細(xì)胞，比main labels的19435多了300多個(gè)細(xì)胞。其他的結(jié)果就不統(tǒng)計(jì)了。上面結(jié)果只展示了部分的細(xì)胞結(jié)果，如果僅僅從這個(gè)結(jié)果來看，這個(gè)參數(shù)數(shù)據(jù)集可以將T細(xì)胞或者其他免疫細(xì)胞分的很細(xì)，比如這里可以分到naive細(xì)胞或者效應(yīng)細(xì)胞等等。因此這個(gè)數(shù)據(jù)集適合哪種sort T細(xì)胞或者其他免疫細(xì)胞。
結(jié)語(yǔ)
以上就是SingleR研究過程中的一些心得，希望可以對(duì)各位小伙伴有所啟發(fā)~~
