図示する

• 久しぶりにSNPデータのR処理の話し
• SNPデータからIBD推定としてz0,z1,z2を算出する(こちら)が、この３つの値は足し合わせると１になるので、自由度が２〜２次元平面に描ける。
• z0,z1,z2をそれぞれ正三角形の頂点にとって描く

# サンプル数
Nsample<-100
# rdiriclhet()はz0,z1,z2からなる３つ組値データ(z0+z1+z2=1, z0,z1,z2>=0)をランダムに発生させる関数
library(MCMCpack) # rdirichlet()を持つパッケージ

# 引数c(10,3,1)はz0の値が重くなりがちになるような重みづけ変数
zdata<-rdirichlet(Nsample,c(1.2,0.3,0.1))

# 原点から３方向の単位ベクトル３つを持つ行列を作成
v3<-matrix(c(0,1,-sqrt(3)/2,-0.5,sqrt(3)/2,-0.5),3,2,byrow=TRUE)

# z0,z1,z2の３つ組値データを３角形座標に変換する
out<-matrix(0,Nsample,2)
for(i in 1:Nsample){
	out[i,]<-t(v3)%*%zdata[i,]
}
# z0 vs. z1, z1 vs. z2, z2 vs.z0のコプロット
plot(as.data.frame(zdata),xlim=c(0,1),ylim=c(0,1))

# 三角形座標へのプロット
plot(out,xlim=c(-1,1),ylim=c(-1,1))
# 三角形の辺を引く
for(i in 1:3){
	pre<-i
	post<-i+1
	if(post>3)post<-1
	segments(v3[pre,1],v3[pre,2],v3[post,1],v3[post,2])

}

• この話は２次元
• 多次元にできる
• ３次元、４要素にすれば、２SNPの４ハプロタイプ頻度が相当する。それを描くのが以下のソース
• ２SNPの４ハプロタイプの状態は正四面体上の位置に対応するので、連鎖平衡面がこの正四面体上に存在する

# ４方向と言えば、２SNPの４ハプロタイプ
CategoryVector <-
function (nc = 3) 
{
    df <- nc - 1
    d <- df + 1
    diagval <- 1:d
    diagval <- sqrt((df + 1)/df) * sqrt((df - diagval + 1)/(df - 
        diagval + 2))
    others <- -diagval/(df - (0:(d - 1)))
    m <- matrix(rep(others, df + 1), nrow = df + 1, byrow = TRUE)
    diag(m) <- diagval
    m[upper.tri(m)] <- 0
    as.matrix(m[, 1:df])
}

# ４ハプロタイプの頻度を作る
# 引数c(10,3,1)はz0の値が重くなりがちになるような重みづけ変数
Nsample<-10000
zdata<-rdirichlet(Nsample,c(0.1,0.01,0.02,0.1))

# 原点から４方向の単位ベクトル４つを持つ行列を作成(３次元空間)
v4<-CategoryVector(nc=4)

# h1,h2,h3,h4の４つ組値データを正四面体に変換する
out<-matrix(0,Nsample,3)
for(i in 1:Nsample){
	out[i,]<-t(v4)%*%zdata[i,]
}
# z0 vs. z1, z1 vs. z2, z2 vs.z0のコプロット
plot(as.data.frame(zdata),xlim=c(0,1),ylim=c(0,1))

# ３次元空間へのプロット

library(rgl)
xlim<-ylim<-zlim<-c(min(v4),max(v4))
plot3d(out[,1],out[,2],out[,3],xlim=xlim,ylim=ylim,zlim=zlim)