2ローカスの組換え

• こちらとこちらの考察をふまえて
• 組換えによってハプロタイプにおいて連鎖平衡に近づく様子をみてみることに
  • このときも時系列を考えていた
• 方法は以前のarrayを使ったもの
  • こちらにあるようにarray D の計算は変更してある
  • 同じくこちらとほとんど一緒のことのような気がするが...
• 仮定
  • ハプロタイプの同時分布:H,H2
  • 連鎖平衡
    • 連鎖平衡に近づく様子のシミュレーションなので
    • 平衡状態のハプロタイプの頻度:H0
  • HW平衡
    • ハプロタイプからディプロタイプの頻度を計算するのにHWEを使った
    • ディプロタイプの頻度:D
• 世代ごとに連鎖平衡のとき(H0)との差を見ている
• 平衡を仮定しないときはどうなるか
  • ミクロ vs マクロ
  • マクロでみれば
    • ハプロタイプの存在頻度の偏り
    • ディプロタイプの存在頻度の偏り
  • ミクロでみれば？
    • セレククション？　正？？負？？
    • そもそも局所平衡？

#1倍体、2倍体の存在頻度の計算

N1<-4    #ローカス1のアレル数
N2<-4    #ローカス2のアレル数

H<-matrix(runif(N1*N2),N1,N2)
H<-H/sum(H)      #これがハプロタイプ上のローカスにおける"同時分布"
H2<-matrix(0,N1,N2)  

#連鎖平衡を使って平衡に達したときの頻度計算する
f1<-apply(H,1,sum)    #ローカス1のアレル頻度
f2<-apply(H,2,sum)    #ローカス2のアレル頻度
H0<-f1%*%t(f2)        #平衡時の２倍体の頻度

persp(H-H0,col="green",theta=30,phi=30,xlab="locus1",ylab="locus2",zlab="probability(-0.1,0.1)",zlim=c(-0.1,0.1))

#2倍体からハプロタイプが発生する確率を要素にしたarrayを作製

s<-0.3     #組換え率

C<-array(0,c(N1,N2,N1,N2,N1,N2))
for (m in 1:N1){
	for(n in 1:N2){
		C[m,n,,,m,n]<-c(1/2)
		for(i in 1:N1){
			for(j in 1:N2){
				if (i != m && j != n){
					C[m,j,i,n,m,n]<-s/2
					C[m,n,i,j,m,n]<-(1-s)/2
					}
				}
			}
		}
	}
C<-C+aperm(C,c(3,4,1,2,5,6))

#次世代のハプロタイプの存在頻度の計算
for (t in 1:15){
D<- array(c(H)%*%t(c(H)),c(N1,N2,N1,N2))  #HWEの仮定をミクロな系で使っている（？）

for (m in 1:N1){
	for (n in 1:N2){
		H2[m,n]<-C[,,,,m,n]%*%D     #ここでH2の要素が計算できる
		}
	}
#H2     #次世代のハプロタイプの頻度をあらわす行列

persp(H2-H0,col="green",theta=30,phi=30,xlab="locus1",ylab="locus2",zlab="probability(-0.1,0.1)",zlim=c(-0.1,0.1))
H<-H2
}