作者：苑晓梅
责编：SXY

单细胞测序对许多复杂组织重新进行分解分析，打破了我们对细胞类型的固有认知。通常情况下，研究人员首先通过无监督聚类，获得细胞簇，然后根据Marker基因手动注释每个簇可能的细胞类型，或者应用"label transfer"比对到已经分型的数据确定自己研究的细胞类型 （这也是单细胞整合分析的一个关键点，具体关注我们的单细胞课程）。然而，对大型数据集根据Marker基因手动注释一来比较繁琐，难以扩展，二来Marker基因具有不唯一性，人为确定有选择困难症 （单细胞分群后，怎么找到Marker基因定义每一类群？）。"Label transfer"的方法需要预先注释的数据，容易受batch effects影响。

那么，就要敲黑板啦！

作者提出了cellassign(https://irrationone.github.io/cellassign/introduction-to-cellassign.html)-应用概率模型综合分析已知的细胞类型标记基因，将单细胞RNA测序数据注释为predefined or de novo cell types。该方法于2019年9月发表在Nature methods,原文是Probabilistic cell-type assignment of single-cell RNA-seq for tumor microenvironment profiling。

流程概览

cellassign基于Marker基因信息将单细胞RNA测序获得的细胞分型匹配到已知细胞类型。与其他从单细胞RNA-seq数据中确定细胞类型的方法不同，cellassign不需要已经标记的单细胞或bulk数据 - 只需要知道每个给定的基因是否是某种细胞类型的marker就好，想获得这些Marker基因，可以看下CellMarker数据库等。

以下为workflow (用户的输入是子图a的上半部分：标准化后的表达矩阵和Marker基因-细胞类型对照表，输出是细胞归属的已知细胞类型或新细胞类型):

包安装

cellassign的安装需要依赖于 tensorflow （机器学习经典包，莫烦Python机器学习）,可以根据以下步骤进行安装：

install.packages("tensorflow")
library(tensorflow)
install_tensorflow()

安装cellassign

BiocManager::install('cellassign')

或者安装开发版

devtools::install_github("Irrationone/cellassign")

具体使用

library(SingleCellExperiment)
library(cellassign)

首先需要构建SingleCellExperiment对象，当然我们做到这步一般是已经成功构建过了，如果没有构建，可以参考以下代码：

读入数据 (Smartseq2)，读入逗号或Tab键分隔的表达矩阵 (原始counts)，存储为数据框，行是基因，列是细胞。

#读入逗号分隔的count matrix，存储为数据框，该数据是单细胞大名鼎鼎的小鼠全器官数据集中的一部分
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

数据库第一列是基因名字，把它移除作为列名字：

dim(dat)
rownames(dat) <- dat[,1]
dat <- dat[,-1]

构建scater对象 （更多操作见 Hemberg-lab单细胞转录组数据分析（七）-导入10X和SmartSeq2数据Tabula Muris）

sceset <- SingleCellExperiment(assays = list(counts = as.matrix(dat)), colData=cell_anns)

# 如果做实验室记录了细胞来源的小鼠、处理等信息，可以整理成表格，采用read.table读入存储到call_anns里面一起构建scater对象`
# sceset <- SingleCellExperiment(assays = list(counts = as.matrix(dat)), colData=cell_anns)
# 查看对象
str(sceset)

作者使用一个由10个标记基因和500个细胞组成的SingleCellExperiment作为示例 (如果自己还没有数据，可以继续用作者提供的数据操作)：

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中注释了真正的细胞类型 （这里是模拟的名字，Group1,2,3等）：

print(head(example_sce$Group))
#> [1] "Group1" "Group2" "Group2" "Group1" "Group1" "Group1"

关于标记基因分组，还提供了一个细胞类型与基因的二元矩阵示例（example_rho），如果基因是给定细胞类型的marker，则标记为1，否则为0：我们先从各种文献、数据库（比如CellMarker）或者直接从PanglaoDB得到先验信息，如

将它以列表的形式读入

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

然后用函数marker_list_to_mat将其转换为二进制标记基因矩阵。

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表示非其中任一类型的细胞，也可用include_other=FALSE去掉该列。

表达矩阵标准化

cellassign识别的是scater对象example_sce的slots部分内容，需要用户提供量化因子用于表达矩阵的标准化。

example_sce已经预先做过运算，操作自己的数据时建议使用scran包的computeSumFactors计算，代码如下。(什么？你做的差异基因方法不合适？中提供了其它的计算方法和计算原理)

同时由于用于cell assign分析的scater对象只是原始表达矩阵的一部分，标准化时建议用原始表达矩阵所有基因进行标准化。

qclust <- quickCluster(sceset, min.size = 30)
sceset <- computeSumFactors(sceset, sizes = 15, clusters = qclust)
sceset <- normalize(sceset)

Scater对象筛选

cellassign函数需要的scater对象是单细胞实验或输入的基因表达矩阵的子集，只包含example_rho中含有的Marker gene的行；如果应用自己的数据，需要一步过滤 (example_sce测试数据已经过滤过，故跳过)，代码如下（注意：Marker gene中基因的命名规则与sceset中基因命名规则一致，比如都为ENSEMBL ID或都为Gene Symbol）：

# 对scater对象进行筛选
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函数最大似然法估计（MLE）细胞类型：

print(head(celltypes(fit)))
#> [1] "Group1" "Group2" "Group2" "Group1" "Group1" "Group1"

以及所有MLE参数使用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函数可视化赋值的概率，该函数返回每个细胞和细胞类型的概率矩阵：

pheatmap::pheatmap(cellprobs(fit))

最后，由于这是模拟数据，我们可以检查与真实组值的一致性：

print(table(example_sce$Group, celltypes(fit)))
#>
#>          Group1 Group2
#>   Group1    236      0
#>   Group2      0    264

肿瘤微环境Marker基因示例

对于人肿瘤微环境中的常见细胞类型，cellassign包中提供了一组示例标记。

标记基因可用于以下细胞类型：

• B cells

• T cells

• Cytotoxic T cells

• Monocyte/Macrophage

• Epithelial cells

• Myofibroblasts

• Vascular smooth muscle cells

• Endothelial cells

这些可以通过调用进行访问:

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一起使用，我们可以通过单元格类型矩阵将它们转换为二进制标记：

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

# 对scater对象进行筛选
sce_marker <- sce[intersect(rownames(marker_mat), rownames(sce)),]

局限性

更多单细胞操作

• 收藏 北大生信平台” 单细胞分析、染色质分析” 视频和PPT分享

• Science: 小鼠肾脏单细胞转录组+突变分析揭示肾病潜在的细胞靶标

• Science：通过单细胞转录组测序揭示玉米减数分裂进程 | 很好的单细胞分析案例

• Nature 首次对阿尔茨海默病进行单细胞转录组分析|详细解读

• Cell 深度 一套普遍适用于各类单细胞测序数据集的锚定整合方案

• 骨髓基质在正常和白血病个体中的细胞图谱 Cell,Nature联袂解析

• 癌中之王：基质微环境塑造胰腺癌瘤内结构|Cell

• Nature系列 整合单细胞转录组学和质谱流式确定类风湿性关节炎滑膜组织中的炎症细胞状态 详细解读

• 10X单细胞测序分析软件:Cell ranger，从拆库到定量

• Hemberg-lab单细胞转录组数据分析（一）- 引言

• Hemberg-lab单细胞转录组数据分析（二）- 实验平台

• Hemberg-lab单细胞转录组数据分析（三）- 原始数据质控

• Hemberg-lab单细胞转录组数据分析（四）- 文库拆分和细胞鉴定

• Hemberg-lab单细胞转录组数据分析（五）- STAR, Kallisto定量

• Hemberg-lab单细胞转录组数据分析（六）- 构建表达矩阵，UMI介绍

• Hemberg-lab单细胞转录组数据分析（七）- 导入10X和SmartSeq2数据Tabula Muris

• Hemberg-lab单细胞转录组数据分析（八）- Scater包输入导入和存储

• Hemberg-lab单细胞转录组数据分析（九）- Scater包单细胞过滤

• Hemberg-lab单细胞转录组数据分析（十）- Scater基因评估和过滤

• Hemberg-lab单细胞转录组数据分析（十一）- Scater单细胞表达谱PCA可视化

• Hemberg-lab单细胞转录组数据分析（十二）- Scater单细胞表达谱tSNE可视化

• 如何火眼金睛鉴定那些单细胞转录组中的混杂因素

• 什么？你做的差异基因方法不合适？

• 单细胞分群后，怎么找到Marker基因定义每一类群？

• 在线平台如何做单细胞测序分析全套？有它so easy!

• 植物单细胞转录组的春天来了，还不上车？Science, PC, PP, MP, bioRxiv各一个

• 三人成虎，概率却不足十分之五？

• 一文掌握GSEA，超详细教程

• 这个只需一步就可做富集分析的网站还未发表就被CNS等引用超过350次

• 什么，你算出的P-value看上去像齐天大圣变的庙？

• GO、GSEA富集分析一网打进

• GSEA富集分析 - 界面操作

• 无需写代码的高颜值富集分析神器

• 去东方，最好用的在线GO富集分析工具

• 跨物种单细胞分析发现胰腺导管癌中一类有免疫原性的抗原呈递成纤维细胞

• NCB|心咽发育多样化的单细胞转录轨迹分析

• 七龙珠|召唤一份单细胞数据库汇总

• 用了这么多年的PCA可视化竟然是错的！！！

• 单细胞预测Doublets软件包汇总-过渡态细胞是真的吗？

• Seurat亮点之细胞周期评分和回归

运行环境

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

对单细胞转录组感兴趣的话，不妨留意一下我们将在北京召开的单细胞转录组数据整合分析专题（11月29号-12月01号）：

10X登录纳斯达克，首日大涨35%，市值49亿的技术你还不会？