RNAseq与负二项分布

RNA-seq中的那些统计学问题（一）为什么是负二项分布？

WHY SEQUENCING DATA IS MODELED AS NEGATIVE BINOMIAL

RNA-seq分析中的dispersion，你知道吗？

DESeq2的建模原理及简单用法

1. 转录组数据统计推断的难题

在RNA-seq中进行两组间的差异分析是最正常不过的了。

我们在其它实验中同样会遇到类似的分析，通常，我们可以用方差分析判定两组“分布”数据间是否存在显著差异。原理是：当组间方差大于组内方差（误差效应），并且统计学显著时，则认为组间处理是可以引起差异的。

有伙伴肯定要问，转录组数据到底有什么了不起的？它们为什么不能用我们熟悉的算法简单地进行计算？

其实统计学家也很无奈啊，看看我们转录组实验得到的这些数据吧：我们的实验只进行少得可怜的生物学重复（n<10），而且，任何基因的表达量都不能是负数，这些数据并不符合正态分布，用于表征表达量的counts是非连续的（芯片信号是连续的），RNA-seq数据的离散通常是高度扭曲的，方差往往会大于均值……，就这些奇怪的特征，使得准确估计方差并没有想象的那么容易。

我们面临两个核心问题：

• 基因表达数据适合用什么统计学分布进行差异显著性检验？

• 如何利用少量生物学重复数据估算基因表达的标准差？

2. 泊松分布 or 负二项分布？

从统计学的角度出发，进行差异分析肯定会需要假设检验，通常对于分布已知的数据，运用参数检验结果的假阳性率会更低。转录组数据中，raw count值符合什么样的分布呢？

count值本质是reads的数目，是一个非零整数，而且是离散的，其分布肯定也是离散型分布。对于转录组数据，学术界常用的分布包括泊松分布 (poisson)和负二项分布 (negative binomial)两种。

2.1. 为什么泊松分布不行？

首先有必要简单地介绍一下泊松分布

泊松分布适合于描述单位时间（或空间）内随机事件发生的次数（事件发生的次数只能是离散的整数）。如某一服务设施在一定时间内到达的人数，电话交换机接到呼叫的次数，汽车站台的候客人数，机器出现的故障数，自然灾害发生的次数，一块产品上的缺陷数，显微镜下单位分区内的细菌分布数等等。

$$P(X=k)=\frac{\lambda^k }{k!}e^{-\lambda},\quad k=0,1,…$$

泊松分布大概长这样：


λ是波松分布所依赖的唯一参数。 λ值愈小分布愈偏倚，随着λ的增大，分布趋于对称。 当λ=20时分布接近于正态分布；当λ=50时， 可以认为波松分布呈正态分布。

在数据分析的早期，确实有学者采用泊松分布进行差异分析，但是发展到现在，几乎全部都是基于负二项分布了，究竟是什么因素导致了这种现象呢？为了解释这个问题，我们必须提到一个概念 overdispersion。

dispersion指的是离散程度，研究一个数据分布的离散程度，我们常用方差这个指标。对于泊松分布而言，其均值和方差是相等的，但是我们的数据确不符合这样的规律。通过计算所有基因的均值和方差，可以绘制如下的图片：

mean <- log10(apply(x, 1, mean))
var <- log10(apply(x, 1, var))
plot(x = mean, y = var, pch = 20)
abline(a = 0, b=1)

library(DESeq2)
library(matrixStats)
library(magrittr)
library(ggplot2)

dds <- readRDS("ddsMat.rds")

mvDat <- data.frame(
  Mean = log10(unname(rowMeans(counts(dds,normalized = TRUE)))),
  Variance = log10(rowVars(counts(dds,normalized = TRUE)))
)

ssDat <- smooth.spline(x = mvDat$Mean, y = mvDat$Variance, spar = 1) %$%
  data.frame(x = x, y = y)

ggplot(mvDat, aes(x = Mean, y = Variance)) +
  geom_point(size = 0.5, alpha = 0.6, colour = "darkblue") +
  geom_line(data = ssDat, aes(x = x, y = y), colour = "red", size = 0.5) +
  theme_bw() +
  xlim(0,10) +
  ylim(0,10) +
  theme(plot.title = element_text(hjust = 0.5))

横坐标为基因在所有样本中的均值，纵坐标为基因在所有样本中的方差，直线的斜率为1，代表泊松分布的均值和方差的分布。可以看到，真实数据的分布是偏离了泊松分布的，方差明显比均值要大。

如果假定总体分布为泊松分布， 根据我们的定量数据是无法估计出一个合理的参数，能够符合上图中所示分布的，这样的现象就称之为overdispersion。

由于真实数据与泊松分布之间的overdispersion，选择泊松分布分布作为总体的分布是不合理。

以上只证明了泊松分布是个不太恰当的分布估计，那怎么证明负二项分布就是合适的分布估计呢？

2.2. 为什么负二项分布行？

主要是从均值与方差之间的关系去证明

同样的，也先简单介绍一下负二项分布：

伯努利试验（Bernoulli experiment）是在同样的条件下重复地、相互独立地进行的一种随机试验，其特点是该随机试验只有两种可能结果：发生或者不发生。 我们假设该项试验独立重复地进行了n次，那么就称这一系列重复独立的随机试验为n重伯努利试验，或称为伯努利概型。

二项分布描述的是n重伯努利实验，在n重贝努利试验中，事件A恰好发生x(0≤x≤n)次的概率为：

$$P_n(x)=C_n^x p^x(1-p)^{n-x}$$

它的概率分布图如下：


负二项分布描述的也是伯努利实验，不过它的目标事件变成了：对于Bernoulli过程，我们设定，当某个结果出现固定次数的时候，整个过程的数量，比如我们生产某个零件，假设每个零件的合格与否都是相互独立的，且分布相同，那么当我们生产出了五个不合格零件时，一共生产了多少合格的零件，这个数量就是一个负二项分布，公式如下：

$$f(k;r;p)=P(x=k)=C_{r+k-1}^k p^k(1-p)^r$$

该公式描述的是，在合格率为p的一堆产品中，进行连续有放回的抽样，当抽到r个次品时，停止抽样，此时抽到的正品正好为k个的概率

它的概率分布如下：


p=0.5, r=5

负二项分布的均值和方差分别为( 如何推导出来的？)：

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

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

将p用μ表示，得到：

$$p=\frac{\mu}{\mu+r},\quad 1-p=\frac{r}{\mu+r}$$

将上一步推出的p和1-p带入到方差的表达式中，得到：

$$\sigma^2=\frac{\mu^2}{r}+\mu$$

记1/r=α，则

$$\sigma^2=\mu+\alpha\mu^2$$

从上面的式子可以看出，均值是方差的二次函数，方差随着均值的增加而进行二次函数形式的递增，正好符合上文 2.1. 为什么泊松分布不行？ 部分均值与方差分布图的情况

其中α和r被称为dispersion parameter

负二项分布与泊松分布的关系，可以用α或r推出：

当 r -> ∞ 时，α -> 0，此时 σ2= μ，为泊松分布；

当 r -> 0 时，α -> ∞，此时overdispersion

3. 方差估计

在生物学重复很少时，我们是很难准确计算每个基因表达的标准差的（相当于这个数据集的离散程度）。我们很可能会低估数据的离散程度。

被逼无奈的科学家提出了一个假设：表达丰度相似的基因，在总体上标准差应该也是相似的。我们把不同生物学重复中表达丰度相同的基因的总标准差取个平均值，低于这个值的都用这个值，高于这个值的就用算出来的值。