生成確率の計算

• 組換えの計算をしていて
• このときに3lociの確率を計算した
• その時のソースは必要でない部分も含まれていたので今回省いた
  • 以下のパラメータも自由に決めることにする
    • ローカスの数：
    • 各ローカスに存在するアレル数：
    • 各ローカス間の組換え確率：
• 計算方法
  • ２倍体(hap1とhap2)と１倍体(hap)が与えらている
  • rsは長さ２のベクトルでfor文の中で値がかわっていく
  • ローカスを順に見ていって、子ハプロタイプ(hap)と一致するようなアレルの組み合わせがつくられる組換えが起こる確率を計算する
  • その確率と、ディプロタイプ(hap1とhap2)のどちらが一致しているか、の情報をrsが持っている
• 
• 
• 最初にhapとhap1とhap2が自由に与えられていると、こちらにあるようにDNA鑑定みたい..

l<-3  #ローカスの数

#各ローカスに存在するアレル数
k<-c()
for(i in 1:l){
k[i]<-sample(2:4,1)
}

hap<-hap1<-hap2<-c()
for (i in 1:l){
hap1[i]<-sample(1:k[i],1)
hap2[i]<-sample(1:k[i],1)
#hap[i]<-sample(1:k[i],1)      #子ハプロタイプも自由に決める場合
hap[i]<-sample(c(hap1[i],hap2[i]),1)    #必ずディプロタイプ(hap1,hap2)から生じている場合
}

v<-rbind(hap1==hap,hap2==hap)


#各ローカス間の組換え確率
#s<-runif(l-1)/2
s<-rep(0.1,length=l-1)
T<-rbind(s,1-s)    #このTは下の計算のためのもの

if(min(v[1,] + v[2,]) > 0){
rs<-c(v[1,1],v[2,1])/2	#最初のローカスが一致している確率
	for ( i in 1:(l-1)){
	rs<-(cbind(1-T[,i],T[,i])*v[,i+1]%*%t(c(1,1)))%*%rs   #順にローカスを増やしていって、一致している確率を計算する
	}
r<-sum(rs)
	
}else{r<-0}

print(v[1,])
print(v[2,])
print(r)