次世代シークエンシングのためのポアソンモデル実習ハンドアウト

• Rmdファイルと、中から参照するshinyアプリの文書
• Rmdファイルと連携しているshinyアプリはこちら
• Rmdファイルからhtmlファイルやepubファイルを自作する方法はこちら
• htmlファイルはこちらにlogin as a guestして"handout"をクリック

Kindle

---
title: "poisson_eLearning"
author: "ryamada"
date: "Monday, February 02, 2015"
output: html_document
---

---
title: "poisson"
runtime: shiny
output: ioslides_presentation
---

## Poisson Processes and Poisson Distribution and Next-generation Sequencing ポアソン過程・ポアソン分布と次世代シークエンシング

## Introduction イントロダクション

Sample papers/documents with poisson processes
ポアソン過程を用いた論文/資料例

- ["Estimating Sequencing Coverage" Illumina](http://res.illumina.com/documents/products/technotes/technote_coverage_calculation.pdf)

- ["Rapid, high-throughput library preparation for next-generation sequencing" Nature methods](http://www.nature.com/nmeth/journal/v7/n8/full/nmeth.f.310.html)

<span style="color:green">**Assingment for e-learners e-ラーニング受講者用課題**</span>

- <span style="color:green">Describe briefly how poisson processe/distribution is used for what in the two documents listed above. 上記２文書の中でポアソン過程・分布が何を対象にどのように使われているかを簡単に記せ。</span>

------

## Contents コンテンツ

- 1 Binomial, poisson and negative-binomial distributions 二項分布・ポアソン分布・負の二項分布
- 2 NGS and poisson processes NGSとポアソン過程
- 3 Poisson processes and statistical tests ポアソン過程と検定

------

## Binomial distribution(1) 二項分布(1)

Formulas and general information on binomial distribution is [here, Wikipedia](http://en.wikipedia.org/wiki/Binomial_distribution).
二項分布に関する一般的事項(関数式を含む)に関する情報は[Wikipedia](http://en.wikipedia.org/wiki/Binomial_distribution)にて確認すること。

The probability of an event to happen is p per 1 trial.
You repeat N trials. 
Probability that A happens k times out of N, P.k.N, is calculated using binomial distribution as;

生起確率がpである試行をN回繰り返すとき、N回中k回起きる確率P.k.Nは二項分布で以下のように計算できる;

```{r}
N <- 5
k <- 2
p <- 0.1
P.k.n <- dbinom(k,N,p)
P.k.n
```
<span style="color:green">**Assingment for e-learners e-ラーニング受講者用課題**</span>

- <span style="color:green"> Calculate the probability of 5 successes out of 12 trials when success rate is 0.2. Report the codes of R to calculate the probability. 成功率0.2のときに12試行中5回成功する確率を計算せよ。Rでの計算コードとともに提出すること。</span>

Answer 解答
```{r,echo=FALSE}
N <- 12
k <- 5
p <- 0.1
dbinom(k,N,p)
```

------


## Binomial distribution(2) 二項分布(2)
Number of successes,k, is one of 0,1,...,N. You can calculate P.k.N for all k and plot the result as; 成功回数kは0,1,...,Nのいずれかになり、そのすべてを計算してプロットすれば：

```{r}
N <- 5
p <- 0.1
k <- 0:N
P.k.N <- dbinom(k,N,p)
plot(k,P.k.N,type="h",xlab="k",ylab="prob")
points(k,P.k.N,pch=20)
```
```{r,echo=FALSE}
P.k.N
```

This is called "binomial distribution". 二項分布です。

## Binomial Distribution(3) 二項分布(3)

Exercises 練習

[Click and try the first interactive plot](https://ryamada.shinyapps.io/shinyapps/poissonInteractive.Rmd)

<span style="color:green">**Assingment for e-learners e-ラーニング受講者用課題**</span>

- <span style="color:green">Change values of N and p in the interactive plot above and describe relation between "left/right symmetricity" and N and p. Nとpの値を変えて分布の形が変わることを確かめ、左右対称性とN,pとの関係について説明せよ</span>

- <span style="color:green">The sum of binomial probability with 0,1,...,N should be 1. Calculate it for N=10 and p=0.31. 二項確率を0,1,...,Nについて足し合わせると1になる。それをN=10, p=0.31のときについて確かめよ。</span>

- <span style="color:green">The mean of distribution is $\sum_{i=0}^{\infty} Pr(i) \times i$. Calculate the mean for N=10 and p=0.31.分布の平均は、$\sum_{i=0}^{\infty} Pr(i) \times i$と計算できる。N=10,p=0.31の場合について計算せよ。The following code would be helpful to calculate weighted sum of vector. ベクトルの重み付き和の計算は以下のコードを参考にせよ。</span>

```{r}
v <- c(1,2,3,4)
w <- c(0.1,0.2,0.3,0.4)
sum(v*w)
```
```{r,echo=FALSE}
N <- 10
p <- 0.31
k <- 0:N
P.k.N <- dbinom(k,N,p)
sum(P.k.N)
sum(P.k.N * k)
```

------

## Poisson distribution(1) ポアソン分布(1)

- Poisson distribution is somewhat similar to binomial distribution.
ポアソン分布は二項分布と似ているところがある。

- Poisson distribution focuses on 
the average number of occurence although the number of trials, N, and probability of occurence, p were by binomial distribution. 二項分布が総試行回数Nと生起確率pとを問題にしたのに対して、平均生起回数を問題にする

- No information on N is necessary for poisson distribution. 試行回数Nはわからなくてもポアソン分布は計算できる

------

## Poisson distribution(2) ポアソン分布(2)

Average number of occurence, Np. 平均生起回数 Np.

If the number of trials is N, then the probability to occur of one trial is $\frac{Np}{p}$.  
Probability of k out of N, follows binomial distribution.

N回の試行であるなら、１回の試行での生起確率は$\frac{Np}{N}$.
N回中k回起きる確率は二項分布に従う。

Poisson distribution gives you the probability of k without N or p. 
ポアソン分布は、Nとpなしにk回起きる確率を与える。


```{r}
Np <- 3.1 # Average number of occurence
N <- 5 # In case o binomial dist
p <- Np/N
N.max <- 10 # k can be infinit but set the max value for calculation purpose
k <- 0:N
k.max <- 0:N.max
P.bin <- dbinom(k,N,p) # binomial
P.pois <- dpois(k.max,Np)   # poisson
round(P.bin,4)
round(P.pois,4)
```

In case of binomial distribution with N=5, probability of k=6 is zero, but in case of poisson distribution with which k can be very large, probability of k=6 is positive. 二項分布ならば、N=5のときにk=6の確率は0だが、ポアソン分布の場合には、kの値に上限はなく、k=6の確率も0より大きい。



## Poisson distribution(3) ポアソン分布(3)

```{r,echo=FALSE,render=TRUE}
ylim <- c(0,1)
plot(0:N.max,rep(0,N.max+1),type="l",ylim=ylim,xlab="k",ylab="prob")
points(k,P.bin,col=1,pch=20,cex=2)
points(k.max,P.pois,col=2,pch=1,cex=2)
legend(9,0.4, c("bin", "pois"), col = c(1,2),
       pch = c(20,1))
```

No dots of binomial for k>5 with dots of pois for k>5. k>5に打点してあるのはポアソン分布であり、二項分布では打点がない。

-----

## Poisson Distribution(4) ポアソン分布(4)

Exercises 練習

<span style="color:green">**Assingment for e-learners e-ラーニング受講者用課題**</span>

- <span style="color:green">The sum of poisson distribution probability with 0,1,...,K(an finite value) should be less than 1. Calculate it for K=10 and p=0.31. Also calculate it for a large K, e.g., 100, that should be one. Report R codes with calculated values. ポアソン分布確率を0,1,...,K(有限値)について足し合わせると1より小さくなる。また、0,1,2,...K=100などと十分に大きくすると、和は1になる。それをK=10,100, p=0.31のときについて確かめよ。Rのコードと計算結果を報告せよ</span>

Answers 解答
```{r,echo=FALSE}
Np <- 3.1 # Average number of occurence
N <- 10
N.max <- 100
k <- 0:N
k.max <- 0:N.max
P.pois <- dpois(k,Np)   # poisson
P.pois.max <- dpois(k.max,Np)   # poisson
sum(P.pois)
sum(P.pois.max)
```

- <span style="color:green">The mean of distribution is $\sum_{i=0}^{\infty} Pr(i) \times i$. Calculate the mean for 0,1,2,...,K={10,100}, and p=0.31. Report R-codes and the calculated values. 分布の平均は、$\sum_{i=0}^{\infty} Pr(i) \times i$と計算できる。0,1,2,...,K={10,100},p=0.31の場合について計算せよ。Rのコードと計算結果を報告せよ</span>

Answers 解答

```{r}
sum(P.pois*k)
sum(P.pois.max*k.max)
```

- <span style="color:green">Random number generation. R's rpois() function generates random values from poisson distribution. Follow the commands below to generate them and calculate their mean and variance. Rのrpois()関数はポアソン分布乱数を生成する。次のコマンドを参考に、乱数を発生し、平均と分散を計算せよ。</span>

```{r}
n.rand <- 100
m <- 5 # mean of poisson dist
R <- rpois(n.rand,m)
```
Answers 解答(Your values should be close to the values below but should not be the same with them. 実行結果は以下の数値に近いはずだが、同一ではないはず)

```{r}
mean(R)
var(R)
```




## Poisson vs. Binomial distribution　ポアソン分布と二項分布の比較

Exercises 練習

[Click and try the second interactive plot](https://ryamada.shinyapps.io/shinyapps/poissonInteractive.Rmd)

<span style="color:green">**Assingment for e-learners e-ラーニング受講者用課題**</span>

- <span style="color:green">Change values of Average and Number of trials for binomial distribution in the interactive plot above and describe the items below.
 平均回数Averageと二項分布の場合の試行回数(Number of trials(binom))の値を変えて分布の形が変わることを確かめ、次のことについて説明せよ。
</span>

-- <span style="color:green">When Average = 3.1, poisson distribution does not change regardless of Number of trials for binom. However binomial distributions with number of trials(binom) 10 and 20 are different. Descrribe deviations of two binomial distributions from poisson distribution. 平均正規回数が3.1のときNumber of trials(binom)の値によらずポアソン分布は変わらないが、二項分布はNumber of trialsを10,20と変えると変化する。２つの二項分布のポアソン分布からのずれについて説明せよ</span>

-- <span style="color:green">Describe differece between poisson distribution and binomial distribution When average = 40 with number of trials(binom) = 100. 平均が40で二項分布の試行回数が100のときのポアソン分布と二項分布違いを説明せよ</span>

## Negative Binomial Distribution(1) 負の二項分布(1)

Poisson distribution is a good model for depth of next-generation sequensing (NGS), where probability to be read is equal throughout the genome. In reality the probability is not equal and negative binomial distribution fits better to depth of NGS than poisson distribution.

ポアソン分布を次世代シークエンサーのデプスの分布にあてはめることがあり、それは比較的よいあてはまりを示すが、その際には、全座位の読まれる確率が等しいこと仮定している。
しかしながら、実際には、座位によって読まれる確率にはばらつきがあり、それにより適合した分布が負の二項分布である。

All three distribution has the same average 3.1.
For binomial distribution, the number of trials is set as 50 and for negative binomial distributio, a parameter to regulate its shape, r, as 2. 
いずれも、平均3.1となるような二項分布・ポアソン分布・負の二項分布を合わせてプロットする。二項分布では総試行回数を50と仮定し、負の二項分布では形状を決めるパラメタrを2とした。

```{r}
Np <- 3.1 # Average number of occurence
N <- 50
N.max <- 20
k <- 0:N
k.max <- 0:N.max
p <- Np/N

r <- 2 # parameter for negative binomial dist

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.7,ylim[2]*0.9, c("bin", "pois","negbinom"), col = c(1,2,3),pch = c(20,1,20))
```

## Negative binomial distribution(2) 負の二項分布(2)

Cumulative distribution 累積分布

Consider you are interested in the fraction of values equal to or less than a threshold.　ある閾値以下の値の割合に興味があるとする。
more than a threshold. ある閾値以上の値に興味があるとする。

```{r}
Np <- 3.1
thres <- 2
N.max <- 10
k.max <- 0:N.max
P.pois <- dpois(k.max,Np)
P.pois
sum(P.pois[1:(thres+1)])
ppois(thres,Np) # Easy way
```
Consider you are interested in the fraction of values equal to or MORE than a threshold.　ある閾値以上の値の割合に興味があるとする。

```{r}
1-sum(P.pois[1:(thres)])
ppois(thres-1,Np,lower.tail=FALSE)
```


## Poisson vs. Negative Binomial distribution　ポアソン分布と負の二項分布の比較

Exercises 練習

[Click and try the third interactive plot](https://ryamada.shinyapps.io/shinyapps/poissonInteractive.Rmd)

<span style="color:green">**Assingment for e-learners e-ラーニング受講者用課題**</span>

- <span style="color:green">Change shape parameter of negative binomial distribution with fixed average.Describe how negative binomial distribution changes along with the increase of the parameter value. 平均値を固定し、負の二項分布の形状パラメタを変化させ、負の二項分布が形状パラメタの増加とともにどのように変化するか説明せよ
</span>

- <span style="color:green">Consider you are interested in the fraction with the value equal to or more than 25. When the average is 30, which of poisson and negative binomial distribution has higher fraction? When the average is 20, which has higher fraction? How do you describe the difference? 25以上の値をとる割合に興味があるとする。平均が20のとき、ポアソン分布と負の二項分布とではどちらが割合が高いか。平均が30のときはどうか。両者の違いをどのように説明するか？
</span>


## 2 NGS and poisson processes NGSとポアソン過程

## Total bps mapped マップされた総塩基数

<span style="color:green">**Assingment for e-learners e-ラーニング受講者用課題**</span>

- <span style="color:green">Consider the whole genome length is $G=10^5$ bp. Assume the average depth is $Np=25$. Calculate the sum of bps of all mapped reads. ゲノム全体で$G=10^5$塩基対と仮定する。今、平均デプス$Np=25$だと言う。マップされた塩基総数はいくつか
</span>

```{r,echo=FALSE}
G <- 10^5
Np <- 25
total.bp <- G*Np
```
- <span style="color:green">Assume poisson model, generate a random value set of depths for $G=10^5$ sites and draw its histogram and calculate its mean and sd. ポアソンモデルを仮定して、ゲノム全体($G=10^5$)のデプスを乱数発生し、ヒストグラムを描き、平均と標準偏差を計算せよ。
</span>

Your mean and sd should be close to the values below but should not be the same with them. 結果は以下の値に近いが同じではない

```{r,echo=FALSE}
R <- rpois(G,Np)
h <- hist(R)
mean(R)
sd(R)
```
- <span style="color:green">Assume the mean and sd of your real data were 25 and 8, that was larger than the variance based on poisson model. Generate random values in negative binomial distribution with various  shape parameters and obtain a simulational sample set with mean 25 and sd somewhere between 7.5 and 8.5. Smaller is the shape parameter, larger the sd.
実観測データでは平均が25で分散が8だったと言う。負の二項分布から乱数を発生させ、平均が25で標準偏差が7.9から8.1の間の値となるようなデプスの値のセットを作成せよ。負の二項分布の形状パラメタは小さいほど標準偏差は大きくなる。
</span>

Example of shape parameter = 0.8 形状パラメタr=0.8の例
```{r}
G <- 10^5
Np <- 25
r <- 0.8
R.sim <- rnbinom(G,r,mu=Np)
mean(R.sim)
sd(R.sim)
hist(R.sim)
```

```{r,echo=FALSE}
r <-　16
R.sim <- rnbinom(G,r,mu=Np)
hist(R.sim)
mean(R.sim)
sd(R.sim)
```

## Increase total read-bases 全リード塩基数を増やす

When average depth = 30, the fraction of depth >= 25 is around 0.85 under the poisson distribution model.
平均デプスが30であるとき、ポアソン分布モデルであれば、デプスが25塩基以上である割合は85%ほどである。

```{r}
Np <- 30
thres <- 25
ppois(thres-1,Np,lower.tail=FALSE)
```

However under the negative binomial distribution model with r = 16, it is about 0.75. しかしながら、負の二項分布(r=16)の場合には70%ほどである。

```{r}
r <- 16
pnbinom(thres-1,r,mu=Np,lower.tail=FALSE)
```
Assume the total read bases and average depth are proportinal and the shape parameter of negative binomial distribution does not change with the increase of total bases, then, 全リード塩基数を増やすとデプス平均が増え、負の二項分布の形状パラメタを変えないとすれば…

<span style="color:green">**Assingment for e-learners e-ラーニング受講者用課題**</span>

- <span style="color:green">Increase the total read bases as $G \times Np \times x$, so that the fraction of depth >= 25 is 0.85 and answer x. デプスが25以上である割合が0.85となるように、リード総塩基数を$G\times Np \times x$のように増やすことにする。xの値を答えよ</span>

Example 例 x=2 -> Too good これは良すぎる
```{r}
Np <- 30
r <- 16
thres <- 25
x <- 2
pnbinom(thres-1,r,mu=Np*x,lower.tail=FALSE)
```

Answer value will return alsomst 0.85 これくらい近くできる
```{r,echo=FALSE}
Np <- 30
r <- 16
thres <- 25
x <- 35.36/30
pnbinom(thres-1,r,mu=Np*x,lower.tail=FALSE)
```

- <span style="color:green">You may find the answer with the fourth plot in [the Interactive plot](https://ryamada.shinyapps.io/shinyapps/poissonInteractive.Rmd)　インターラクティブプロットの４番目を使って求めることもできる</span>

## Mutations and Polymorphisms 変異と多型

Some loci are heterogeneous due to gametic or somatic mutations. 遺伝子多型のためや体細胞変異のために、一様でない場所がある。

In case of SNP, fractions of two alleles are 0.5 and 0.5. In case of cancer cell populations, number of alleles can be two or more and the fractions vary.
遺伝子多型の場合は、その割合は0.5, 0.5 であり、癌細胞集団などの場合には、2またはそれより多くのアレルがあり、その割合はまちまちである。

We handle biallelic cases with the fraction u and 1-u. ここでは2アレルの場合を扱い、その割合をu,1-uとする。

We adopted negative binomial distribution for depth because of heterogeneity in readability among genomic loci. Under the assumption two alleles at the same loci have the same readability, the number of reads per alleles at particular loci is in binomial distribution.

ゲノム上の場所により読まれやすさに違いがあるためにデプスを負の二項分布で近似したが、同一座位の2アレルの読まれやすさは等しいとみなすならば、各座位のデプスの2アレルへの振り分けは二項分布になる。

SNPであれば、u=0.5,1-u=0.5での二項分布、それ以外であれば、u,(1-u)に応じた二項分布となる。

```{r,echo=FALSE}
# 二項分布になることの確認
#n <- 10^6
#m <- 10
#R1 <- rpois(n,m)
#R2 <- rpois(n,m)
#R12 <- R1+R2
#s <- list()
#for(i in 0:max(R12)){
#  s[[i+1]] <- which(R12==i)
#}

#par(ask=TRUE)

#for(i in 2:length(s)){
#	if(length(s[[i]])>0){
#	tab <- c(tabulate(R1[s[[i]]]+1,i))
#	d <- dbinom(0:(i-1),i-1,0.5)
#	a <- (tab/sum(tab))
#	
#	plot(a,d)
#	abline(0,1,col=2)
#	}
#
#}
```

## Binomial Distribution Revisited もう一度二項分布

Assume a SNP locus with depth =30. Draw probability of reads of one allele. デプスが30のSNP座位を考える。片方のアレルの本数が何本になるかを描いてみる。

```{r}
d <- 30
k <- 0:d
u <- 0.5
P.binom <- dbinom(k,d,u)
plot(k,P.binom,type="h")
```


Probability to have more than 10 reads for both alleles is ; どちらのアレルのリード数も10より大きくなる確率は;

```{r}
thres <- 10
k2 <- d-k
k.min <- apply(cbind(k,k2),1,min)
s <- which(k.min>thres)
col <- rep(1,d+1)
col[s] <- 2
plot(k,P.binom,type="h",col=col)
points(k,P.binom,pch=20,col=col)
sum(P.binom[s])
```

Probability to have loci with the ratio of two alleles are less than 0.75 or more than 1/0.75 is; 2アレルのリード数の比が0.75未満または1/0.75より大の確率は;

```{r}
frac.A <- k/d
frac.B <- 1-frac.A
thres <- 0.75
frac.ratio <- frac.A/frac.B
s <- which(frac.ratio < thres | frac.ratio > 1/thres)
col <- rep(1,d+1)
col[s] <- 2
plot(k,P.binom,type="h",col=col)
points(k,P.binom,pch=20,col=col)
sum(P.binom[s])
```

## Two Alleles May Not Be EVEN 2アレルは不均等かもしれない

Above we assumed readability of two alleles are equal, but allele-itself, something close to allele might affect on readbility and u might deviate from 0.5. Also copy number of chromosomes in sample may somewhat deviate from eveness. 上の例では、2アレルの読まれ方が均等であることを仮定したが、アレル自体やアレル近傍の影響でuが0.5からずれるかもしれない。また、サンプル調整の段階でコピー数に均等からのずれが生じるかもしれない。

The following assumes normal distribution around 0.5 for u and its effect on the probability to have more than 10 reads for both alleles with Monte-carlo.
以下では、uに0.5を中心とする正規分布を仮定し、それが、「両アレルとも10リードより多くなる確率がどのように影響されるかをシミュレーショナルに調べる

Distribution of u. uの分布
```{r}
n.iter <- 10000
us <- rnorm(n.iter,0.5,0.05)
hist(us)
```

Distribution of probability with >10 reads for both alleles. 両アレル >10 readsの確率の分布
```{r}
res <- rep(NA,n.iter)
d <- 30
k <- 0:d
thres <- 10
k2 <- d-k
k.min <- apply(cbind(k,k2),1,min)
s <- which(k.min>thres)
for(i in 1:n.iter){
  u <- us[i]
  P.binom <- dbinom(k,d,u)
  res[i] <- sum(P.binom[s])
}
hist(res)
```

## Combine negative binomial and binomial distributions 負の二項分布と二項分布とを組み合わせる

For a given depth we calculated the probability to have more than 10 reads for both alleles of a SNP. We calculated in two ways; (i) u=0.5 and fixed and (ii) u ~ Normal(0.5). あるデプスにおいて、両アレルが10リードより多く観測される確率を計算した。二つの方法(i) uを0.5に固定、(ii) u~Normal(0.5)で実施した。

The depth vary in the negative binomial distribution. Now we calculate the probability to have more than 10 with the negative binomial distribution. デプスは負の二項分布に従って変わるので、それを考慮してこの確率を計算することにする。

Average depth = 30, shape parameter =16
```{r}
thres <- 10
n.iter <- 10^5
Np <- 30
r <- 16
r.depth.pois <- rpois(n.iter,Np)
r.depth.nb <- rnbinom(n.iter,r,mu=Np)
us <- rnorm(n.iter,0.5,0.05) # same as above
res.pois <- res.nb <- rep(NA,n.iter)
for(i in 1:n.iter){
  d <- r.depth.pois[i]
  k <- 0:d
  k2 <- d-k
  k.min <- apply(cbind(k,k2),1,min)
  s <- which(k.min>thres)
  u <- us[i]
  P.binom <- dbinom(k,d,u)
  res.pois[i] <- sum(P.binom[s])
  d <- r.depth.nb[i]
  k <- 0:d
  k2 <- d-k
  k.min <- apply(cbind(k,k2),1,min)
  s <- which(k.min>thres)
  u <- us[i]
  P.binom <- dbinom(k,d,u)
  res.nb[i] <- sum(P.binom[s])
}
h <- hist(c(res.pois,res.nb),plot=FALSE)
hist(res.pois,breaks=h$breaks,density=17,ylim=c(0,max(h$counts)))
hist(res.nb,breaks=h$breaks,add=TRUE,col=2,density=15)
legend(0,max(h$counts)/2,c("pois","negbinom"),text.col=c(1,2))
```

Some fraction with 0 probability because sum of depths of both alleles is low. 両方のアレルのデプスの和が小さい場合に確率が0になる

```{r}
matplot(cbind(sort(res.pois),sort(res.nb)),type="l",ylim=c(0,1))
legend(n.iter*0.8,0.3,c("pois","negbinom"),text.col=c(1,2))
```

<span style="color:green">**Assingment for e-learners e-ラーニング受講者用課題**</span>

- <span style="color:green">Consider a loci of somatic mutation in cancer cells, where the fraction of mutant chromosomes is assumed u = 0.1. Answer the depth ncessary to have both the mutant and wild reads more than 5 with the probabililty more than 0.8. がん細胞サンプルにおいて、変異染色体の割合が0.1と予想したローカスがある。変異型リード数が5より大きくなる確率が0.8より大きくなるために必要なデプスを答えよ</span>

-- <span style="color:green">First generate binomial distribution with u=0.1 for a depth of your choice. 適当にデプスを定め、u=0.1での二項分布をまず作成せよ</span>

Example 例 depth = 30

```{r,echo=FALSE}
d <- 30
k <- 0:d
u <- 0.1
P.binom <- dbinom(k,d,u)
plot(k,P.binom,type="h")
```

-- <span style="color:green">Next calculate the probability to have more than 5 reads for both alleles.　次に両方のアレルのリード数が5より大きくなる確率を計算せよ</span>

Example 例 depth =30 Too low 低すぎる

```{r,echo=FALSE}
thres <- 5
k2 <- d-k
k.min <- apply(cbind(k,k2),1,min)
s <- which(k.min>thres)
col <- rep(1,d+1)
col[s] <- 2
plot(k,P.binom,type="h",col=col)
points(k,P.binom,pch=20,col=col)
sum(P.binom[s])
```

-- <span style="color:green">Change depth value as you like. デプスの値を適当に変えてみる</span>

```{r,echo=FALSE}
d <- 150
k <- 0:d
u <- 0.1
P.binom <- dbinom(k,d,u)
thres <- 5
k2 <- d-k
k.min <- apply(cbind(k,k2),1,min)
s <- which(k.min>thres)
col <- rep(1,d+1)
col[s] <- 2
plot(k,P.binom,type="h",col=col)
points(k,P.binom,pch=20,col=col)
sum(P.binom[s])
```

## Error calls エラーコール

Error call of bases are unavoidable. 塩基のコールにエラーはつきもの

Assume depth is in negative binomial distribution and base call errors happen with the fixed probability, p.err, throughout the genome. リードデプスが負の二項分布に従い、ベースコールエラーがゲノム全体で同一の確率 p.err で起きるとする。

## Two binomial distributions of variants and errors

For a given depth, number of error alleles is in binomial distribution.
特定のデプスに関しては、エラーアレル本数は二項分布

Consider depth = 30 and u for a variant is 0.5 and u.err for an error is 0.01.
デプス=30とし、バリアントのuを0.5、エラーのそれをu.err = 0.01とする

```{r}
u <- 0.5
u.err <- 0.01
d <- 30
k <- 0:d
P.binom <- dbinom(k,d,u)
P.binom.err <- dbinom(k,d,u.err)
matplot(k,cbind(P.binom,P.binom.err),type="h",xlab="k",ylab="prob")
legend(d*0.7,max(c(P.binom,P.binom.err)*0.8),c("variant","error"),text.col=c(1,2))
```

## Whole picture of two binomial distibutions for various depths 異なるデプスの全体について、二つの二項分布を眺める


```{r,echo=FALSE}
plot2Ddist <- function(mat,r=5,theta=20,phi=35,shade=0.5,xlab="depth",ylab="No. errors",zlab="probability",main=""){
  persp(z=mat[,rep(seq(ncol(mat)), each=2)], r=r, theta=theta, phi=phi, shade=shade,border=NULL, col=rep(c("#808080FE","#00000000"), each=nrow(mat)-1),xlab=xlab,ylab=ylab,zlab=zlab,main=main)
}
```

```{r}
u <- 0.5
u.err <- 0.01 # p.err
Np <- 30
r <- 16
k.max <- 0:60
P.negbinom <- dnbinom(k.max,r,mu=Np)
res <- res.err <- matrix(NA,length(k.max),length(k.max))
for(i in 1:length(k.max)){
  P.binom <- dbinom(0:(i-1),i-1,u) * P.negbinom[i]
  P.binom.err <- dbinom(0:(i-1),i-1,u.err) * P.negbinom[i]
  res[i,1:length(P.binom)] <- P.binom
  res.err[i,1:length(P.binom.err)] <- P.binom.err
}
no.err <- apply(res.err,2,sum)
```

The followings show the depth in one axis and the number of variant allele (or error call allele) in the other axis.

The peak for variants is at the mode of depth distribution and at the half of depth.
The peak for errors is at the mode of depth and 0 of depth because of low error rate.

以下の図はデプスの分布とバリアントアレル(またはエラーアレル)の本数を描いたもの。
バリアントのピークはデプス分布の最大値でデプスの半分のところ。
エラーの場合は、デプス分布の最大値で、0のところ。エラー率が低いから。

```{r}
par(mfcol=c(1,2))
persp(k.max,k.max,res,phi=20,theta=50,xlab="depth",ylab="no.allele",main="variant",shade=0.3,col=rgb(0.0,1,1))
persp(k.max,k.max,res.err,phi=20,theta=50,xlab="depth",ylab="no.allele",main="error",shade=0.3,col=rgb(0.0,1,1))
```
```{r,echo=FALSE}
par(mfcol=c(1,1))
```

From another angle.
別のアングルから見たところ。

```{r}
par(mfcol=c(1,2))
persp(k.max,k.max,res,phi=20,theta=150,xlab="depth",ylab="no.allele",main="variant",shade=0.3,col=rgb(0.0,1,1))
persp(k.max,k.max,res.err,phi=20,theta=150,xlab="depth",ylab="no.allele",main="error",shade=0.3,col=rgb(0.0,1,1))
```
```{r,echo=FALSE}
par(mfcol=c(1,1))
```


## Variants are Hidden in Errors バリアントはエラーの中に隠れる

Throughout the genome, onle few sites are variants, therefore the distribution of variant/error alleles looks like below.

ゲノム全体ではバリアントサイトは少ないので、バリアントとエラーとを併せた分布は以下のようになる

Find a peak of variants. バリアントのピークはどこにあるか？

```{r}
par(mfcol=c(1,2))
frac <- 0.1
res.comb <- res*frac+ res.err*(1-frac)
persp(k.max,k.max,res.comb,phi=20,theta=50,xlab="depth",ylab="no.allele",main="whole",shade=0.3,col=rgb(0.0,1,1))

persp(k.max,k.max,res.comb,phi=20,theta=150,xlab="depth",ylab="no.allele",main="whole from another angle",shade=0.3,col=rgb(0.0,1,1))
```
```{r,echo=FALSE}
par(mfcol=c(1,1))
```

## Identify Variants Among Errors エラーの中からバリアントを見つける

When we find variants from a pile of errors, we compare two binomial distributions of variants and errors for given depths.
バリアントをエラーの山から見つけるとき、デプスごとにバリアントとエラーが作る二項分布を比較する。

Assume depth = 30, u =0.5 for variants, u.err =0.01 for errors, fraction of variant sites among genome f =0.01.

デプス=30,u=0.5 (バリアント), u.err = 0.01 (エラー)、ゲノム全体でのバリアント箇所の割合f=0.01とする。


```{r}
u <- 0.5
u.err <- 0.01
f <- 0.01
d <- 30
k <- 0:d
P.binom.w <- dbinom(k,d,u)*f
P.binom.err.w <- dbinom(k,d,u.err)*(1-f)
matplot(k,cbind(P.binom.w,P.binom.err.w),type="h",xlab="k",ylab="prob")
legend(d*0.7,max(c(P.binom.w,P.binom.err.w)*0.8),c("variant","error"),text.col=c(1,2))
```

## Likelihood ratio 尤度比

Assume depth = 30 and k reads fit to a variant allele. The likelihood of this observation based on the hypothesis that the locus is variant and another likelihood of this observation based on the hypothesis that call error made this are to be calculated. Their ratio is likelihood ratio.

デプスが30で、バリアントアレルと思われるアレルリードがk本であったとき、バリアントであったがためにk本のアレルリードが得られた尤度と、コールエラーのためにk本の「似せのアレルリードが得られた尤度」との比をとる

Likelihood ratio is the ratio of probabilities based on tow binomial distributions.
尤度比は二項分布から計算した確率の比である。

The fraction of variants among genome is a prior. ゲノム全体におけるバリアントの割合は事前確率である。

```{r}
par(mfcol=c(1,2))
matplot(k,cbind(P.binom.w,P.binom.err.w),type="h",xlab="k",ylab="prob")
legend(d*0.5,max(c(P.binom.w,P.binom.err.w)*0.8),c("variant","error"),text.col=c(1,2))
log.like.ratio <- log10(P.binom.w/P.binom.err.w)
plot(k,log.like.ratio)
abline(h=2,col=2)
legend(d*0.5,2-0.5,"x100")
par(mfcol=c(1,1))
```

<span style="color:green">**Assingment for e-learners e-ラーニング受講者用課題**</span>

[Click and try the fifth interactive plot](https://ryamada.shinyapps.io/shinyapps/poissonInteractive.Rmd)


- <span style="color:green">You can change depth, variant's u, mean error rate, fraction of variants. "no of minor allele reads to be focused" change the y-axsi scale of the right panel. The likelihood ratio is the ratio of black/red bars in the right panel. Change the parameters as you like and see the ratio and calculate the likelihood ratio with R yourself.
デプス、バリアントのu、エラー率、ゲノムにおけるバリアントの割合を変えることができる。"no of minor allele reads to be focused"は右のパネルのy軸のスケールを変える。尤度比は、黒と赤のバーの長さの比である。自分でパラメタを適当に設定し、尤度比を目で確認せよ、そのうえで、Rでそれを計算せよ</span>

## Ditecting novel variants 新規バリアント検出

When we calculated log likelihood in the previous section, we KNEW the variant allele, which might be the case of know variant sites.

前の例で対数尤度比を計算したとき、バリアントアレルがどちらであるかを『知って』いた。それは既知バリアントサイトなどの場合である。

When we are to detect novel variants, we should calculate likelihood ratio for "allele observed fewer".

新規バリアントの検出の場合には、『観測されたアレルのうち少ない方が』という計算をする必要がある。

```{r}
u <- 0.5
u.err <- 0.01
f <- 0.01
d <- 30
k <- 0:d


P.binom.w <- dbinom(k,d,u)*f
P.binom.err.w <- dbinom(k,d,u.err)*(1-f)

P.binom.w.2 <- P.binom.w + P.binom.w[(d+1):1]
P.binom.err.w.2 <- P.binom.err.w + P.binom.err.w[(d+1):1]

k.turn <- ceiling(d/2)+1
k.minor <- k[1:k.turn]

P.binom.w.2 <- P.binom.w.2[1:k.turn]
P.binom.err.w.2 <- P.binom.err.w.2[1:k.turn]

log.like.ratio <- log10(P.binom.w/P.binom.err.w)
log.like.ratio.2 <- log10(P.binom.w.2/P.binom.err.w.2)
```

```{r}
matplot(k.minor,cbind(P.binom.w.2,P.binom.err.w.2),type="h",xlab="k",ylab="prob")
legend(d*0.7,max(c(P.binom.w.2,P.binom.err.w.2)*0.8),c("variant","error"),text.col=c(1,2))
```


```{r}
par(mfcol=c(1,2))
matplot(k.minor,cbind(P.binom.w.2,P.binom.err.w.2),type="h",xlab="k",ylab="prob")
legend(k.turn*0.5,max(c(P.binom.w.2,P.binom.err.w.2)*0.8),c("variant","error"),text.col=c(1,2))

plot(k.minor,log.like.ratio.2)
abline(h=2,col=2)
legend(k.turn*0.5,2-0.5,"x100")
par(mfcol=c(1,1))
```

## Comparison of know and novel variants 既知・新規バリアントの比較

Change conditions
条件を変えてみよう

```{r}
u <- 0.3
u.err <- 0.03
f <- 0.05
d <- 7
k <- 0:d


P.binom.w <- dbinom(k,d,u)*f
P.binom.err.w <- dbinom(k,d,u.err)*(1-f)

P.binom.w.2 <- P.binom.w + P.binom.w[(d+1):1]
P.binom.err.w.2 <- P.binom.err.w + P.binom.err.w[(d+1):1]

k.turn <- ceiling(d/2)+1
k.minor <- k[1:k.turn]

P.binom.w.2 <- P.binom.w.2[1:k.turn]
P.binom.err.w.2 <- P.binom.err.w.2[1:k.turn]


log.like.ratio <- log10(P.binom.w/P.binom.err.w)
log.like.ratio.2 <- log10(P.binom.w.2/P.binom.err.w.2)
```

```{r}
par(mfcol=c(1,2))
matplot(k,cbind(P.binom.w,P.binom.err.w),type="h",xlab="k",ylab="prob",main="known")
legend(k.turn*0.5,max(c(P.binom.w.2,P.binom.err.w.2)*0.8),c("variant","error"),text.col=c(1,2))
matplot(k.minor,cbind(P.binom.w.2,P.binom.err.w.2),type="h",xlab="k",ylab="prob",main="novel")
legend(k.turn*0.5,max(c(P.binom.w.2,P.binom.err.w.2)*0.8),c("variant","error"),text.col=c(1,2))
```


```{r}
plot(k,log.like.ratio,main="known vs. novel",col=1)
points(k.minor,log.like.ratio.2,col=2,pch=20)
abline(h=2,col=2)
legend(k.turn*0.5,min(log.like.ratio)+1,c("known","novel"),text.col=c(1,2))
```

• shiny文書

---
title: "binom-poisson-negbinom"
author: "ryamada"
date: "Monday, February 02, 2015"
output: html_document
runtime: shiny
---

## Binomial distribution

```{r, echo=FALSE}
inputPanel(
  sliderInput("N_adjust", label = "Number of trials:",
              min = 0, max = 50, value = 5, step = 1),
  
  sliderInput("p_adjust", label = "Probability to occur:",
              min = 0, max = 1, value = 0.1, step = 0.05)
)

renderPlot({
  N <- input$N_adjust
  k <- 0:N
  p <- input$p_adjust
  plot(k,dbinom(k,N,p),type="h")
})
```

## Binomial and poisson distributions

Number of trials(binom) should be more than Average.

```{r, echo=FALSE}
N_maxs <- seq(from=0,to=100,by=10)
inputPanel(
  selectInput("N_max", label = "Horizontal scale:",
              choices = N_maxs, selected = 20),
  sliderInput("Np_adjust2", label = "Average:",
              min = 0, max = max(N_maxs), value = 3.1, step = 0.01),
  sliderInput("N_adjust2", label = "Number of trials(binom):",
              min = 0, max = max(N_maxs), value = 10, step = 1)  

)


renderPlot({
  N.max <- as.numeric(input$N_max)
  N <- input$N_adjust2
  Np <- input$Np_adjust2
  k <- 0:N
  k.max <- 0:N.max
  p <- Np/N
  
  P.bin <- dbinom(k,N,p)
  P.pois <- dpois(k.max,Np)
  ylim = range(c(P.bin,P.pois))
  plot(0:N.max,rep(0,N.max+1),type="l",ylim=ylim,xlab="k",ylab="prob")
  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)
  legend(N.max*0.9,ylim[2]*0.9, c("bin", "pois"), col = c(1,2),
       pch = c(20,1))
})
```

## Binomial, poisson and negative-binominal distributions.

Number of trials(binom) should be more than Average.

```{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("Np_adjust3", label = "Average:",
              min = 0, max = max(N_maxs), value = 3.1, step = 0.01),
  sliderInput("N_adjust3", label = "Number of trials:",
              min = 0, max = max(N_maxs), value = 5, step = 1),

  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
  Np <- input$Np_adjust3
  k <- 0:N
  k.max <- 0:N.max
  p <- Np/N

  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))
})
```

## Poisson and negative-binominal cumulative distributions.

Fraction of equal to or more than threshold is indicated.
```{r, echo=FALSE}
N_maxs <- seq(from=0,to=100,by=10)
inputPanel(
  selectInput("N_max4", label = "Horizontal scale:",
              choices = N_maxs, selected = 20),
  sliderInput("Np_adjust4", label = "Average:",
              min = 0, max = max(N_maxs), value = 3.1, step = 0.01),

  sliderInput("negbin_param4", label = "Shape param of negbinom:",
              min = 0, max = max(N_maxs), value = 2, step = 0.05),
  sliderInput("thres4", label = "Threshold:",
              min = 0, max = max(N_maxs), value = 5, step = 1)
)



renderPlot({
  r <- input$negbin_param4
  N.max <- as.numeric(input$N_max4)
  Np <- input$Np_adjust4
  k.max <- 0:N.max
  thres <- input$thres4
  P.pois <- dpois(k.max,Np)
  P.negbinom <- dnbinom(k.max,r,mu=Np)
  Q.pois <- ppois(k.max,Np)
  Q.negbinom <- pnbinom(k.max,r,mu=Np)

  ylim = c(-0.2,1)
  plot(0:N.max,rep(0,N.max+1),type="l",ylim=ylim,xlab="k",ylab="cummulative prob")
  abline(v=thres-0.5)
  cex <- 2/log10(N.max+1)
  points(k.max,Q.pois,col=2,pch=1,cex=cex,type="s")
  points(k.max,Q.negbinom,col=3,pch=20,cex=cex,type="s")
  legend(N.max*0.9,ylim[2]*0.3, c("pois","negbinom"), col = c(2,3),
       lty = c(1,1))
  q.pois <- ppois(thres-1,Np,lower.tail=FALSE)
  q.negbinom <- pnbinom(thres-1,r,mu=Np,lower.tail=FALSE)
  y.leg <- c(max(1-q.pois,1-q.negbinom))
  if(q.pois > q.negbinom){
    y.leg <- c(y.leg-0.05,y.leg)
  }else{
    y.leg <- c(y.leg,y.leg-0.05)
  }
  legend(thres +1,y.leg[1]-0.1,round(q.pois,3),text.col=2,bty="n")
  legend(thres +1,max(y.leg)-0.05,"frac >= thres",text.col=1,bty="n")
  legend(thres +1,y.leg[2]-0.1,round(q.negbinom,3),text.col=3,bty="n")

})
```


## Binomial vs. Binomial 


```{r, echo=FALSE}

inputPanel(
  sliderInput("N5", label = "N or depth:",
              min = 0, max = 100, value = 20, step = 1),

  sliderInput("m1", label = "mean 1 for variant:",
              min = 0, max = 1, value = 0.5, step = 0.01),
  sliderInput("m2", label = "mean 2 for error:",
              min = 0, max = 1, value = 0.01, step = 0.001),
   sliderInput("frVar", label = "fraction of variants:",
              min = 0, max = 0.2, value = 0.01, step = 0.001),
  sliderInput("ylim", label = "no of minor allele reads to be focused:",
              min = 0, max = 100, value = 4, step = 1)

)



renderPlot({
  N <- input$N5
  k <- 0:N
  m1 <- input$m1
  m2 <- input$m2
  d1 <- dbinom(k,N,m1)
  d2 <- dbinom(k,N,m2)
  f <- input$frVar
  
  par(mfrow=c(1,2))
  plot(k,d1,type="h",xlab="no.minor allele reads",ylab="prob",ylim=c(0,1),main="binom dists")
  points(k+0.1,d2,type="h",col=2)
  legend(N*0.8,0.8,c("variant","error"),text.col=c(1,2))
  kid <- input$ylim
  ylim <- c(0,max(d1[kid+1]*f,d2[kid+1]*(1-f))*1.2)
  plot(k,d1*f,type="h",xlab="no.minor allele reads",ylab="prob",ylim=ylim,main="frac of var considered")
  points(k+0.1,d2*(1-f),type="h",col=2)
  legend(N*0.8,0.8,c("variant","error"),text.col=c(1,2))
  par(mfcol=c(1,1))
})
```