2014-02-19

配列し損ねるとき

論文を紹介していただく機会があった。こちらの論文。
RNA配列上に遺伝子多型があって、そのアレル別の発現量比推定をしたい、というときに、RNA-seqのデータをレファレンスゲノムにマッピングする過程をはさむと、ゲノムレファレンスアレルを持つ配列がうまく張り付く確率は、非レファレンスアレルを持つ配列が張り付く確率よりも高い、というバイアスが入ってしまうね、それを避けるために、マッピングのときにゲノム配列を２つ用意して、片方はいわゆるレファレンスアレルを持たせ、もう片方は別のアレルを持たせる、というやり方や、多型箇所をNなどの「不確定塩基文字」にした配列にマッピングさせるというやり方があるけど、そうするとどうなるの？という論文でした。

こちらの図にあるように、「普通のレファレンス」にマップする場合を水平軸に、そうでなくて「比レファレンスアレルの配列もマッピングされやすくするようなレファレンス」にマップする場合を鉛直軸にとって、アレル別の発現量比の対数をプロットすると、おおまかには、よい正相関をとりますが、「普通のレファレンス」を使ったときには、レファレンスアレルの割合が大きめになることが、『プロット雲がy=xのラインより下・右に来る』ことで示されている。
そのほかに見て取れることは、相関の雲の下・右の方向にはばらつきが大きめながら、上・左の方向は『硬い』傾向を示しています。

このプロットの形の説明として、「非レファレンスアレルを持つ配列」が、単にマップされ損ねることがある、ということで説明がつくかな、と議論になったので確かめてみます。
この簡単なモデルで大まかにはプロットの様子が再現できるようです。
簡単なモデルでのあてはまりがよいということは、モデルを構成しているパラメタ推定をして適当な補正をすることもそれほど悪いアイディアではない、ということですね。

# 遺伝子数
n <- 10000
# ２アレルの量比のログを乱数発生
logr <- rnorm(n)
# 比そのもの
r <- exp(logr)

# 各遺伝子でどれくらいのデプスで読まれるかを、平均デプスkのポアソン分布とする（ポアソン分布をここで使うのは、この分野でよくやること)
k <- 50
# 非レファレンスアレルの配列がマップされる確率を１より小さくして、「マップしそこね」を想定する
P <- 0.8
# 各遺伝子が２アレル合わせて何デプスになるかを乱数発生
q <- rpois(n,k)
# 読まれたリードの内訳をレファレンスと非レファレンスとに分ける
q.ref <- rep(0,n)
q.alt <- rep(0,n)
# そのうち、レファレンスゲノムにマップしそこねやすい場合にきちんとマップされた本数を別途乱数で与える
q.alt.2 <- q.alt
# ループしなくてもできそうだけれど、R初心者にはループの方がよいか…
for(i in 1:n){
# 発現量比から、発現割合を計算
	tmp.prob <- 1/(1+r[i])
# 発現量比に応じて、デプス本数を２アレル別の本数に分ける
	q.ref[i] <- sum(sample(0:1,q[i],replace=TRUE,prob=c(tmp.prob,1-tmp.prob)))
	q.alt[i] <- q[i]-q.ref[i]
# 読まれた非レファレンスアレル配列を確率Pでマップ成功させ、その本数を数える
	q.alt.2[i] <- sum(sample(0:1,q.alt[i],replace=TRUE,prob=c(1-P,P)))
}
# プロット用に発現量比の対数を計算
log.r.obs <- log(q.alt/q.ref)
log.r.obs.2 <- log(q.alt.2/q.ref)
# コプロット
plot(log.r.obs,log.r.obs.2)
# 補助線
abline(v=0)
abline(h=0)


#hist(log.r.obs)