2015-02-18

DNA鑑定にポアソン対数正規分布

昨日の記事や一昨日の記事に対数正規分布・ポアソン対数正規分布のことをメモした
常染色体上のSTRをPCR増幅をかませた実験によってピークシグナル型データを取得し、それに基づいてDNA型タイピングをするやり方がある
そのデータについて、対数正規分布・ポアソン対数正規分布を応用してみる話としてまとめてみる
話の骨格
- このSTR-DNA型タイピングは、ベースラインにノイズがあり、真のシグナルピークがその中に立ち上がる形をとる
- 真のシグナルは、DNA分子からPCR増幅によって指数関数的に増幅させて形成されるもので、蛍光を光学的検出したものである
- ノイズは、「種」がないところから、光学系によくあるショット雑音として拾われるものである(と考えられている)
- ノイズも真のシグナルも正の値のみを取り、対数正規分布様の形を取るらしいことは実験データ的に知られている。きれいに対数正規分布であるかどうかは、はてな
- PCR増幅を使っているので、対数正規分布的になることは妥当では、ある
- ポアソン的かどうかははっきりしないが、すくなくともqRT-PCRのように、特定の遺伝子のmRNA分子は「存在０」であったりする場合のことなどを考えると、「0を含めた、0,1,2,...」を扱えることは適切っぽく、それが理由で、単なる対数正規分布ではなく、ポアソン対数正規分布を使うメリットはあるらしい
勉強会(2015/04/25)のための資料として以下のshinyを作成
- 以下の囲みファイルを"hoge.Rmd"として保存し、Rstudioで開き、shinyアプリ扱いにすると、インターラクティブプロットが現れます(文字のエンコードに注意が必要かもしれません)(Rstudio-shinyサイト)

---
title: "ポアソン対数正規分布"
author: "ryamada"
date: "Monday, February 16, 2015"
output: html_document
runtime: shiny
---

# ノイズ
```{r,echo=FALSE}
my.dlnorm <- function(x,mea=1,mod=0.5){
  S <- sqrt(2/3*log(mea/mod))
	M <- log(mod) + S^2
	dlnorm(x,M,S)
}
my.plnorm <- function(q,mea=1,mod=0.5,lower.tail=TRUE){
  S <- sqrt(2/3*log(mea/mod))
	M <- log(mod) + S^2
	plnorm(q,M,S,lower.tail=lower.tail)
}
my.qlnorm <- function(p,mea=1,mod=0.5,lower.tail=TRUE){
  S <- sqrt(2/3*log(mea/mod))
	M <- log(mod) + S^2
	qlnorm(p,M,S,lower.tail=lower.tail)
}
my.rlnorm <- function(n,mea=1,mod=0.5){
  S <- sqrt(2/3*log(mea/mod))
	M <- log(mod) + S^2
	rlnorm(n,M,S)
}
library(poilog)
my.dpoislnorm <- function(x,mea=1,mod=0.5){
  S <- sqrt(2/3*log(mea/mod))
	M <- log(mod) + S^2
	dpoilog(x,M,S)
}

my.rpoislnorm <- function(n,mea=1,mod=0.5){
	S <- sqrt(2/3*log(mea/mod))
	M <- log(mod) + S^2
	rpoilog(n,M,S,keep0=TRUE)
}
```

正規分布ノイズ
```{r}
n <- 10^3
r1 <- rnorm(n)
plot(r1,type="h",ylim = c(min(r1)-1,max(r1)+20))
```

対数正規分布ノイズ
```{r}
n <- 10^3
r2 <- my.rlnorm(n)
plot(r2,type="h",ylim=c(0,max(r2)+20))
```

RFU noise : 
Analytical Thresholds and Sensitivity: Establishing RFU Thresholds for Forensic DNA Analysis
J Forensic Sci, January 2013, Vol. 58, No. 1
doi: 10.1111/1556-4029.12008
http://www.az-forensics.com/docs/pdfs/Analytical-Thresholds-and-Sensitivity.pdf p.123

# Log-normal distirbution 対数正規分布

```{r, echo=FALSE}
inputPanel(
  sliderInput("mod2", label = "mode:",
              min = 0, max = 2000, value = 500, step = 1)
,
sliderInput("mean2", label = "mean:",
              min = 0, max = 10000, value = 2000, step = 1)
,
sliderInput("xmax2", label = "xmax:",
              min = 0, max = 10000, value = 3000, step = 1)
)

renderPlot({
xmax <- input$xmax2
mea <- input$mean2
mod <- input$mod2
med <- (mea^2*mod)^(1/3)

S <- sqrt(2/3*log(mea/mod))
M <- log(mod) + S^2

t <- seq(from=0,to=xmax,length=101)
d <- dlnorm(t,M,S)
plot(t,d,type="l")
abline(v=mea,col=2)
abline(v=mod,col=3)
abline(v=med,col=4)

})
```

# 対数正規分布の考え方

[Log-normal Distributions
across the Sciences:
Keys and Clues Bioscience. 2001; 51:341-352](http://stat.ethz.ch/~stahel/lognormal/bioscience.pdf)

p. 343



# Shot noise ショット・ノイズ

光刺激の離散性と関連するノイズ

[ショット雑音](http://en.wikipedia.org/wiki/Shot_noise)

Improving signal-to-noise ratio in fluorescence detection for medical purposes

http://essay.utwente.nl/62160/1/BSc_S_van_Binsbergens_en_R_Fentrop.pdf 

離散性があって、でもバラつきがある。

# ポアソン分布

小さい期待値でのポアソン分布と
大きい値でのポアソン分布とでは

バラツキの大きさと真のシグナルとの比(ノイズ・シグナル比)が問題になる。

```{r, echo=FALSE}
inputPanel(
  sliderInput("mean3", label = "mean:",
              min = 0, max = 100, value = 5, step = 0.01)

,
sliderInput("xmax3", label = "xmax:",
              min = 0, max = 1000, value = 30, step = 1)
)

renderPlot({
xmax <- input$xmax3
mea <- input$mean3

t <- seq(from=0,to=xmax,by=1)
d <- dpois(t,mea)
plot(t,d,type="h")
points(t,d,pch=20)
abline(v=mea,col=2)
})
```

# ポアソン分布に従う「種」とそこからの指数関数的増加

[No Control Genes Required: Bayesian Analysis of qRT-PCR Data」(http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0071448)

種：標本中のmRNA分子の個数
qRT-PCR : PCRを用いた指数関数的増幅による検出

「ゼロ」に対応する値

「ノイズ」

「観測シグナルは指数関数的増幅数」

「指数関数的増幅数はmRNA分子個数に対数正規乱数が掛け合わさった値」→に基づくポアソン分布」

# ポアソン対数正規分布にできてただの対数正規分布にできないこと

「離散的」であること

「０」であること

# Poisson Log-normal distribution ポアソン対数正規分布

平均(mean)と中央値(mode)とを指定して、ポアソン対数正規分布を見てみよう。

表示範囲はxmaxで調整する。

```{r, echo=FALSE}
inputPanel(
  sliderInput("mod4", label = "mode:",
              min = 0, max = 2000, value = 500, step = 1)
,
sliderInput("mean4", label = "mean:",
              min = 0, max = 10000, value = 2000, step = 1)
,
sliderInput("xmax4", label = "xmax:",
              min = 0, max = 10000, value = 3000, step = 1)
)

renderPlot({
xmax <- input$xmax4
mea <- input$mean4
mod <- input$mod4
med <- (mea^2*mod)^(1/3)


t <- seq(from=0,to=xmax,by=1)
d <- my.dpoislnorm (t,mea,mod)
plot(t,d,type="h")
abline(v=mea,col=2)
abline(v=mod,col=3)

})
```


## ノイズか真のシグナルか

シグナルの中央値(SignalMode)とシグナルの平均(SignalMean)を指定し、さらに、ノイズの中央値と平均(NoiseMode,NoiseMean)を指定して、「シグナルがノイズ」から際立っているかどうかを視覚的に見てみる。

シグナルは赤。
ノイズは黒。


このモデルで考えるのであれば、「ある特定のシグナル強度」が得られたときに、それが「真のシグナルなのか、ノイズなのか」の尤度が計算できる。

```{r, echo=FALSE}
nS = 10
nN = 10^3
ord = rep(1,nS+nN)
sid = sample(1:(nS+nN),nS)
nid = (1:(nS+nN))[-sid]

inputPanel(
  sliderInput("mod5", label = "SignalMode:",
              min = 0, max = 2000, value = 500, step = 1)
,
sliderInput("mean5", label = "SignalMean:",
              min = 0, max = 10000, value = 2000, step = 1)
,
  sliderInput("mod6", label = "NoiseMode:",
              min = 0, max = 100, value = 5, step = 0.1)
,
sliderInput("mean6", label = "NoiseMean:",
              min = 0, max = 200, value = 15, step = 0.1)
,
sliderInput("ymax5", label = "ymax:",
              min = 0, max = 10000, value = 3000, step = 1)
)

renderPlot({
ymax <- input$ymax5
meaS <- input$mean5
modS <- input$mod5
meaN <- input$mean6
modN <- input$mod6

rS <- my.rpoislnorm (nS,meaS,modS)
rN <- my.rpoislnorm (nN,meaN,modN)
col <- rep(1,nS+nN)
col[sid] <- 2
val <- rep(0,nS+nN)
val[sid] <- rS
val[nid] <- rN
plot(val,type="h",col=col,ylim=c(-ymax*0.1,ymax))

})
```