2ローカスの組換え

• きのうの続き
• 仮定の条件を弱める
  • ２倍体の頻度を計算するときにHWEを仮定していた
  • HWDにする
    • つまり２倍体の存在頻度に偏りをもたせる
• そのためにHWDというarrayを作る
  • これは 0と1 でできたarray
  • HWDをDにかけることで特定のディプロタイプが存在しないようにした
    • 本当は{0,1}でなくてもよいはずだが..

HWD<-array(sample(c(0,1),size=N1^2*N2^2,replace=T,prob=c(0.2,0.8)),c(N1,N2,N1,N2))

• これを使ってみると...
  • 前回と異なる分布が得られた
• これがこの場合の安定分布なのだろう
  • ただし分かっていないのが
    • 安定分布が何種類あるか
    • どの安定分布に近づくか
    • それを決める初期（？）条件は何か

#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.5,0.5)",zlim=c(-0.5,0.5))

#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)

HWD<-array(sample(c(0,1),size=N1^2*N2^2,replace=T,prob=c(0.2,0.8)),c(N1,N2,N1,N2))        #このarrayが 0と1でできている
#次世代のハプロタイプの存在頻度の計算
for (t in 1:50){
D<- array(c(H)%*%t(c(H)),c(N1,N2,N1,N2)) 
D<- D * HWD       #特定のディプロタイプが存在しなくなる
D<-D/sum(D)      #和を１に直す
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.5,0.5)",zlim=c(-0.5,0.5))
H<-H2
}