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

• 1. アレル、ジェノタイプ
  • アレルとアレル頻度

# あるマーカー
# アレルの数
Na<-4
# アレル頻度
Pa<-c(0.3,0.2,0.1,0.4)
# Pa<-c(3,2,1,4)のようにして、あとで頻度の和が１になるようにしてもよい
# 念のため、アレル頻度の和が１になるように補正する
Pa<-Pa/sum(Pa)

• • • アレル頻度を適当に決めたいとき

# 一様乱数からの発生
Pa<-runif(Na)
# アレル頻度が1/Na付近に固まるので…
# そうならないように道具をウェブ上から導入する
install.packages("MCMCpack")
library(MCMCpack)
# Na個の足して１になる乱数を１セット得る
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