泊松分布的分布函数_RNAseq的count矩阵数据分布类型

微信公众号：生信小知识
关注可了解更多的生信教程及技巧。问题或建议，请公众号留言

RNA-seq的count矩阵数据分布类型

1. 转录组数据统计推断的难题2. 泊松分布 or 负二项分布？2.1. 为什么泊松分布不行？2.1.1 泊松分布2.1.2 泊松分布的方差与均值关系2.2 为什么负二项分布行？2.2.1 二项分布2.2.2负二项分布3. 方差估计4. 后记

其实很早就在简书上看过有人解释了，这篇文章基本上也是在copy下面这个简书的内容，只不过，现在我是已经完全理解了这篇文章的内容了。

https://www.jianshu.com/p/ad24bb90b972

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

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

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

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

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

我们面临两个核心问题：

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

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

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

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

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

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

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

2.1.1 泊松分布

https://wenku.baidu.com/view/7bf1eeb62bf90242a8956bec0975f46526d3a7ca.html?fr=search-1-aladdin

https://zhuanlan.zhihu.com/p/24711669

泊松分布的公式如下：

满足泊松分布的3个要求：

• 小概率时间

• 相互独立

• 概率稳定

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

λ是所有值的均值。

经典例子如下：

已知某家小杂货店，平均每周售出2个水果罐头。请问该店水果罐头的最佳库存量是多少？

假定不存在季节因素，可以近似认为，这个问题满足以下三个条件：

(1)顾客购买水果罐头是小概率事件。

(2)购买水果罐头的顾客是独立的，不会互相影响。

(3)顾客购买水果罐头的概率是稳定的。

在统计学上，只要某类事件满足上面三个条件，它就服从"泊松分布"。

在这个例子中，泊松分布中各个参数的含义如下：

• P——每周销售k个罐头的概率。

• X——水果罐头的销售变量。

• λ——每周水果罐头的平均销售量，是一个常数，本例为2

根据公式，计算出对应的p值：

每周买出的罐头数目k	概率p	累计概率
0	0.135	0.135
1	0.271	0.406
2	0.271	0.677
3	0.180	0.857
4	0.090	0.947
5	0.036	0.983
6	0.017	1.000

从上表可见，如果存货4个罐头，94.7%的概率不会缺货(平均每19周发生一次)；

如果存货5个罐头，98%的概率不会缺货(平均59周发生一次)。

2.1.2 泊松分布的方差与均值关系

对于泊松分布而言，其均值和方差是相等的，但是我们的数据确不符合这样的规律。通过计算所有基因的均值和方差，可以绘制如下的图片：

图中的点纵坐标为基因在所有样本中的均值，横坐标为基因在所有样本中的方差。

图中画了一条斜率为1的直线，代表了泊松分布的均值和方差的分布。

可以看到，真实数据的分布是偏离了泊松分布的，方差明显比均值要大。

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

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

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

2.2 为什么负二项分布行？

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

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

2.2.1 二项分布

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

它的概率分布图如下：

代码如下：

# 二项分布
if (T) {
  x 100000, size=100, prob=0.9)
  hist(x,
       xlim = c(min(x), max(x)),
       probability = T,  #概率函数，默认为频率函数
       nclass =  max(x) - min(x) +1,
       col =  'lightblue',
       main = "Binomial Distribution, n=100, p=.9" )
  lines(density(x,bw=1,kernel = 'gaussian'),col = 'blue',lwd = 3)
}

实例如下：

生产产品，结果只有正品和次品(二项)，求5次中2次是正面的概率。

2.2.2负二项分布

负二项分布描述的也是伯努利实验，不过它的目标事件变成了：对于Bernoulli过程，我们设定，当某个结果出现固定次数的时候，整个过程的数量。

比如我们生产某个零件，假设每个零件的合格与否都是相互独立的，且分布相同，那么当我们生产出了五个不合格零件时，一共生产了多少合格的零件，这个数量就是一个负二项分布，公式如下：

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

它的概率分布如下：

代码如下：

# 负二项分布
if (T) {
  x 100000, size=100, prob=0.9)
  hist(x,
       xlim = c(min(x), max(x)),
       probability = T,
       nclass =  max(x) - min(x) +1,
       col =  'lightblue',
       main = "Negative Binomial Distribution, n=100, p=.9" )
  lines(density(x,bw=1,kernel = 'gaussian'),col = 'blue',lwd = 3)
}

经过一堆数学神奇世界的计算后，我们可以得到：

其中u为mean，σ为sd，σ²为variance。

他们之间的关系经过一堆数学神奇世界的计算后，我们可以得到：

其中的符号这里再次注明下：

• 合格率为p

• 抽到r个次品时停止抽样

• 正品正好为k个

从上面的式子可以看出，均值是方差(σ^2)的二次函数，方差随着均值的增加而进行二次函数形式的递增，正好符合上文均值与方差分布图的情况

其中r被称为dispersion parameter。

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

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

当 r -> 0 时，此时overdispersion

3. 方差估计

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

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

4. 后记

其实关于到最后你会慢慢发现，很多生信分析的背后，其实是在玩统计，你需要了解自己手上数据的分布符合什么模型，根据自己数据的类型去选择合适的工具罢了。当然，我也不能直接这么草率的就下这个结论，毕竟我现在也只是个菜鸡，上面那句话也只是我现在的体会 罢了，大家看着玩玩就好啦~

有兴趣的可以沟通交流下哟~