Learn-SciBet

1. 背景介绍

细胞类型注释是scRNA-seq里非常基础的一步，常见的策略是将基于无监督分群结果的基因差异表达分析结果作为marker,结合先验知识，判断该分群为何种类型。

近年来，随着scRNA-seq的数据集的逐渐积累，也有很多研究组提出了一些基于已有分群结果进行有监督的细胞类型注释的方法，比如scmap，或者Seurat里的TransferData。此外，也有评估相关方法的benchmark研究：

SciBet是张泽民老师组最近新发表的一个快捷（实测真的比Seurat的TransferData快，而且效果确实差不多），方便的进行有监督的细胞类型注释的方法。论文见SciBet as a portable and fast single cell type identifier，相关报道见Nature Communications | 张泽民课题组发表单细胞转录组数据快速注释新方法

2. 方法原理简介

如果读者只对软件使用感兴趣的本节可以略过。

根据上面提到的那篇报道，如果用一句话概括SciBet的原理，应该是这样的：

张泽民实验室的博士生李辰威、刘宝琳联合任仙文副研究员开发的SciBet则有效地解决了这一问题：他们从“同一类型的单细胞表达谱服从同一多项分布”这一基本假设出发，对训练集数据中不同细胞类型分别进行建模，进而通过极大似然估计来对测试集细胞进行有监督注释。

接下来我们根据作者论文的方法部分和补充材料，学习一下这个巧妙地思路。分为如下几个部分：

• 预备知识
• 单细胞表达谱的统计模型
• 根据表达谱计算信息熵
• E-test
• 构建有监督的细胞类型分类模型

我们假定读者已经修过概率论等相关课程。

2.1 预备知识

预备知识1：

我们需要向读者回顾关于负二项分布的知识，如果一个离散随机变量Y的pmf(probability mass function)为:

那么我们称其服从负二项分布。记为Y ∼ NB(r, p)。

预备知识2：

对于两个随机变量X和Y，在给定X = x下，Y的条件pmf为：

预备知识3：

关于负二项分布，有如下结论：

对于随机变量Y, Λ，若Y|Λ ∼ Poisson(Λ), Λ ∼ Gamma(α, β), 则 其中

证明：

因为：

所以

根据预备知识3，我们有：

预备知识4：

负二项分布是泊松分布和伽马分布的混合分布。

预备知识5(该结论为概率论和数理统计的常识)：

对于随机变量的pdf为f_X(x)，做变换Y = g(X),其中g为单调函数，X和Y各自样本空间的𝒳和𝒴分别满足：
𝒳 = {x|f_X(x) > 0}, 𝒴 = {y|y = g(x), x ∈ 𝒳}
若f_X(x)在𝒳中连续且逆变换g^− 1(x)在𝒴中连续可微，则新的随机变量Y的pdf为：

该结论为Statistical Inference(Casella Berger著)中的定理2.1.5，证明见该书。

预备知识6（只有这个结论作者在补充材料里给了证明）：

The scaling property of the Gamma distribution

若 X ∼ Gamma(α, β), 则

证明: 由题设：

根据题设，由预备知识5,有：

故

预备知识6：

若独立同分布样本X₁, …, X_n服从X_i ∼ Poisson(λ), 且观测值分别为x₁, …, x_n，则参数λ的矩估计量和极大似然估计量相等，均为

证明很简单，在很多教材也能找到，故从略。

预备知识7：

若独立同分布样本X₁, …, X_n服从X_i ∼ Gamma(α, β), 且观测值分别为x₁, …, x_n，则参数α, β的极大似然估计满足：

证明：

由题设：

故取对数后极大似然函数为：

由

解得

预备知识8：

若随机变量X的pdf为f_X(x),样本空间为𝒳, 定义information generating function:
T(u) = ∫_𝒳f_X^u(x)dx
则：
H(X) =  − T^′(1)

证明：(严格证明的话，得要考虑函数f_X(x)和𝒳本身的性质，我们这里假定它们可以满足积分号下求导，如果要严格的话，请参考测度论相关教材)

而随机变量X的信息熵的定义为：

H(X) =  − ∫_𝒳f_X(x)ln f_X(x)dx

故H(X) =  − T^′(1)

预备知识9：

若X ∼ Gamma(α, β)，则其信息熵S为

由预备知识8经过计算即可证明，从略。

预备知识10：

多项分布

注：更多关于伽马分布的知识可见：理解Gamma分布、Beta分布与Dirichlet分布

2.2 单细胞表达谱的统计模型

经过大量的统计分析和后续的实验验证（相关证据可参考这篇文献）有这样一个经验性的结论：

观察到单细胞基因表达的count(比如UMI count)的分布可以用负二项分布很好的拟合,且相同细胞类型的单细胞表达谱服从同一个分布。

结合2.1中的预备知识，我们可以将单细胞基因表达的count表示为泊松分布的伽马分布的混合分布。所以作者参考了SAVER可以进行如下建模：

假设我们现在有C个细胞，m个基因的原始表达谱数据，里面的数值为reads count或者是UMI count。如果我们记观察得到细胞c的某个基因i的UMI count为Z_ic, i = 1, …, m, c = 1, …, C，那么对于同一类型的细胞而言，有：

Z_ic ∼ Poisson(λ_is_c), λ_i ∼ Gamma(α_i, β_i)

其中λ_i表示基因i在细胞c中的真实表达量，s_c = ∑_cZ_ic表示这个细胞中的UMI总数，与测序深度有关。而α_i, β_i则是两个参数表征某个细胞类型中的基因的真实表达分布的参数。

接下来的可以根据预备知识里的结论进行参数估计：

由预备知识6，我们很容易得到这也为我们常用的进行normalized的策略提供了一种依据。而由预备知识7，有。

2.3 根据表达数据计算信息熵

不是所有的基因都是对后续的统计学习有用的，需要进行特征选择，也就是说，挑选出那些能表示不同群细胞之间表达差异的基因。本文的新意是基于信息熵(也就是香农熵)的概念引入了新的进行特征选择的方法：E-test。在讲E-test之前，我们需要看看作者是如何实现从表达量中计算信息熵的。

根据2.2我们知道可以将观测得到单细胞表达gene countZ_ic表示成Z_ic与真实表达量λ_ic的混合分布，真正能反映不同细胞之间表达差异的是λ_ic的分布。所以接下来要计算λ_ic的分布的信息熵。

对于相同类型的细胞而言，根据2.1中的结论9，真实表达量λ_i的信息熵为：

用代入（1）可得：

记：

并且记C个细胞的平均normalized表达量为X_i, 显然有

根据上述记号，（1）最终化简为：

接下来，我们考虑道不同的细胞类型。若细胞c, c = 1, …, C_j属于细胞类型j，定义所有属于j的细胞的平均标准化的表达量为，根据上面的结论，可得：

其中

X_ij是直接可以通过实验数据计算的，而h_ij则需要估计。作者假设h_ij是只是基因特异的参数，也就是h_ij = h_i。理由如下：

若λ_{i, c ∈ j}为随机变量λ_ij ∼ Γ(α_ij, β_ij)的观测值，如果我们记基因i从细胞类型j到细胞类型j^′的表达量的fold change为F_{i, j → j^′}, 并且假设λ_ij^′ = F_{i, j → j^′}λ_ij，那么由2.1中的预备知识6，可知 。 故α_ij = α_i, 结合（6），h_ij = h_i

最终，我们可以得到

2.4 E-test

首先，零假设为所有不同类型细胞都是从同一个细胞类型（记为 group 0）中均匀随机采样，那么基因的平均表达量为, 则group 0 的信息熵可以计算为：

接着计算基因i在所有细胞类型j中与group 0 的信息熵的差之和：

(8)-(7)并求和，得：

其中

利用Jesen不等式可以证明GM_i ≤ AM_i，故ΔS ≥ 0

若要进行假设检验，还需要计算ΔS的显著性，作者的策略是基于置换检验的：

若所有预先定义分群的细胞类型的细胞均来自同一个样本，则对任意的标签为细胞类型j的细胞的size-factor normalized的表达量，根据中心极限定理，单细胞数目足够多的时候，X_ij ∼ N(μ_i, σ_i)，很容易得到参数的无偏估计，所以置换被简化成了每次从分布n个不同的随机数X_i。接下来就是这么生成一个ΔS的分布，然后计算这个分布中大于从真实数据中测得的ΔS比例，作为显著性。

默认的Feature为500个基因。

2.5 构建有监督学习的模型

根据2.4，在训练集中，我们可以进行特征选择选出一些“informative gene”进行模型训练。

作者假设从相同的转录出的mRNA是不可区分的，而且每个mRNA的产生是相互独立的，记录对于细胞类型为j的细胞，基因i产生一个mRNA的概率为p_ij,若有m个informative gene 则对细胞类型j,我们得到了一个随机向量p_j = (p_1j, ⋯, p_mj)其中∑_ip_ij = 1且其服从多项式分布,则 对于属于细胞类型j的细胞y，其后验表达谱y = (y₁, …, y_n)可以计算为：

其中概率p_ij可估计为

同样的，在测试集中，未知细胞类型的细胞y属于细胞类型j的概率也为

其中p_ij是训练集中学习到的参数。

如果，我们引入,由极大似然的原则可知，细胞y最有可能的细胞类型ĵ为

2.6 方法总结

可以用原文献中的流图对SciBet进行总结：

3. 软件使用

3.1 在线版

在线版使用请按照官网的Online classification教程。

值得一提的是，作者提供了很从不同组织，不同实验条件的单细胞测序数据中训练好的Signature：

这些不同的signature可以供不同研究者结合自己的兴趣使用。

3.2 本地版

3.2.1 安装

if (!requireNamespace("devtools", quietly = TRUE)) install.packages("devtools")
devtools::install_github("PaulingLiu/scibet")

如果出现错误，请看相关issue

3.2.2 作者提供的测试数据

按照官网的教程E-test and SciBet；下载所需的数据；然后后按照其说明文档进行即可。

接下来，我们看如何用使用SciBet结合Seurat进行单细胞分析。

3.2.3 SciBet结合Seurat

为了说明问题，我们选择Seurat自带的pbmcsca数据集, 该数据集已经提供了预先定义好的细胞标签，代码如下：

library(Seurat)
library(pbmcsca.SeuratData)
library(ggplot2)
library(scibet)
library(tidyverse)
library(viridis)

####0.---define utilized function--------
####---Export expr data from 10x to tibble----
#' @param seuv3 a seuv3 object 
#' 
myGetExpr <- function(seuv3,...){
	expr <- GetAssayData(object = seuv3, slot = "data")
	expr <- as_tibble(t(as.matrix(expr)),rownames = NA)
	return(expr)
}

###1.load data and preprocessing
data("pbmcsca")
###---avoid warning-----
pbmcsca <- UpdateSeuratObject(pbmcsca)
###---show predifined cell type------
table(pbmcsca$CellType)

# B cell              CD14+ monocyte 
# 5020                        5550 
# CD16+ monocyte                 CD4+ T cell 
# 804                        7391 
# Cytotoxic T cell              Dendritic cell 
# 9071                         433 
# Megakaryocyte         Natural killer cell 
# 977                        1565 
# Plasmacytoid dendritic cell                  Unassigned 
# 164                          46 

###---split data-----
pbmc.list <- SplitObject(pbmcsca, split.by = "Method")
###---normalize-------
for (i in names(pbmc.list)) {
	pbmc.list[[i]] <- NormalizeData(pbmc.list[[i]], verbose = FALSE)
}
names(pbmc.list)

[1] "Smart-seq2"          "CEL-Seq2"            "10x Chromium (v2) A"
[4] "10x Chromium (v2) B" "10x Chromium (v3)"   "Drop-seq"           
[7] "Seq-Well"            "inDrops"             "10x Chromium (v2)"  

测试 reference base模式的效果

###----2. test scibet----
###----define reference and query-----
reference <- pbmc.list[[1]]
query <- pbmc.list[[2]]


###----test query reference mode----
reference.expr <- myGetExpr(reference)
query.expr <- myGetExpr(query)


reference.label <- as.character(reference$CellType)
test.label <- as.character(query$CellType)

reference.expr <- cbind(reference.expr,label=reference.label)

prd.label <- SciBet(train = reference.expr, test = query.expr)
Confusion_heatmap(test.label, prd.label)
ggsave(filename = "res/fig/learn_scibet_confusionheatmap_refmode.pdf",
	   width = 6,height = 6)

准确率为

num1 <- length(test.label)
num2 <- tibble(
	ori = test.label,
	prd = prd.label
) %>%
	dplyr::filter(ori == prd) %>%
	nrow(.)

num2/num1

0.851711

测试用作者提供的训练好的模型

###----test load_model mode-----

###using 30 major cell types signature----

model <- readr::read_csv(file = "http://scibet.cancer-pku.cn/major_human_cell_types.csv")
model <- pro.core(model)

prd <- LoadModel(model)
prd.label <- prd(query.expr)

Confusion_heatmap(test.label,prd.label)
ggsave(filename = "res/fig/learn_scibet_confusionheatmap_signaturemode.pdf",
	   width = 6,height = 6)