次世代シークエンシングデータと二項分布・ポアソン分布・負の二項分布

次世代シークエンシングでは、たくさんの座位が何本のリードでカバーされたか〜デプス〜を問題にする
このデプスがいくつになるかを、座位ごとに分布をとるときに、シンプルには、ポアソン分布を仮定する
ポアソン分布を仮定するということは、全座位の読まれ方が等確率であることを仮定している
実際には座位ごとに読まれ方の多寡があるので、現実データにより適合するのが負の二項分布である
以下のソースは、全部でN回の試行(最大数が有限)をして、そのうち、何回成功したか、として得られる二項分布と、平均何回成功するかを問題にしつつ、最大成功回数は無限大まで飛ばしたポアソン分布、ポアソン分布よりも分散が大きくなるように１パラメタを加えた負の二項分布をプロットする
平均成功回数をそろえてプロットしている

r <- 2
N.max <- 20
N <- 50
k <- 0:N
k.max <- 0:N.max
p <- 0.1
Np <- N*p
P.bin <- dbinom(k,N,p)
P.pois <- dpois(k.max,Np)
P.negbinom <- dnbinom(k.max,r,mu=Np)
ylim = range(c(P.bin,P.pois,P.negbinom))
plot(0:N.max,rep(0,N.max+1),type="l",ylim=ylim)
abline(v=N)
cex <- 2/log10(N.max+1)
points(k,P.bin,col=1,pch=20,cex=cex)
points(k.max,P.pois,col=2,pch=1,cex=cex)
points(k.max,P.negbinom,col=3,pch=20,cex=cex,type="h")
legend(N.max*0.8,ylim[2]*0.9, c("bin", "pois","negbinom"), col = c(1,2,3),pch = c(20,1,20))

Rstudioのshinyを使えば、インターラクティブにパラメタ変更できるプロットも描ける(shinyの使い方はこちら)

```{r, echo=FALSE}
N_maxs <- seq(from=0,to=100,by=10)
inputPanel(
  selectInput("N_max3", label = "Horizontal scale:",
              choices = N_maxs, selected = 20),
  sliderInput("N_adjust3", label = "Number of trials:",
              min = 0, max = max(N_maxs), value = 5, step = 1),
  
  sliderInput("p_adjust3", label = "Probability to occur:",
              min = 0, max = 1, value = 0.1, step = 0.05),
  sliderInput("negbin_param3", label = "Shape param of negbinom:",
              min = 0, max = max(N_maxs), value = 2, step = 0.05)
)



renderPlot({
  r <- input$negbin_param3
  N.max <- as.numeric(input$N_max3)
  N <- input$N_adjust3
  k <- 0:N
  k.max <- 0:N.max
  p <- input$p_adjust3
  Np <- N*p
  P.bin <- dbinom(k,N,p)
  P.pois <- dpois(k.max,Np)
  P.negbinom <- dnbinom(k.max,r,mu=Np)
  ylim = range(c(P.bin,P.pois,P.negbinom))
  plot(0:N.max,rep(0,N.max+1),type="l",ylim=ylim)
  abline(v=N)
  cex <- 2/log10(N.max+1)
  points(k,P.bin,col=1,pch=20,cex=cex)
  points(k.max,P.pois,col=2,pch=1,cex=cex)
  points(k.max,P.negbinom,col=3,pch=20,cex=cex,type="h")
  legend(N.max*0.9,ylim[2]*0.9, c("bin", "pois","negbinom"), col = c(1,2,3),
       pch = c(20,1,20))
})
```