遺伝的浮動

遺伝統計学の基礎―Rによる遺伝因子解析・遺伝子機能解析―

• こちらのブログの内容をまとめた上の本参照
• 遺伝的浮動のシミュレーション
  • 集団のサイズを固定:N
  • 次の世代へ残す遺伝子の最大数:k
  • 世代:T
  • 最初の変異アレル数:initial

• 行列Aは推移確率を要素にしている
  • 　　
  • 　　　
• 確率の計算結果を行列Dに入れてある
  • 　（アインシュタインの記法）
• 最後はエントロピーを計算してみた
  • 拡散してはじめは増加するが途中から減少する

N<-50 #集団のサイズ
k<-3    #残せる遺伝子の最大数
#k<-"infinity"  #残せる遺伝子が無限大の場合
T<-10   #世代数
initial <- 5    #最初に存在する変異アレル数


#推移確率の行列を作る
A<-matrix(0,N+1,N+1)
for (i in 0:N){
	for(j in 0:N){
		if (k == "infinity"){
			A[j+1,i+1]<-choose(N,j)*i^j*(N-i)^(N-j)/N^N
			}else{
				A[j+1,i+1]<- choose(i*k,j)*choose((N-i)*k,N-j)/choose(N*k,N)
				}
		
		}
	}
#apply(A,2,sum)  #確認

#各世代ごとの変異アレル数の存在確率のベクトル
v<-c(rep(0,N+1))
v[initial+1]<-1   #第１世代の分布

#各世代のベクトルvを並べて行列Dを作る
D<-c()
D<-cbind(D,v)
for (t in 2:T){
	v<-A%*%v   #次世代のベクトル
	D<-cbind(D,v)
	}
#apply(D,2,sum)  #確認
persp(D,theta=150,phi=30,shade=0.2,xlab="mutations",ylab="generations",zlab="probability",col="green")

#エントロピーを計算してみる
S<-apply(-D*log(D),2,sum)
plot(S)