図示する

  • 久しぶりにSNPデータのR処理の話し
  • SNPデータからIBD推定としてz0,z1,z2を算出する(こちら)が、この3つの値は足し合わせると1になるので、自由度が2〜2次元平面に描ける。
  • z0,z1,z2をそれぞれ正三角形の頂点にとって描く
# サンプル数
Nsample<-100
# rdiriclhet()はz0,z1,z2からなる3つ組値データ(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))

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

# z0,z1,z2の3つ組値データを3角形座標に変換する
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])

}
  • この話は2次元
  • 多次元にできる
  • 3次元、4要素にすれば、2SNPの4ハプロタイプ頻度が相当する。それを描くのが以下のソース
  • 2SNPの4ハプロタイプの状態は正四面体上の位置に対応するので、連鎖平衡面がこの正四面体上に存在する
# 4方向と言えば、2SNPの4ハプロタイプ
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])
}

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

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

# h1,h2,h3,h4の4つ組値データを正四面体に変換する
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))

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

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