数学・コンピュータ関連の姉妹ブログ『ryamadaのコンピュータ・数学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

RFUシグナル強度と 正規・対数正規・ポアソン対数正規・二項分布 〜ノイズとシグナル、ホモとヘテロ〜

RFU PCR DNA鑑定 ショットノイズ 対数正規分布 ポアソン対数正規分布 ヘテロインバランス

RFUシグナル強度と正規・対数正規・ポアソン対数正規・二項分布: ノイズとシグナル、ホモとヘテロ

RFUシグナル強度と正規・対数正規・ポアソン対数正規・二項分布: ノイズとシグナル、ホモとヘテロ

  • スライドの図表作成のためのRmd。スライドはデータ作成方法についての説明が少ないので、データ作成については、以下を。
---
title: "RFUのばらつき"
author: "ryamada"
date: "Thursday, May 07, 2015"
output: html_document
---

# はじめに
この文書は日本法医学会(2015年6月於高知)に合わせて開かれた、京都大学医学部法医学教室主催の方数学勉強会６月特別回での話題提供『RFUシグナル強度と^K正規・対数正規・ポアソン対数正規・二項分布^K〜ノイズとシグナル、ホモとヘテロ〜』の補助資料です。


DNA鑑定に用いるSTRのPCR増幅シグナルであるRFUのバラツキを、そのノイズ・シグナル識別と、ヘテロ接合体の２アレルのRFUの違いについての考察をするために、いくつかの仮定をしてシミュレーショナルなデータ作成をRで行って図を作成しました。

そのコード記録です。

学術論文として厳密性を目指していませんし、図が描ければOK、というスタンスで書いた走り書きですので、乱数試行をすると不具合のある図が描かれたりするかもしれませんが、勉強会の内容を復習するためなどの役には立つと思いますので、文書公開することにしました。
以上についてご容赦いただいたうえでお使いください。

なお、話題提供時のパワーポイントスライドは、http://www.genome.med.kyoto-u.ac.jp/wiki_tokyo/index.php/法数学勉強会RFUシグナル強度と正規・対数正規・ポアソン対数正規・二項分布〜ノイズとシグナル、ホモとヘテロ〜　からアクセスできます。


# PCRとRFU (relative fluorescence unit)

サンプル中にあるDNA(もしくはRNA)分子をテンプレートとしてPCR増幅をすると、その初期分子コピー数の多寡に応じて、あるサイクル数から、増幅分子数を表すRFUシグナルの測定値が立ち上がり、やがてプラトーに達する。

単純化すると
$$
RFU = \frac{Max}{1+e^{-k(Cycle-g)}}
$$
という式で表せるようなシグモイドカーブになる。

Maxはプラトー値であり、gはテンプレート分子数によって決まる値である。
実験から、gの値は、テンプレート分子数の対数の一次線形式で表されることが知られている。
kはカーブの傾斜の程度を決めるパラメタである。

# サンプルの希釈系列

サンプルの倍々希釈系列でサイクルごとのRFUを描くと、等間隔のシグモイドカーブが並ぶ。

以下のコードは、コード内のコメントに書いたような条件で、カーブを描く。

```{r}
# x はサイクル数
# pは濃度
# 濃度の対数に比例して、シグモイドカーブはサイクル数軸に平行移動する
# 濃度がp0のときに、サイクル数x0でRFU値の中間値を取るような設定とする
# aは傾きを決めるパラメタ。aが大きいとシグモイドカーブの傾きが大きくなる
my.sig.curve.2 <- function(x,p,m=0,M=2000,x0=27,p0=2^(-2),a=1){
  b <- x0 + log10(p0)*a
	gp <- -a*log10(p)+b
	m + M/(1+exp(-a*(x-gp)))
}
# DNA鑑定のサイクル数は27くらいらしい
x <- seq(from=0,to=40,length=121)
# x <- 0:40 # No. cycle
# 濃度は実濃度ではなく、基準濃度の何倍か、のような扱いを考慮
# 倍々希釈系列を７段階作る
ps <- 2^(0:(-6)) # conc
a <- 2.2
min.rfu <- 0
max.rfu = 1800
outs <- matrix(0,length(x),length(ps))
for(i in 1:length(ps)){
  outs[,i] <- my.sig.curve.2(x,ps[i],p0 = 2^(-2.5),a=a,M=max.rfu)
}

# 濃度別のカーブ。サイクル数とRFUの関係
matplot(x,outs,type="l",xlab="No. cycles",ylab="RFU",xlim=c(20,35))
abline(v=27)
```

# サンプル濃度とシグモイドカーブ

サンプル濃度を希釈系列ではなく、実濃度を等間隔で変化させて同様にシグモイドカーブを描く。
```{r}
ps <- seq(from=0,to=1,length=100)
ps <- ps[-1]
a <- 2.2
min.rfu <- 0
max.rfu = 1800
outs <- matrix(0,length(x),length(ps))
for(i in 1:length(ps)){
	outs[,i] <- my.sig.curve.2(x,ps[i],p0 = 2^(-2.5),a=a,M=max.rfu)
}

# 濃度別のカーブ。サイクル数とRFUの関係
matplot(x,outs,type="l",xlab="No. cycles",ylab="RFU")
abline(v=27)
```

# 特定のサイクル数のときのテンプレート濃度とRFUの関係

DNA鑑定では、特定のサイクルにてRFU値を測定する。

横軸に試料濃度を、通常の軸または対数軸で取って、RFUをプロットする。

横軸を対数軸にした場合には、対数濃度とRFUがシグモイドカーブを描くことがわかる。


```{r}
# サイクル数が27のときの、濃度ごとのRFU
plot(ps,outs[which(x==27),],type="b",pch=20,cex=0.5,xlab="Relative Concentration",ylab="RFU",main="27 cycles")

# conc = 1,1/2,1/4,1/8,1/16
abline(v=2^(-6:0))

# 横軸を対数とした片対数グラフ
plot(log2(ps),outs[which(x==27),],type="b",pch=20,cex=0.5,xlab="log2(conc)",ylab="RFU",main="27 cycles")
abline(v=log2(2^(-6:0)))
```

# ホモ型とヘテロ型とのRFUの比

ヘテロ接合体の場合、あるアレルの試料濃度は、ホモ接合体の半分になる(すべての調整条件が等しい場合)。

したがって、ヘテロ接合体の場合のRFU値はホモ接合体より低くなる。
しかし、そのRFUピーク値は、半分になるわけではない。それは、濃度とRFUとの関係が比例関係ではなくシグモイダルだからである。

それを具体的に計算する。

上の例でもシミュレーションした７段階の倍々希釈系列における、濃度の倍変化をホモ/ヘテロ関係と見立てて、６対の架空の「ホモ・ヘテロ」ペアについて、そのRFU値の比を算出する。

RFU値が相対的に高く、プラトー部分に近いときは、ホモ・ヘテロ比は１に近く、逆に立ち上がり部分に近いときは、ヘテロはホモの半分よりも小さくなることがわかる。

```{r}
ps <- 2^(0:(-6)) # conc
a <- 2.2
min.rfu <- 0
max.rfu = 1800
outs <- matrix(0,length(x),length(ps))
for(i in 1:length(ps)){
  outs[,i] <- my.sig.curve.2(x,ps[i],p0 = 2^(-2.5),a=a,M=max.rfu)
}

# 濃度別のカーブ。サイクル数とRFUの関係
#matplot(x,outs,type="l",xlab="No. cycles",ylab="RFU",xlim=c(20,35))
#abline(v=27)
RFU.27 <- outs[which(x==27),]
RFU.27[2:length(RFU.27)]/RFU.27[1:(length(RFU.27)-1)]
```


# シグモイドカーブの傾きがばらつく

実際に同じテンプレート量で増幅しながらRFUを測定すると、カーブの立ち上がり部分はそれほどぶれないが、カーブの傾きはばらつく。
そのバラツキの具合は中心となるある角度から角度が増える方、角度が減る方について対称ではなく、角度が減る方についてバラツキが大きくなる傾向がある。

これをシミュレーションするために、特定のシグモイドカーブ傾斜のパラメタ値に対して、指数分布に従う乱数を与えると、傾きのばらつきが実観測データに似せることができる。

以下では、そのようにasに値を与えている。

この分子実験的説明は、増幅効率が悪い場合があって、その場合には、その遅れが次のサイクルにも影響を引き継ぎ、サイクル数の進行とともに拡大する、というようにするのがよいのかもしれないが、その正確なところは不明である。

ここでは、あくまでも実験データにある程度合致するカーブを描くための方便としてasを用いる。



```{r}
ps <- 2^(-2) # conc
a <- 2.2
as <- a - exp(rnorm(50,0,0.2))
as <- a - qexp(seq(from=0,to=1,length=10),2)
#as <- a - rnorm(50,0,0.3)
as <- sort(as,decreasing=TRUE)
min.rfu <- 0
max.rfu = 1800
outs <- matrix(0,length(x),length(as))
for(i in 1:length(as)){
  outs[,i] <- my.sig.curve.2(x,ps,p0 = 2^(-20),a=as[i],M=max.rfu)
}

# 濃度別のカーブ。サイクル数とRFUの関係
matplot(x,outs,type="l",xlab="No. cycles",ylab="RFU",xlim=c(10,35))

```


# 指数分布的な傾き遅れがあるときのRFU分布

asに指数分布乱数を与えることで、実験データを模すことに成功した。

その方法を用いて、実際にある特定のサイクル数におけるRFUがどのような分布をするかをシミュレーションする。

上の例よりは、指数乱数を小さくしないと、RFU分布が大きすぎるので、その点を調整して分布を取る。

そうすると、RFUが小さい方に裾を引く分布となることがわかる。

この条件でのシグモイドカーブと、あるサイクル数でのRFU分布とをそれぞれ示す。

```{r}
ps <- 2^(-2) # conc
a <- 2.2
#as <- a - exp(rnorm(50,0,0.2))
#as <- a - qexp(seq(from=0,to=1,length=10),2)
as <- a - rexp(1000,50)
as[which(as<0)] <- 0
#as <- a - rnorm(50,0,0.3)
as <- sort(as,decreasing=TRUE)
min.rfu <- 0
max.rfu = 1800
outs <- matrix(0,length(x),length(as))
for(i in 1:length(as)){
  outs[,i] <- my.sig.curve.2(x,ps,p0 = 2^(-20),a=as[i],M=max.rfu)
}

# 濃度別のカーブ。サイクル数とRFUの関係
matplot(x,outs,type="l",xlab="No. cycles",ylab="RFU",xlim=c(10,35))
```

```{r}
hist(outs[40,],xlab="RFU",main="")
```

# 試料のアレル・インバランス

通常、ゲノムDNAのヘテロ型サイトのDNA分子数は細胞回収の段階では１：１で回収されると考えられる。
それは各細胞が１：１で保有しているからである。

それをDNA保存液とする段階、保存液から実験用に取り出す過程で２種類のアレルに相当するDNA分子数にインバランスが生じうる。

以下のシミュレーションでは、保存液を作成するまでは、１：１が維持され、そこから実験用に取り出すときにのみ、アレル・インバランスが生じるという単純な設定を採用した。

実験が好条件で行なわれる場合というのは、25-100ngのヒトゲノムを用いる。
それは10,000強の遺伝子コピー(DNA分子コピー)に相当するとされる。

それは、大まかに言えば、1mlの血液に5000個の白血球があるとして、その細胞がDNAの供給源であると考えると、全血0.1ml分に相当する500,000個の細胞から、DNA分子として
1,000,000コピー分が抽出・保存されたとして、その保存液の1/100を実験用に取り出した、と言うような設定に相当する。

10,000DNA分子に占める片方のアレルの割合がどのようにばらつくかをヒストグラムにする。
次いで、1000DNA分子、100DNA分子が実験用に取り出されたとした場合も同様にヒストグラム化する。

分子数が少ないほどバラツキは大きいが、２アレルの比が1:2(33% : 66%)となるような差は、かなり稀であることがわかる。

```{r}
n <- 5*10^5
g <- rep(0:1,n)

n.iter <- 10000
N <- 10000
n.a <- rep(0,n.iter)
for(i in 1:n.iter){
  n.a[i] <- sum(sample(g,N))
}
par(mfcol=c(1,3))
hist(n.a/N,xlab="Fraction of 1 allele",main="N=10000")

n <- 5*10^5
g <- rep(0:1,n)

n.iter <- 10000
N <- 1000
n.a <- rep(0,n.iter)
for(i in 1:n.iter){
  n.a[i] <- sum(sample(g,N))
}
hist(n.a/N,xlab="Fraction of 1 allele",main="N=1000")

n <- 5*10^5
g <- rep(0:1,n)

n.iter <- 10000
N <- 100
n.a <- rep(0,n.iter)
for(i in 1:n.iter){
  n.a[i] <- sum(sample(g,N))
}
hist(n.a/N,xlab="Fraction of 1 allele",main="N=100")
abline(v=c(1/3,2/3),col=2)
par(mfcol=c(1,1))
```



初期回収量を減らすしてもあまり変わらない。

```{r}
n <- 5*10^4
g <- rep(0:1,n)

n.iter <- 10000
N <- 10000
n.a <- rep(0,n.iter)
for(i in 1:n.iter){
  n.a[i] <- sum(sample(g,N))
}
par(mfcol=c(1,3))
hist(n.a/N,xlab="Fraction of 1 allele",main="N=10000")

n <- 5*10^4
g <- rep(0:1,n)

n.iter <- 10000
N <- 1000
n.a <- rep(0,n.iter)
for(i in 1:n.iter){
  n.a[i] <- sum(sample(g,N))
}
hist(n.a/N,xlab="Fraction of 1 allele",main="N=1000")

n <- 5*10^4
g <- rep(0:1,n)

n.iter <- 10000
N <- 100
n.a <- rep(0,n.iter)
for(i in 1:n.iter){
  n.a[i] <- sum(sample(g,N))
}
hist(n.a/N,xlab="Fraction of 1 allele",main="N=100")
abline(v=c(1/3,2/3),col=2)
par(mfcol=c(1,1))
```


# まとめ

ノイズは、ポアソン対数正規分布を仮定した。

真のシグナルには、シグモイドカーブを仮定し、そこに傾き係数のバラツキとして指数乱数を仮定するとともに、サンプルのアレル・インバランスに二項分布を仮定した。


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


以下では、ノイズにポアソン対数正規分布を使い、
真のシグナルには、シグモイドカーブと傾き係数の指数乱雑項を用い、ヘテロ接合体のアレルインバランスは二項分布の代わりに正規分布を用いた。

ホモ接合体のRFUがいわゆるカットオフ値150よりは大きいがそれほど大きくはないような条件での、ヘテロ接合体の２アレルRFU値がどれくらい低く、違いがあるか、そして低い方のアレルがノイズに埋もれるのかをシミュレーションしている。


黒がノイズで、赤はホモ接合体のRFU。ホモ接合体のRFU(赤)がばらつくのは、傾斜係数のバラツキを反映している。

そして、赤のシグナルの右隣りに２本の緑のシグナルをヘテロ接合体のそれとして表示した。
この緑のシグナルは、左隣りの赤のシグナルの場合と傾斜条件をそろえている(実際には、ホモとヘテロが相並ぶことはないので、あくまでも人工的な表示設定である)。
そのうえで、二つのアレルのサンプル濃度は正規分布を用いて変えている。(片方のアレルが多ければもう片方のアレルが少なくなるという関係)。


```{r}
# 対数正規分布を扱うための諸関数と
# ポアソン対数正規分布を扱うための諸関数
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)
}

# ノイズシグナル数
nN <- 1100
# ホモシグナル数
nHomo <- 10
# ヘテロ箇所数はホモシグナル数と同じ
nHetero <- nHomo
#nS = 10
#nN = 10^3
# ノイズとホモとヘテロのすべてのシグナルを横一列に並べるために
# 番号をつける
ord = rep(1,nHomo+nHetero*2+nN)
# ホモシグナルは100番地おきとする
sid = (1:nHomo)*100
# ヘテロの２つのシグナルは、２つのホモのシグナルの中間に
# 少しの間隔をあけて確保する
heteroid.1 <- sid + 40
heteroid.2 <- sid + 60
# ノイズシグナルであるべき番地
nid = (1:(nN))[-c(sid,heteroid.1,heteroid.2)]

# 描図設定
ymax <- 3000
#meaS <- 2000
#modS <- 500

# ノイズの分布を決めるパラメタ。平均値と中央値
meaN <- 8
modN <- 7

# 真のシグナルのためのシグモイドカーブパラメタ設定
ps <- 2^(-2) # conc
a <- 2.2
# 傾斜のバラツキ設定(指数乱数)
#as <- a - exp(rnorm(50,0,0.2))
#as <- a - qexp(seq(from=0,to=1,length=10),2)
as <- a - rexp(1000,50)
as[which(as<0)] <- 0
#as <- a - rnorm(50,0,0.3)
as <- sort(as,decreasing=TRUE)
min.rfu <- 0
max.rfu = 1800
# outsはホモ用、outs.hetero.1,2はヘテロ用
outs <- matrix(0,length(x),length(as))
outs.hetero.1 <- matrix(0,length(x),length(as))
outs.hetero.2 <- matrix(0,length(x),length(as))
ps.hetero <- ps/2
for(i in 1:length(as)){
  outs[,i] <- my.sig.curve.2(x,ps,p0 = 2^(-20),a=as[i],M=max.rfu)
  # 正規分布を使ってアレル別の濃度を与えている
  # 正規分布のsdを調整するパラメタ
  sd.param <- 4
  tmp.p.1 <- rnorm(1,ps.hetero,ps.hetero/sd.param)
  tmp.p.2 <- ps.hetero*2-tmp.p.1
  outs.hetero.1[,i] <- my.sig.curve.2(x,tmp.p.1,p0 = 2^(-20),a=as[i],M=max.rfu)
  outs.hetero.2[,i] <- my.sig.curve.2(x,tmp.p.2,p0 = 2^(-20),a=as[i],M=max.rfu)
}

# 濃度別のカーブ。サイクル数とRFUの関係
# matplot(x,outs,type="l",xlab="No. cycles",ylab="RFU",xlim=c(10,35))
#rS <- my.rpoislnorm (nS,meaS,modS)

# ホモ・ヘテロのRFUを上で作成した中から、特定のサイクル数相当として取り出す
# サイクル数:これを変えるとRFUを高くしたり低くしたりできる
cycle <- 44
rHomo <- sample(outs[cycle,],nHomo)
# ヘテロの場合は、ペアで値を生成しているので、そのようにする
heteroid <- sample(1:length(as),nHetero)
rHetero <- c(rbind(outs.hetero.1[cycle,heteroid],outs.hetero.2[cycle,heteroid]))
# ノイズのRFUを乱数発生(ポアソン対数正規分布)
rN <- my.rpoislnorm (nN,meaN,modN)

# ノイズは黒、ホモは赤、ヘテロは緑、と色指定する
col <- rep(1,nHomo+nHetero*2+nN)
col[sid] <- 2
col[c(heteroid.1,heteroid.2)] <- 3
val <- rep(0,nHomo+nN)
val[sid] <- rHomo
val[c(heteroid.1,heteroid.2)] <- rHetero
val[nid] <- rN[nid]
ymax <- max(val)
plot(val,type="h",col=col,ylim=c(-ymax*0.1,ymax))
```

サンプルの
ヘテロ・インバランスの具合は、片方のアレルの濃度が以下のような分布になるものとして行っている。
```{r}
sd.param <- 4
hist(rnorm(10^4,ps.hetero,ps.hetero/sd.param))
```

logo

図書室

日記全体の目次
(記事一覧)

疾患遺伝子解析を学ぶためのカタログ(H18 3月改訂)の目次

 教科書

京都大学大学院医学研究科集中講義 全８限Part I(2007)
京都大学大学院医学研究科集中講義 全８限Part II(2008)
京都大学大学院医学研究科統計遺伝学講義 全８限(2009)

以下の関連サイトへのリンク道案内ページ
遺伝子多型解析Wiki
東大医科研HGCゲノム機能解析分野(旧)
理研リウマチチーム解析資料ページ


疾患遺伝子解析を学ぶためのカタログ(H18 3月改訂)
遺伝統計学のための統計学基礎 全９限
Haploviewを用いた連鎖不平衡マッピング解析実習 全９限
Arlequinを用いた集団遺伝学実習 全７限
一から始める遺伝子多型・形質関連解析
多数のGWAS/GWASメタアナリシスとその周辺解析用のアプリケーションへのリンク

複合遺伝性疾患とは

----------News----------
20万アクセス
2008/07/20
40万アクセス
2010/09/12


-------------------------