作者：苑曉梅
責(zé)編：SXY
單細(xì)胞測序?qū)υS多復(fù)雜組織重新進(jìn)行分解分析，打破了我們對細(xì)胞類型的固有認(rèn)知。通常情況下，研究人員首先通過無監(jiān)督聚類，獲得細(xì)胞簇，然后根據(jù)Marker基因手動注釋每個簇可能的細(xì)胞類型，或者應(yīng)用"label transfer"比對到已經(jīng)分型的數(shù)據(jù)確定自己研究的細(xì)胞類型 （這也是單細(xì)胞整合分析的一個關(guān)鍵點(diǎn)，具體關(guān)注我們的單細(xì)胞課程）。然而，對大型數(shù)據(jù)集根據(jù)Marker基因手動注釋一來比較繁瑣，難以擴(kuò)展，二來Marker基因具有不唯一性，人為確定有選擇困難癥 （單細(xì)胞分群后，怎么找到Marker基因定義每一類群？）。"Label transfer"的方法需要預(yù)先注釋的數(shù)據(jù)，容易受batch effects影響。
那么，就要敲黑板啦！
作者提出了cellassign(https://irrationone.github.io/cellassign/introduction-to-cellassign.html)-應(yīng)用概率模型綜合分析已知的細(xì)胞類型標(biāo)記基因，將單細(xì)胞RNA測序數(shù)據(jù)注釋為predefined or de novo cell types。該方法于2019年9月發(fā)表在Nature methods,原文是Probabilistic cell-type assignment of single-cell RNA-seq for tumor microenvironment profiling。
流程概覽
cellassign基于Marker基因信息將單細(xì)胞RNA測序獲得的細(xì)胞分型匹配到已知細(xì)胞類型。與其他從單細(xì)胞RNA-seq數(shù)據(jù)中確定細(xì)胞類型的方法不同，cellassign不需要已經(jīng)標(biāo)記的單細(xì)胞或bulk數(shù)據(jù) - 只需要知道每個給定的基因是否是某種細(xì)胞類型的marker就好，想獲得這些Marker基因，可以看下CellMarker數(shù)據(jù)庫等。
以下為workflow (用戶的輸入是子圖a的上半部分：標(biāo)準(zhǔn)化后的表達(dá)矩陣和Marker基因-細(xì)胞類型對照表，輸出是細(xì)胞歸屬的已知細(xì)胞類型或新細(xì)胞類型):
包安裝
cellassign的安裝需要依賴于 tensorflow （機(jī)器學(xué)習(xí)經(jīng)典包，莫煩Python機(jī)器學(xué)習(xí)）,可以根據(jù)以下步驟進(jìn)行安裝：
install.packages("tensorflow")library(tensorflow)install_tensorflow()
安裝cellassign
BiocManager::install('cellassign')
或者安裝開發(fā)版
devtools::install_github("Irrationone/cellassign")
具體使用
library(SingleCellExperiment)
library(cellassign)
首先需要構(gòu)建SingleCellExperiment對象，當(dāng)然我們做到這步一般是已經(jīng)成功構(gòu)建過了，如果沒有構(gòu)建，可以參考以下代碼：
讀入數(shù)據(jù) (Smartseq2)，讀入逗號或Tab鍵分隔的表達(dá)矩陣 (原始counts)，存儲為數(shù)據(jù)框，行是基因，列是細(xì)胞。
#讀入逗號分隔的count matrix，存儲為數(shù)據(jù)框，該數(shù)據(jù)是單細(xì)胞大名鼎鼎的小鼠全器官數(shù)據(jù)集中的一部分
dat = read.delim("FACS/Kidney-counts.csv", sep=",", header=TRUE)dat[1:5,1:5]
## X A14.MAA000545.3_8_M.1.1 E1.MAA000545.3_8_M.1.1
## 1 0610005C13Rik 0 0
## 2 0610007C21Rik 1 0
## 3 0610007L01Rik 0 0
## 4 0610007N19Rik 0 0
## 5 0610007P08Rik 0 0
## M4.MAA000545.3_8_M.1.1 O21.MAA000545.3_8_M.1.1
## 1 0 0
## 2 0 0
## 3 0 0
## 4 0 0
## 5 0 0
數(shù)據(jù)庫第一列是基因名字，把它移除作為列名字：
dim(dat)
rownames(dat) <- dat[,1]
dat <- dat[,-1]
構(gòu)建scater對象 （更多操作見 Hemberg-lab單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析（七）-導(dǎo)入10X和SmartSeq2數(shù)據(jù)Tabula Muris）
sceset <- SingleCellExperiment(assays = list(counts = as.matrix(dat)), colData=cell_anns)
# 如果做實(shí)驗(yàn)室記錄了細(xì)胞來源的小鼠、處理等信息，可以整理成表格，采用read.table讀入存儲到call_anns里面一起構(gòu)建scater對象`
# sceset <- SingleCellExperiment(assays = list(counts = as.matrix(dat)), colData=cell_anns)
# 查看對象
str(sceset)
作者使用一個由10個標(biāo)記基因和500個細(xì)胞組成的SingleCellExperiment作為示例 (如果自己還沒有數(shù)據(jù)，可以繼續(xù)用作者提供的數(shù)據(jù)操作)：
data(example_sce)
print(example_sce)
#> class: SingleCellExperiment
#> dim: 10 500
#> metadata(1): params
#> assays(6): BatchCellMeans BaseCellMeans ... TrueCounts counts
#> rownames(10): Gene186 Gene205 ... Gene949 Gene994
#> rowData names(6): Gene BaseGeneMean ... DEFacGroup1 DEFacGroup2
#> colnames(500): Cell1 Cell2 ... Cell499 Cell500
#> colData names(5): Cell Batch Group ExpLibSize EM_group
#> reducedDimNames(0):
#> spikeNames(0):
為了方便起見，在SingleCellExperiment的Group slot中注釋了真正的細(xì)胞類型 （這里是模擬的名字，Group1,2,3等）：
print(head(example_sce$Group))
#> [1] "Group1" "Group2" "Group2" "Group1" "Group1" "Group1"
關(guān)于標(biāo)記基因分組，還提供了一個細(xì)胞類型與基因的二元矩陣示例（example_rho），如果基因是給定細(xì)胞類型的marker，則標(biāo)記為1，否則為0：我們先從各種文獻(xiàn)、數(shù)據(jù)庫（比如CellMarker）或者直接從PanglaoDB得到先驗(yàn)信息，如
將它以列表的形式讀入
example_rho<- list(B_cell = c("Gene186", "Gene269", "Gene526", "Gene536", "Gene994"),
T_cell = c("Gene205", "Gene575", "Gene754", "Gene773", "Gene949"))
print(str(example_rho))
#> List of 2
#> $ B_cell: chr [1:5] "Gene186" "Gene269" "Gene526" "Gene536" ...
#> $ T_cell: chr [1:5] "Gene205" "Gene575" "Gene754" "Gene773" ...
#> NULL
然后用函數(shù)marker_list_to_mat將其轉(zhuǎn)換為二進(jìn)制標(biāo)記基因矩陣。
example_rho <- marker_list_to_mat(example_rho)
print(example_rho)
#> B_cell T_cell other
#> Gene186 1 0 0
#> Gene205 0 1 0
#> Gene269 1 0 0
#> Gene526 1 0 0
#> Gene536 1 0 0
#> Gene575 0 1 0
#> Gene754 0 1 0
#> Gene773 0 1 0
#> Gene949 0 1 0
#> Gene994 1 0 0
其中other表示非其中任一類型的細(xì)胞，也可用include_other=FALSE去掉該列。
表達(dá)矩陣標(biāo)準(zhǔn)化
cellassign識別的是scater對象example_sce的slots部分內(nèi)容，需要用戶提供量化因子用于表達(dá)矩陣的標(biāo)準(zhǔn)化。
example_sce已經(jīng)預(yù)先做過運(yùn)算，操作自己的數(shù)據(jù)時(shí)建議使用scran包的computeSumFactors計(jì)算，代碼如下。(什么？你做的差異基因方法不合適？中提供了其它的計(jì)算方法和計(jì)算原理)
同時(shí)由于用于cell assign分析的scater對象只是原始表達(dá)矩陣的一部分，標(biāo)準(zhǔn)化時(shí)建議用原始表達(dá)矩陣所有基因進(jìn)行標(biāo)準(zhǔn)化。
qclust <- quickCluster(sceset, min.size = 30)
sceset <- computeSumFactors(sceset, sizes = 15, clusters = qclust)
sceset <- normalize(sceset)
Scater對象篩選
cellassign函數(shù)需要的scater對象是單細(xì)胞實(shí)驗(yàn)或輸入的基因表達(dá)矩陣的子集，只包含example_rho中含有的Marker gene的行；如果應(yīng)用自己的數(shù)據(jù)，需要一步過濾 (example_sce測試數(shù)據(jù)已經(jīng)過濾過，故跳過)，代碼如下（注意：Marker gene中基因的命名規(guī)則與sceset中基因命名規(guī)則一致，比如都為ENSEMBL ID或都為Gene Symbol）：
# 對scater對象進(jìn)行篩選
sceset <- sceset[intersect(rownames(example_rho), rownames(sceset)),]
cellassign核心步驟
# cellassign object
# 獲取sizefactor
s <- sizeFactors(example_sce)
fit <- cellassign(exprs_obj = example_sce,
marker_gene_info = example_rho,
s = s,
learning_rate = 1e-2,
shrinkage = TRUE,
verbose = FALSE)
print(fit)
#> A cellassign fit for 500 cells, 10 genes, 2 cell types with 0 covariates
#> To access cell types, call celltypes(x)
#> To access cell type probabilities, call cellprobs(x)
使用celltypes函數(shù)最大似然法估計(jì)（MLE）細(xì)胞類型：
print(head(celltypes(fit)))
#> [1] "Group1" "Group2" "Group2" "Group1" "Group1" "Group1"
以及所有MLE參數(shù)使用mleparams：
print(str(mleparams(fit)))
#> List of 9
#> $ delta : num [1:10, 1:2] 2.32 0 2.5 2.9 2.89 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : chr [1:10] "Gene186" "Gene205" "Gene269" "Gene526" ...
#> .. ..$ : chr [1:2] "Group1" "Group2"
#> $ beta : num [1:10, 1] 0.487 -0.255 -1.016 1.195 1.476 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : chr [1:10] "Gene186" "Gene205" "Gene269" "Gene526" ...
#> .. ..$ : NULL
#> $ phi : num [1:500, 1:10, 1:2] 1.86 1.87 1.86 1.86 1.86 ...
#> $ gamma : num [1:500, 1:2] 1.00 1.56e-145 1.01e-43 1.00 1.00 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : NULL
#> .. ..$ : chr [1:2] "Group1" "Group2"
#> $ mu : num [1:500, 1:10, 1:2] 22.6 80.9 11.5 15.5 15.8 ...
#> $ a : num [1:10(1d)] 1.03 1.08 1.13 1.19 1.26 ...
#> $ theta : num [1:2(1d)] 0.472 0.528
#> ..- attr(*, "dimnames")=List of 1
#> .. ..$ : chr [1:2] "Group1" "Group2"
#> $ ld_mean: num 1
#> $ ld_var : num 0.531
#> NULL
我們還可以使用cellprobs函數(shù)可視化賦值的概率，該函數(shù)返回每個細(xì)胞和細(xì)胞類型的概率矩陣：
pheatmap::pheatmap(cellprobs(fit))
最后，由于這是模擬數(shù)據(jù)，我們可以檢查與真實(shí)組值的一致性：
print(table(example_sce$Group, celltypes(fit)))
#>
#> Group1 Group2
#> Group1 236 0
#> Group2 0 264
腫瘤微環(huán)境Marker基因示例
對于人腫瘤微環(huán)境中的常見細(xì)胞類型，cellassign包中提供了一組示例標(biāo)記。
標(biāo)記基因可用于以下細(xì)胞類型：
B cells
T cells
Cytotoxic T cells
Monocyte/Macrophage
Epithelial cells
Myofibroblasts
Vascular smooth muscle cells
Endothelial cells
這些可以通過調(diào)用進(jìn)行訪問:
data(example_TME_markers)
注意這是兩列marker的列表：
names(example_TME_markers)
#> [1] "symbol" "ensembl"
ensembl 包含基因符號：
lapply(head(example_TME_markers$ensembl, n = 4), head, n = 4)
#> $`B cells`
#> VIM MS4A1 CD79A PTPRC
#> "ENSG00000026025" "ENSG00000156738" "ENSG00000105369" "ENSG00000081237"
#>
#> $`T cells`
#> VIM CD2 CD3D CD3E
#> "ENSG00000026025" "ENSG00000116824" "ENSG00000167286" "ENSG00000198851"
#>
#> $`Cytotoxic T cells`
#> VIM CD2 CD3D CD3E
#> "ENSG00000026025" "ENSG00000116824" "ENSG00000167286" "ENSG00000198851"
#>
#> $`Monocyte/Macrophage`
#> VIM CD14 FCGR3A CD33
#> "ENSG00000026025" "ENSG00000170458" "ENSG00000203747" "ENSG00000105383"
要將這些與cellassign一起使用，我們可以通過單元格類型矩陣將它們轉(zhuǎn)換為二進(jìn)制標(biāo)記：
marker_mat <- marker_list_to_mat(example_TME_markers$ensembl)
marker_mat[1:3, 1:3]
#> B cells T cells Cytotoxic T cells
#> ENSG00000010610 0 1 0
#> ENSG00000026025 1 1 1
#> ENSG00000039068 0 0 0
重要的是：單細(xì)胞實(shí)驗(yàn)或輸入基因表達(dá)矩陣應(yīng)該是相應(yīng)的子集，以匹配標(biāo)記輸入矩陣的行，例如， 如果sce是具有ensembl ID作為rownames的SingleCellExperiment，則調(diào)用：
# 對scater對象進(jìn)行篩選
sce_marker <- sce[intersect(rownames(marker_mat), rownames(sce)),]
局限性
1.CellAssign適用于標(biāo)記基因已經(jīng)非常明確的條件下，未知的細(xì)胞類型或細(xì)胞狀態(tài)可能是不可見的；2.作者在兩種不同細(xì)胞類型中的相同標(biāo)記的中等或高表達(dá)之間沒有先驗(yàn)區(qū)分。
更多單細(xì)胞操作
Science：通過單細(xì)胞轉(zhuǎn)錄組測序揭示玉米減數(shù)分裂進(jìn)程 | 很好的單細(xì)胞分析案例
Nature系列 整合單細(xì)胞轉(zhuǎn)錄組學(xué)和質(zhì)譜流式確定類風(fēng)濕性關(guān)節(jié)炎滑膜組織中的炎癥細(xì)胞狀態(tài) 詳細(xì)解讀
Hemberg-lab單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析（二）- 實(shí)驗(yàn)平臺
Hemberg-lab單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析（三）- 原始數(shù)據(jù)質(zhì)控
Hemberg-lab單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析（五）- STAR, Kallisto定量
Hemberg-lab單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析（六）- 構(gòu)建表達(dá)矩陣，UMI介紹
Hemberg-lab單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析（七）- 導(dǎo)入10X和SmartSeq2數(shù)據(jù)Tabula Muris
Hemberg-lab單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析（八）- Scater包輸入導(dǎo)入和存儲
Hemberg-lab單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析（九）- Scater包單細(xì)胞過濾
Hemberg-lab單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析（十一）- Scater單細(xì)胞表達(dá)譜PCA可視化
Hemberg-lab單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析（十二）- Scater單細(xì)胞表達(dá)譜tSNE可視化
植物單細(xì)胞轉(zhuǎn)錄組的春天來了，還不上車？Science, PC, PP, MP, bioRxiv各一個
運(yùn)行環(huán)境
sessionInfo()
#> R version 3.6.0 (2019-04-26)
#> Platform: x86_64-apple-darwin15.6.0 (64-bit)
#> Running under: macOS Mojave 10.14
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
#>
#> locale:
#> [1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8
#>
#> attached base packages:
#> [1] parallel stats4 stats graphics grDevices utils datasets
#> [8] methods base
#>
#> other attached packages:
#> [1] SingleCellExperiment_1.6.0 SummarizedExperiment_1.14.0
#> [3] DelayedArray_0.10.0 BiocParallel_1.18.0
#> [5] GenomicRanges_1.36.0 GenomeInfoDb_1.20.0
#> [7] cellassign_0.99.11 pheatmap_1.0.12
#> [9] matrixStats_0.54.0 edgeR_3.26.5
#> [11] org.Hs.eg.db_3.8.2 AnnotationDbi_1.46.0
#> [13] IRanges_2.18.1 S4Vectors_0.22.0
#> [15] Biobase_2.44.0 BiocGenerics_0.30.0
#> [17] limma_3.40.2 magrittr_1.5
#> [19] BiocStyle_2.12.0
#>
#> loaded via a namespace (and not attached):
#> [1] pkgload_1.0.2 bit64_0.9-7 jsonlite_1.6
#> [4] assertthat_0.2.1 BiocManager_1.30.4 blob_1.1.1
#> [7] GenomeInfoDbData_1.2.1 yaml_2.2.0 remotes_2.1.0
#> [10] sessioninfo_1.1.1 pillar_1.4.2 RSQLite_2.1.1
#> [13] backports_1.1.4 lattice_0.20-38 glue_1.3.1
#> [16] reticulate_1.12 digest_0.6.20 RColorBrewer_1.1-2
#> [19] XVector_0.24.0 colorspace_1.4-1 plyr_1.8.4
#> [22] htmltools_0.3.6 Matrix_1.2-17 pkgconfig_2.0.2
#> [25] devtools_2.0.2 bookdown_0.11 zlibbioc_1.30.0
#> [28] purrr_0.3.2 scales_1.0.0 processx_3.4.0
#> [31] whisker_0.3-2 tibble_2.1.3 usethis_1.5.1
#> [34] withr_2.1.2 cli_1.1.0 crayon_1.3.4
#> [37] memoise_1.1.0 evaluate_0.14 ps_1.3.0
#> [40] fs_1.3.1 pkgbuild_1.0.3 tools_3.6.0
#> [43] prettyunits_1.0.2 stringr_1.4.0 munsell_0.5.0
#> [46] locfit_1.5-9.1 callr_3.3.0 compiler_3.6.0
#> [49] rlang_0.4.0 grid_3.6.0 RCurl_1.95-4.12
#> [52] rstudioapi_0.10 bitops_1.0-6 base64enc_0.1-3
#> [55] rmarkdown_1.13 testthat_2.1.1 gtable_0.3.0
#> [58] DBI_1.0.0 R6_2.4.0 tfruns_1.4
#> [61] dplyr_0.8.3 knitr_1.23 tensorflow_1.13.1
#> [64] bit_1.1-14 rprojroot_1.3-2 desc_1.2.0
#> [67] stringi_1.4.3 Rcpp_1.0.1 tidyselect_0.2.5
#> [70] xfun_0.8
對單細(xì)胞轉(zhuǎn)錄組感興趣的話，不妨留意一下我們將在北京召開的單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)整合分析專題（11月29號-12月01號）：
