Md ハプロタイプ頻度の最尤推定

提供: StatGenetKyotoU
移動: 案内検索
ハプロタイプ頻度の最尤推定
カテゴリ 名称 説明
医学生物学トピック SNP ハプロタイプ
数学トピック 最尤法 もっともあてはまりの良い値の近似値を出す
R関数・Rスキル

極値と微分=0

  • f(x)=x^3-2x^2+x-1を微分するとf'(x) = 3x^2-4x+1
    • 黒の曲線がf(x)、赤の破線がf'(x)。緑のラインはf'(x)=0の根がx=1/3,1であることを示し、水色のラインはx=1/3,1のときにf(x)が極大値・極小値を取ることを示す

微分.jpeg

x <- seq(from=-0.5,to=1.5,length=100)
f<-x^3-2*x^2+x^1-1
fp <- 3*x^2-4*x+1
matplot(x,cbind(f,fp),type="l")
abline(h=0,col=3)
abline(v=c(1/3,1),col=3)
abline(h=c(-23/27,-1),col=5)

アレル頻度の最尤推定

  • Aがn本Bがm本観測された。Aの割合をp、Bの割合を1-pのとき、この観察をする確率(観察n,mにおけるpの尤度)は
    • f(p)=(n+m)!/(n!m!) p^n(1-p)^m
  • g(p) = log(f(p))を取ると、n乗・m乗の部分が楽になるのでそのようにすると
    • g(p) = log((n+m)!/(n!m!))+ n log(p) + m log(1-p)
  • どちらも n/(n+m) = 0.25 のところで最大値になっていることがわかる

尤度対数尤度.jpeg

p <- seq(from=0,to=1,length=100)
n <- 10; m <- 30;
fp <- factorial(n+m)/(factorial(n)*factorial(m))*p^n*(1-p)^(m)
par(mfcol=c(1,2))
plot(p,fp,type="l")
abline(v=n/(n+m),col=2)
gp <- log(fp)
plot(p,gp,type="l")
abline(v=n/(n+m),col=2)
  • f(p),g(p)を微分してみる
    • f'(p) = (n+m)!/(n!m!)(np^(n-1)(1-p)^m - mp^n(1-p)^(m-1)) = (n+m)!/(n!m!)p^(n-1)(1-p)^(m-1)(n(1-p) - mp)
    • g'(p) = n/p-m/(1-p) = (n(1-p)-mp)/(p(1-p))
  • と、どちらもn(1-p)-mp=0がf'(p)=0,g'(p)=0から出てくることがわかる
    • これはp=n/(n+m)のことで、これがグラフの極大値を与えるpの値で、これがpの最尤推定値

2SNPのハプロタイプの最尤推定

4ハプロタイプの頻度の自由度は3

  • 2SNPの作るハプロタイプは4種類、その頻度をh1,h2,h3,h4とすれば、h1+h2+h3+h4=1なので、自由度は3
  • 連鎖平衡状態では、4ハプロタイプ頻度は、2SNPのアレル頻度p,qのみから以下のように定まるから、自由度は(変数の数であって)2
    • h1=pq,h2=p(1-q),h3=(1-p)q,h4=(1-p)(1-q)
  • 連鎖不平衡状態は、連鎖平衡状態から、ずれているので
    • h1=pq+d1, h2=p(1-q)+d2, h3=(1-p)q +d3, h4=(1-p)(1-q)+d4とおける
    • 新たに4変数 d1,d2,d3,d4が導入されたが、h1=h2=p,h1+h3=qという制約があるから、d1+d2=0,d1+d3=0という制約があることになり
    • 結局、d1=d4,d2=d3=-d1となって
    • h1=pq+d1,h2=p(1-q)-d1,h3=(1-p)q-d1,h4=(1-p)(1-q)+d1のように、p,q,d1の3変数で表される…自由度は3だから3変数で表されるのが道理
  • こんな表し方にすると
    • h1 = pq + r(pq(1-p)(1-q))^0.5
    • h2 = p(1-q) - r(pq(1-p)(1-q))^0.5
    • h3 = p(1-q) - r(pq(1-p)(1-q))^0.5
    • h4 = p(1-q) + r(pq(1-p)(1-q))^0.5
  • p=qのときr=1で以下の通り、"absolute LD"
    • h1 = p
    • h2 = 0
    • h3 = 0
    • h4 = (1-p)

4ハプロタイプの尤度関数

  • 2SNPのディプロタイプは3x3=9種類
  • 3x3行列状にg11,g12,...,g33とする
  • 対数尤度関数は
    • g11 log(h1^2) + g12 log(2h1h2) + g13 log(h2^2) + g21 log(2h1h3) + h22 log(2h1h4 +2h2h3) + g23 log(2h2h4) + g31 log(h3^2) + g32 log(2h3h4) + g33 log(h4^2)
  • 式を整理して
    • (2g11+g12+g21)log(h1)+(2g13+g12+g23)log(h2)+(2g31+g21+g32)log(h3)+(2g33+g23+g32)log(h4) + g22 log(2(h1h4+h2h3))
  • h1,h2,h3,h4は4変数、自由度3なので、少々面倒くさいので次のようにしてみる
  • H=(h1,h2,h3,h4)の値セットを適当に作って、その尤度を計算して、それが最大になる値を調べる
    • Hはなるべくいろいろな値を取らせたい(ただしh1+h3+h3+h4=1)
    • このような乱数セットはディリクレ分布からの乱数として取り出せる(ディリクレ分布はこちら)
  • まず、9ジェノタイプカウントとハプロタイプ頻度から対数尤度を計算する関数を作る
my.loglike <- function(h,g){
	(2*g[1,1]+g[1,2]+g[2,1])*log(h[1])+(2*g[1,3]+g[1,2]+g[2,3])*log(h[2])+(2*g[3,1]+g[2,1]+g[3,2])*log(h[3])+(2*g[3,3]+g[2,3]+g[3,2])*log(h[4]) + g[2,2]*log(2*(h[1]*h[4]+h[2]*h[3]))
}

  • ディリクレ分布からの乱数を発生させる
n.sets <- 1000
install.packages("MCMCpack")
library(MCMCpack)
R <- rdirichlet(n.sets,rep(1,4))
# ばらつき具合を確認する
plot(as.data.frame(R))
# 3次元なら…
R.3 <- rdirichlet(n.sets,rep(1,3))
library(rgl)
plot3d(R.3)

ディリクレ3次元.png

  • 4ハプロタイプ頻度の行列Rにmy.loglike()関数を作用させて、対数尤度が大きいハプロタイプにはどんなものがあるかを見てみます
g <- matrix(c(25,50,25,50,100,50,25,50,25),3,3,byrow=TRUE)
# apply関数はRの行(第2引数:1)ごとに関数my.loglike()を作用させています、ただし、my.loglike()関数の第2引数gには別途、ディプロタイプ行列gを与えています
out <- apply(R,1,my.log.like,g=g)
# 対数尤度が出たので、その値の大きい順を確認して
ord <- order(out,decreasing=TRUE)
# 対数尤度が大きい方から、4ハプロタイプの頻度を見てみます
R[ord[1:10],]
# 対数尤度も見てみます
out[ord[1:10]]

その他関連