２ローカスの組換え その２
- こちらの改良版
- コメントにいただいたようにハプロタイプをベクトルの要素にするとわかりにくいので行列にすることに
- やはり２ローカス
- 改良点はハプロタイプを行列表記にしたこと
- 行列とarrayで計算する
- 仮定など（前回と文字が異なる）
- 目的はHからH2を計算すること
- C:このarrayを使って計算する
#1倍体、2倍体の存在頻度の計算 N1<-3 #ローカス1のアレル数 N2<-3 #ローカス2のアレル数 H<-matrix(0,N1,N2) #ハプロタイプの存在頻度：これを行列表記にしてある(N1×N2 2次元) H[1,1]<-1/2 #最初のハプロタイプの頻度 H[2,2]<-1/4 H[3,3]<-1/4 B<-array(0,c(N1,N2,N1,N2)) #二倍体の存在頻度B:これはarray (N1×N2×N1×N2の4次元) for (i in 1:N1){ for(j in 1:N2){ for (k in 1:N1){ for (l in 1:N2){ B[i,j,k,l]<-H[i,j]*H[k,l] } } } }
#2倍体からハプロタイプが発生する確率を要素にしたarrayを作製 s<-0.3 #組換え率 C<-array(0,c(N1,N2,N1,N2,N1,N2)) #このarrayは前の記事のB_kをひとまとめにしたもの （6次元） 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)) #これでCができる 転置を使って書く必要はないが使うと短く書ける
#次世代のハプロタイプの存在頻度の計算 H2<-matrix(0,N1,N2) #H2は次世代のハプロタイプの存在頻度 (2次元) for (m in 1:N1){ for (n in 1:N2){ H2[m,n]<-C[,,,,m,n]%*%B #ここでH2の要素が計算できる } } H2 #次世代のハプロタイプの頻度をあらわす行列 persp(H2,col="green",theta=30,phi=30,xlab="locus1",ylab="locus2",zlab="probability",zlim=c(0,1))
- Cを作るところ
- Cは6次元
- 前の４元が次世代のハプロタイプ(m,n)が二倍体からできるときの確率をならべた4次元のarray
- 後ろの2元は前の4次元のarrayを番号付けするためのもの
- 前回の方法でも転置を使った(こちらの記事のBを作るところ)
- 使わない方が分かりやすいので以下に掲載
- 上の３つのスクリプトの真ん中と同じこと
C<-array(0,c(N1,N2,N1,N2,N1,N2)) s<-0.3 for (m in 1:N1){ #次世代のハプロタイプが(m,n) for(n in 1:N2){ C[m,n,m,n,m,n]<-1 #２倍体が(m,n),(m,n)のときハプロタイプ(m,n)ができる確率は1 for(i in 1:N1){ for(j in 1:N2){ if (i != m && j != n){ #上の場合以外について(i != m, j !=n としておいて) C[i,n,m,n,m,n]<-1/2 #2倍体が(i,n),(m,n) なので ハプロタイプ(m,n)ができる確率は1/2 C[m,j,m,n,m,n]<-1/2 C[m,n,i,n,m,n]<-1/2 C[m,n,m,j,m,n]<-1/2 C[m,j,i,n,m,n]<-s/2 #mとnが異なる染色体にあるので、組換えが起きると(m,n)が作られる C[i,n,m,j,m,n]<-s/2 C[m,n,i,j,m,n]<-(1-s)/2 #mとnが同じ染色体にあるので、組換えがおこらなければ(m,n)が作られる C[i,j,m,n,m,n]<-(1-s)/2 } } } } }
