座標で表す - ryamadaの遺伝学・遺伝統計学メモ

ハプロイド・ディプロイド
- あるマーカーのアレル数がkのとき、アレル頻度に関する情報はk個の値で(も)与えられる。ただし、総和が１という制約があるので、自由度はk-1
- この自由度k-1の情報を幾何的にどう配置するか、はいくつかやり方がある
  - 例えば、アレルに勝手に順序をつけて、１列に並べてしまう
  - k次元空間に配置する
  - k-1次元空間に配置する
- このハイチの仕方は、ひとまず、棚上げする(k-1次元空間に置くのがおそらく良い)
- ディプロイドにするとは、ハプロイドの組み合わせ〜積のこと
  - 幾何的にも積にする
  - 独立な積は幾何的には直積
  - HWEは独立
- ただしディプロイドの観測は、アレルの組み合わせとしてしか観測できないので、情報が減じた状態
  - したがって、ハプロイドの直積を減次元状態に射影する必要がある
以上のようなことをもっとも簡単な1SNPのHWE状態の観測ディプロタイプに対してやってみる
そのルールが一般化した形で作れれば…
2SNPが4ハプロタイプを作る
4ハプロタイプのディプロイドのタイプは、父方・母方由来を区別すれば $4^2=16$ 通り
父母由来を区別しないで、 $(AB,ab),(Ab,aB)$ を区別すれば10通り
$(AB,ab),(Ab,aB)$ を区別しなければ、9通り
4ハプロタイプの頻度を与えよう

p1 <- c(0.3,0.7)
p2 <- c(0.4,0.6)
f <- 0.6 # r^2の値
# 連鎖平衡の下でのハプロタイプ頻度
h.LE <- p1 %*% t(p2)
# 連鎖不平衡の下でのハプロタイプ頻度
delta <- f^(1/2) * prod(h.LE)^(1/4)
h.LD <- h.LE + delta * matrix(c(1,-1,-1,1),2,2)
# r^2は以下で求まる統計量
chisq.test(h.LE,correct=FALSE)
chisq.test(h.LD,correct=FALSE)

16種類のディプロタイプ頻度はハプロタイプ頻度の掛け合わせ。ただし、後述するように、３ｘ３ディプロタイプの扱いを「回転対称」で扱いたいので、ハプロタイプも「回転対称」な順序に並べることとする

h.LD.2 <- c(h.LD[1],h.LD[2],h.LD[4],h.LD[3])
h.LE.2 <- c(h.LE[1],h.LE[2],h.LE[4],h.LE[3])
diplotype.16 <- c(h.LD.2) %*% t(c(h.LD.2))

16タイプの頻度が4x4行列に入っている
これを「３ｘ３行列」の番地に対応付けたい
- 4つのハプロタイプを(1,0),(0,1),(-1,0),(0,-1)と2次元平面の４つの単位ベクトルに対応付けると、ディプロタイプは対応するハプロタイプのベクトルの和で対応づけられる。
- ホモ接合体に対応する２点が、正方形の頂点に対応するが、この正方形の辺がｘ軸、ｙ軸と並行になるように回転して、「番地付け」がうまく行く様に平行移動することにする

mm <- matrix(c(1,0,0,1,-1,0,0,-1),byrow=TRUE,ncol=2)
x<-c(t(outer(mm[,1],mm[,2],"+")))
x<-c(t(outer(mm[,1],mm[,1],"+")))
y<-c(t(outer(mm[,2],mm[,2],"+")))
xy <- cbind(x,y)
rot <- matrix(c(cos(5*pi/4),-sin(5*pi/4),sin(5*pi/4),cos(5*pi/4)),byrow=TRUE,2,2)

xy.rot <- t(rot %*% t(xy))

xy.rot.2 <- xy.rot - min(xy.rot)
xy.rot.3 <- xy.rot.2/max(xy.rot.2)*2 +1

xy.rot.3

xy.rot.4 <- zapsmall(xy.rot.3)
xy.rot.4

ハプロタイプペアとディプロタイプ番地との対応を取って、ディプロタイプ頻度を計算する

diplotype.freq <- matrix(0,3,3)

for(i in 1:3){
	for(j in 1:3){
		tmp <- which(xy.rot.4[,1] == i & xy.rot.4[,2] ==j)
		diplotype.freq[i,j] <- sum(diplotype.16[tmp])
	}
}

diplotype.freq

ちなみに、HWEが成り立っているときは、ディプロタイプの３ｘ３表のカイ自乗値は、ハプロタイプの２ｘ２表のそれのカイ自乗値fに対して $f+f^2$ になるらしい

dip.test.out <- chisq.test(diplotype.freq,correct=FALSE)
dip.test.out$statistic- f*(1+f)