dnbinom(x, size, prob)：返回发生x次失败事件的概率
rnbinom(n, size, prob)：返回每组发生失败事件的次数
pnbinom(q, size, prob)：返回累积概率
qnbinom(p, size, prob)：返回相应分位点x

参数解释：

size：成功实验次数
mu：有size次成功平均实验次数
P：成功率为p （mu与p只需要一个参数）
目的：成功实验次数为r，成功率为p，平均每次需要进行k次失败实验的分布，总实验次数为n=k+r

理论分布：

k <- seq(0,10,by = 1)
#每次有1次成功实验且成功率为1/2，所需要的实验次数为k+1次，k为横坐标。k表示在r次成功（r=1）之前有k次失败，k=0表示第一次（0+1=1）实验就成功，概率为50%；k=1表示第二次实验（1+1=2）才成功，概率为50%*50%=25%；以此类推。 
y <- dnbinom(k, size = 1, prob = 1/2)
plot(k,y)

image.png

模拟:

成功实验次数为1时，平均实验次数为4的分布。生成100000个数据
k1 <- rnbinom(100000, mu = 4, size = 1) 
h1 <- hist(k1, breaks = 20)
mean(k1) #4.00653
成功实验次数为1时，成功概率为0.4的分布。生成100000个数据
k2 <- rnbinom(100000, p = 0.4, size = 1) 
h2 <- hist(k2,breaks=20)
mean(k2)#1.50545

image.png
image.png

当成功率p=0.5,成功次数size=5，失败次数为k时,进行n组实验，方差与均值的关系

n <- seq(1,10000)
mu <- c()
var <- c()
for (x in n){
  k3<- rnbinom(x, p = 0.5, size = 5)
  mu <- c(mu,mean(k3))
  var <- c(var,var(k3))
}
data <- data.frame(n,diff=var-mu)
library(ggplot2)
ggplot(data,aes(x=n,y=diff))+geom_point()

可以发现当n逐渐增大时，均值与方差的差值逐渐稳定，接近于一个常数

image.png

验证当p,size不变，定义符合p为0.5，size为5,均值为5的1000个数据，随后均值变大到1，2，3....，方差有什么变化

fold <- seq(1,500)
mu3 <- c()
var3 <- c()
for (f in fold){
  k3 <- rnbinom(10000, mu = f, size = 5) 
  mu3 = c(mu3, f)
  var3 <- c(var3, var(k3))
}
mean.var33= data.frame(mu3,var3)
ggplot(mean.var33,aes(x=mu3,y=var3))+geom_point()+geom_abline(intercept = 0, slope = 1)

发现均值和方差近似一元二次方程，至于是不是得后面验证

image.png

将均值和方差对数化，发现均值和方差近似一元二次方程

fold <- seq(1,500)
mu3 <- c()
var3 <- c()
for (f in fold){
  k3 <- rnbinom(10000, mu = f, size = 5) 
  mu3 = c(mu3, f)
  var3 <- c(var3, var(k3))
}
mean.var33= data.frame(mu3,var3)
log10.mean3.var3 = log10(mean3.var3)
ggplot(log10.mean3.var3,aes(x=mu3,y=var3))+geom_point()+geom_abline(intercept = 0, slope = 1)

image.png

检验一组数据是否符合负二项分布

计算所有基因的均值和方差，发现不相等，存在过离散现象
mean <- log10(apply(data, 1, mean))
var <- log10(apply(data, 1, var))
plot(x = mean, y = var, pch = 20)
abline(a = 0, b=1)
检测这组数据是否确实存在过离散现象，显著性检验的p值小于0.05，表明确实存在过度离势，符合负二项分布
install.packages("qcc")
library(qcc)
qcc.overdispersion.test(data,type="poisson")
Overdispersion test Obs.Var/Theor.Var Statistic p-value
       poisson data          2.493859  249383.4       0
可以假设每一个reads计数是一个事件，计到该基因则为成功，计到其他基因则为失败。

dispersion指的是离散程度，研究一个数据分布的离散程度，我们常用方差这个指标。对于泊松分布而言，其均值和方差是相等的，但是我们的数据确不符合这样的规律。通过计算所有基因的均值和方差，可以绘制如下的图片,横坐标为基因在所有样本中的均值，纵坐标为基因在所有样本中的方差，直线的斜率为1，代表泊松分布的均值和方差的分布。可以看到，真实数据的分布是偏离了泊松分布的，方差明显比均值要大。

image.png