负二项分布

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
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 216,240评论 6 498
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 92,328评论 3 392
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 162,182评论 0 353
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,121评论 1 292
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,135评论 6 388
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,093评论 1 295
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,013评论 3 417
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,854评论 0 273
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,295评论 1 310
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,513评论 2 332
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,678评论 1 348
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,398评论 5 343
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,989评论 3 325
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,636评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,801评论 1 268
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,657评论 2 368
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,558评论 2 352

推荐阅读更多精彩内容