单细胞RNA-seq数据分布的选择

最新推荐文章于 2024-06-09 10:27:49 发布

goodljq

最新推荐文章于 2024-06-09 10:27:49 发布

阅读量1.7k

点赞数 6

分类专栏：知识总结文章标签：统计模型生物信息学

本文链接：https://blog.csdn.net/goodljq/article/details/117114849

版权

知识总结专栏收录该内容

6 篇文章 6 订阅

订阅专栏

在单细胞RNA测序数据（single-cell RNA sequencing，以下简称scRNA-seq）的分析中，对数据分布的模型假设是很重要的一步，合理的数据分布模型直接影响了后续的分析流程或者是方法开发。例如单细胞差异表达分析方法就是基于数据的分布，通过对分布参数的假设检验，进而判断两组数据是否存在差异。那么人们是如何一步步探索适合的scRNA-seq数据的分布的呢？现阶段的发展又是怎样的？有哪些新产生的问题？本文将就这这些问题加以梳理。这里的分布指的是某一基因转录出RNA的量在一组细胞中（gene expression across cells）的分布。

正态分布

说起分布，人们的第一反应大概都是正态分布。因为简单分布中正态分布最能符合现实生活中很多变量的观测。另外根据中心极限定理，如果一个特定事件受多个因素的影响，而每个因素对结果的影响都很小的时候，各种因素作用的和服从正态分布。如下图所示：

图片来源：https://www.zhihu.com/question/26854682

但是细胞中RNA数量的值是离散的，而正态分布是连续分布。另外，scRNA-seq数据往往不是对称的，这与正态分布也不相吻合。因此，正态分布不适合用作scRNA-seq数据分布。

不过不同细胞转录出的RNA的量往往存在数量级上的差异，因此在数据分析时常常将数据做对数变换： $x^{'} = log_{10}(x+1)$ ，而人们发现变换后的结果近似服从正态分布，因此scRNA-seq数据也被认为是服从对数正态（log-normal）分布的。

泊松分布

人们也尝试从测序的机理上来建模scRNA-seq数据的分布。如下图所示，一个细胞中部分基因各自转录出若干RNA，假设每条RNA被测序工具捕捉到的概率为 p ，那么这一事件服从Bernoulli分布；而从总数为n的所有RNA中捕捉到Gene1对应的RNA数量 $n_{g1}$ 就服从二项分布：

$n_{g1}\sim Binomial(n,p) \\$

而一个细胞中转录出的RNA数量 n 非常多，捕捉到某一条特定RNA的概率 p 也相当小，因此二项分布就近似成为了泊松分布：

$n_{g1}\sim Poisson(\lambda =np) \\$

图片来自组会ppt，浩翔师兄手绘，在此鸣谢~

但泊松分布中均值始终等于方差，这一点不符合scRNA-seq数据的实际情况。下图展示了scRNA-seq数据的实际分布情况（散点，每个点代表一个基因对应RNA数量的统计量）和泊松分布的理论分布（直线）。可以看到，二者差别非常大，随着RNA均值的增加，其方差和均值之间的差距越来越大，这一现象称为"Over-dispersion"。因此我们仍然需要寻找更合理的分布来拟合scRNA-seq数据。

图片来源：https://www.biostars.org/p/84445/

负二项分布

进一步思考，在泊松分布中均值等于方差，而唯一参数 $\lambda$ 的值是不变的，如果它是变化的呢？经过一番探索后，人们发现如果 $\lambda$ 的先验分布取伽马分布的时候，即 $\lambda \sim Gamma(\alpha,\beta)$ ，后验分布满足负二项分布，因此负二项分布也称为Gamma-Possion分布。负二项分布包含两个参数： $NB(r,p)$ ，其均值为：

$\mu=\frac{pr}{1-p} \\$

方差为：

$\sigma^2=\frac{pr}{(1-p)^2}=\mu+\frac{\mu^2}{r} \\$

均值不等于方差，因此可以解决scRNA-seq数据中over-dispersion的问题。同时负二项分布可以更好地拟合多种数据分布，如下图所示：

负二项分布随参数的变化，图片来自wikipedia

负二项分布算是scRNA-seq数据分析的模型中广泛应用的数据分布了，由此开发出的差异表达与缺失值填补（单细胞领域习惯称为imputation）方法实用性也更强，比如edgeR[1], DESeq[2] ,SAVER[3] 等等。

零膨胀负二项分布

在广泛应用负二项分布的同时，人们也发现scRNA-seq数据还有一个特点，那就是零值非常多，下图给出了一个真实scRNA-seq数据中零表达基因比例的直方图 [4]：

scRNA-seq数据中完全不表达的基因占比也是很高的

由于基因表达数据中的零值既可能来自生物过程中不表达的基因（称为True Zero），还可能来自测序过程中由于技术原因导致的丢失（称为False Zero 或者Dropout Zero），因此人们尝试在NB模型中加入一个零膨胀因子，用零膨胀负二项分布（Zero-Inflated Negative Binomial）来建模scRNA-seq数据：

$f_{ZINB}(m|\theta,r,p)=\theta\cdot I_0(m)+(1-\theta)\cdot f_{NB}(m|r,p) \\$

其中 $\theta$ 表征零值的比例， $r,p$ 为负二项分布的参数， $I_0(\cdot)$ 为示性函数，当自变量为0时值为1，否则为0。

对于上面提到的真零假零的问题，ZINB模型也给予了回答，假设测序过程中RNA的捕获率为 $\beta$ ，那么不同 $\beta$ 对应的分布如下图所示 [4]：

图片来源：Miao et al., Bioinformatics, 2018.

红色bar表示真实不表达基因产生的0值，蓝色bar为负二项分布，当捕获率不断降低时，NB分布会左移向0值靠拢，产生“假零”。而分析ZINB中两个部分的比例可以计算一个基因零值是假零的概率（Dropout Rate）。

$d=\frac{(1-\theta)\cdot f_{NB}(0|r,p)}{\theta+(1-\theta)\cdot f_{NB}(0|r,p)} \\$

混合分布

当然ZINB本身就是一种混合分布。只是NB和ZINB在scRNA-seq数据中应用很多，为了保证逻辑性单独介绍了一下。人们还提出了很多混合分布模型来拟合scRNA-seq数据，下面简单介绍三种方法：

Gamma-Normal 分布

混合分布的形式如下：

$f(x)=\lambda\cdot Gamma(x;\alpha,\beta)+(1-\lambda)\cdot Normal(x;\mu, \sigma) \\$

其中Gamma分布用来表示Dropout现象，正态分布用来表示真实基因分布。当然需要注意这里需要对原始的scRNA-seq数据进行对数变换（ $log_{10}$ transform）。Gamma-Normal分布的一个应用是单细胞imputation方法scImpute [5]。

Beta-Poisson分布

类似上面提到的Gamma-Poisson分布，一种观点 [6] 认为利用Beta分布的两个参数来表示转录动力学中的Burst Frequency和Burst Size这两个参数，可以更好地模拟转录过程。因此假设Poisson分布中的参数 $\lambda$ 满足Beta的先验分布，scRNA-seq数据的分布可以表示为：

$BP(x|\alpha,\beta,\lambda)=Poisson(x|\lambda\sim Beta(\alpha,\beta)) \\$

Logistic回归 + 正态分布

在这种混合分布模型中，先用logistic回归判断一个基因是否表达：

$logit(Pr(Z_{ig}=1))=X_i\beta^D_g \\$

之后再用正态分布拟合表达的基因

$Pr(Y_{ig}=y|Z_{ig}=1)=N(X_i\beta^C_g, \sigma^2_g) \\$

这种分布的一个应用是单细胞差异表达基因检测方法MAST [7]。

补充说明

ZINB分布在单细胞RNA数据分析领域越来越受到研究者的青睐，但也出现了一些声音反对这种假设，反对者认为仅使用NB而不增加零膨胀项仍然可以反映数据分布。对于时下常用的两套测序方法——Plane-based sequencing & Droplet-based sequencing来说，前者计数对象是reads而后者计数对象是UMI (Unique molecular identifiers)，而有研究表明UMI计数的分布并不是ZINB的，如下图 [8] ：

图片来源：Chen et al., BMC Genome Biology, 2018.

图中x、y两个axis指的是两个表达量（reads or UMI）相近的细胞。SMART-seq测量的是reads，CEL-seq测量的是UMI，中间图是CEL-seq测量结果decomplex成reads的计数结果（这里涉及测序方法的讨论，感兴趣的朋友可以多查阅测序方法的文献，这里基于篇幅就不深入讨论了）。可以看到reads的分布是双峰的，UMI的分布是单峰的，作者用NB很好的拟合了UMI的分布结果，因此认为不必要增加零膨胀项来拟合分布。

个人浅见是，对于scRNA-seq数据，可以用ZINB来统一建模，因为ZINB可以同时拟合双峰和单峰的数据分布，对于上图右侧的UMI分布，也可以通过ZINB来拟合。当然实际应用中的选择还有待探索，也欢迎大家来一起讨论为什么同样是测序，测reads得到的结果分布和测UMI的分布差别很大，是方法机理上的不同导致的，还是测序技术本身的不同带来的偏差，或是其他。

另外，本文中所有的分布都只是“假设”，现在领域并没有金标准（gold standard）。大家在使用某种特定模型时，关注其背后对数据分布的假设是必要的，对于数据预处理（比如是否做对数变换）也有指导意义。