家系データをシミュレーションで作る1

  • 1. アレル、ジェノタイプ
    • アレルとアレル頻度
# あるマーカー
# アレルの数
Na<-4
# アレル頻度
Pa<-c(0.3,0.2,0.1,0.4)
# Pa<-c(3,2,1,4)のようにして、あとで頻度の和が1になるようにしてもよい
# 念のため、アレル頻度の和が1になるように補正する
Pa<-Pa/sum(Pa)
      • アレル頻度を適当に決めたいとき
# 一様乱数からの発生
Pa<-runif(Na)
# アレル頻度が1/Na付近に固まるので…
# そうならないように道具をウェブ上から導入する
install.packages("MCMCpack")
library(MCMCpack)
# Na個の足して1になる乱数を1セット得る
Pa<-rdirichlet(1,rep(1,Na))
# できたPaの中身を見る
Pa
# 和が1担っていることを確認する
sum(Pa)

# 同一染色体上にある多型の数Ns
Ns<-4
# Ns個のマーカーのそれぞれのアレルの数をNaS
NaS<-c(3,4,5,2)
# アレル頻度をPaS
PaS<-list()

for(i in 1:Ns){
  PaS[[i]]<-c(rdirichlet(1,rep(1,NaS[i])))
  #PaS[[i]]<-c(runif(NaS[i]))
  #PaS[[i]]<-PaS[[i]]/sum(PaS[[i]])
}
NaS
PaS

hap1<-hap2<-rep(0,Ns)
for(i in 1:Ns){
  hap1[i]<-sample(1:NaS[i],1,prob=PaS[[i]])
  hap2[i]<-sample(1:NaS[i],1,prob=PaS[[i]])
}
hap1
hap2