家系データシュミレーション

乱数 R 家系図 遺伝学

YN先生のサイト

アレル頻度とは
ハプロタイプは、アレルの数珠つなぎ(tandem: 縦につないだもの)

詳しくはこちら

N###############
### 1 
###############

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

###############
### 2 
###############

# 一様乱数からの発生
Pa<-runif(Na)
# アレル頻度が1/Na付近に固まるので…
# そうならないように道具をウェブ上から導入する
install.packages("MCMCpack")
library(MCMCpack)
# Na個の足して１になる乱数を１セット得る
Pa<-rdirichlet(1,rep(1,Na))
# できたPaの中身を見る
Pa
# 和が1担っていることを確認する
sum(Pa)

###############
### 3 
###############

# 同一染色体上にある多型の数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])))
}
NaS
PaS

hap1<-hap2<-rep(0,Ns) #父親由来、母親由来のハプロタイプ２つ用意する
for(i in 1:Ns){
  hap1[i]<-sample(1:NaS[i],1,prob=PaS[[i]])　#1〜NaS[i]番目までのアレルの中から１つ選ぶ。選び出される確率はPaS[i]
  hap2[i]<-sample(1:NaS[i],1,prob=PaS[[i]])
}
hap1
hap2

###############
### 4 
###############
#関数にする(毎回上を使うのはめんどくさい)

# 入力はマーカーのアレルNaS、アレル頻度PaSと作成するハプロタイプの本数N
MakeHaplotype<-function(NaS,PaS,N){
  # haplotypeはN x length(NaS)の行列に納める
  H<-matrix(0,N,length(NaS))
  for(i in 1:length(NaS)){
    H[,i]<-sample(1:NaS[i],N,prob=PaS[[i]])
  }
  H
}
# 使ってみる
N<-2 # 一人分の２本のハプロタイプ
Haps<-MakeHaplotype(NaS,PaS,N)
Haps

###############
### 5 
###############

# 組換えをシミュレーションするために、２本のハプロタイプを作る
# Ns個のマーカーの(Ns-1)個の「間」について、モルガン単位で遺伝的距離を定める
Dg<-runif(Ns-1)
# 組換えを起こして２つの組換え体を作ろう
# マーカー間ごとに、交叉がおきるかどうかを確率的に決める
rands.crossing<-runif(length(Dg))
# 交差がおきた箇所を取り出す
Cross<-which((Dg-rands.crossing)>0)
Cross
# 交差が起きたかどうかの情報を用いて組換え体を作る
hap1<-Haps[1,]
hap2<-Haps[2,]
rec.hap1<-hap1
rec.hap2<-hap2
if(length(Cross)>0){
  for(i in 1:length(Cross)){
    # 交差箇所より尾部について、rec.hap1,rec.hap2で交換する
    tmp.tail1<-rec.hap1[(Cross[i]+1):Ns]
    tmp.tail2<-rec.hap2[(Cross[i]+1):Ns]
    rec.hap1[(Cross[i]+1):Ns]<-tmp.tail2
    rec.hap2[(Cross[i]+1):Ns]<-tmp.tail1
  }
}
Cross
hap1
hap2
rec.hap1
rec.hap2

###############
### 6 
###############
#やっぱり関数化

# 2本のハプロタイプと、多型間遺伝的距離から、組換え体２本を返す
Rec.Haplotype<-function(hap1,hap2,d){
	Ns<-length(hap1)
	rands.crossing<-runif(length(d))
	# 交差がおきた箇所を取り出す
	Cross<-which((d-rands.crossing)>0)
	# 交差が起きたかどうかの情報を用いて組換え体を作る
	rec.hap1<-hap1
	rec.hap2<-hap2
	if(length(Cross)>0){
		for(i in 1:length(Cross)){
			# 交差箇所より尾部について、rec.hap1,rec.hap2で交換する
			tmp.tail1<-rec.hap1[(Cross[i]+1):Ns]
			tmp.tail2<-rec.hap2[(Cross[i]+1):Ns]
			rec.hap1[(Cross[i]+1):Ns]<-tmp.tail2
			rec.hap2[(Cross[i]+1):Ns]<-tmp.tail1
		}
	}
	matrix(c(rec.hap1,rec.hap2),byrow=TRUE,nrow=2)
}
# 使ってみる

rec.haps<-Rec.Haplotype(Haps[1,],Haps[2,],Dg)
rec.haps

###############
### 7 
###############


# 父・母のジェノタイプをランダムに作る
father<-MakeHaplotype(NaS,PaS,2)
mother<-MakeHaplotype(NaS,PaS,2)
father
mother
# 配偶子を作る
father.gamates<-Rec.Haplotype(father[1,],father[2,],Dg)
mother.gamates<-Rec.Haplotype(mother[1,],mother[2,],Dg)
# 父由来・母由来の配偶子をランダムに選んで子のジェノタイプを決める
from.father<-father.gamates[sample(1:2,1),]
from.mother<-mother.gamates[sample(1:2,1),]
offspring<-matrix(c(from.father,from.mother),byrow=TRUE,nrow=2)
offspring

father
mother
offspring

###############
### 8 
###############
#関数化

MakeOffspring<-function(father,mother,d){
  # 配偶子を作る
  father.gamates<-Rec.Haplotype(father[1,],father[2,],d)
  mother.gamates<-Rec.Haplotype(mother[1,],mother[2,],d)
  # 父由来・母由来の配偶子をランダムに選んで子のジェノタイプを決める
  from.father<-father.gamates[sample(1:2,1),]
  from.mother<-mother.gamates[sample(1:2,1),]
  offspring<-matrix(c(from.father,from.mother),byrow=TRUE,nrow=2)
  offspring  
}

# 使ってみる
offspring<-MakeOffspring(father,mother,Dg)
father
mother
offspring

###############
### 9 
###############

#「子は親よりも大きいID」をつけることに責任を持つ(よほど複雑な家系図でない限り)←プログラムを書きやすくするために

# 両親と一人っ子
M.Oyako.Trio<-matrix(c(0,0,0,0,1,2),byrow=TRUE,ncol=2)
M.Oyako.Trio
# １人とその両親とそのさらに両親
M.3gen<-matrix(c(0,0,0,0,0,0,0,0,1,2,3,4,5,6),byrow=TRUE,ncol=2)
M.3gen
# 一組のいとこ
# いとこペアの１人、その両親、その片方の両親(共有される祖父母)
# 共有された祖父母のもう一人の子、その配偶者、その子(いとこペアのもう一人)
M.cousins<-matrix(c(0,0,0,0,1,2,0,0,1,2,0,0,3,4,5,6),byrow=TRUE,ncol=2)
M.cousins 

###############
### 10 
###############

# M.cousins[,1]はM.cousings[,1]の第１列のこと
n.rows<-length(M.cousins[,1]) # nrows<-nrow(M.cousins)でも同じこと
# 人数=n.rows 分のジェノタイプは　マーカー数 x 2 行列が人数分あるので、それを納めるリストを作る
sim.G<-list()

# 第1行から第n.rows行まで、順番に処理するループ
for(i in 1:n.rows){
  # ここで処理される人はi行目の人
  # その人の両親は
  tmp.father<-M.cousins[i,1]
  tmp.mother<-M.cousins[i,2]
  # tmp.father,tmp.motherのジェノタイプは
  if(tmp.father==0){
    father.genotype<-MakeHaplotype(NaS,PaS,2)
  }else{
    # 親のジェノタイプはこのループ処理ですでに決まっている
    father.genotype<-sim.G[[tmp.father]]
  }
  if(tmp.mother==0){
    mother.genotype<-MakeHaplotype(NaS,PaS,2)
  }else{
    mother.genotype<-sim.G[[tmp.mother]]
  }
  # 両親のジェノタイプが決まったので、この人のジェノタイプを決める
  sim.G[[i]]<-MakeOffspring(father.genotype,mother.genotype,Dg)
}
sim.G
M.cousins

###############
### 11 
###############

# 引数はM.cousins,NaS,PaS,Dg
Simulate.Genotype<-function(M,NaS,PaS,Dg){
  # M.cousins[,1]はM.cousings[,1]の第１列のこと
  n.rows<-length(M.cousins[,1]) # nrows<-nrow(M.cousins)でも同じこと
  # 人数=n.rows 分のジェノタイプは　マーカー数 x 2 行列が人数分あるので、それを納めるリストを作る
  sim.G<-list()

  # 第1行から第n.rows行まで、順番に処理するループ
  for(i in 1:n.rows){
    # ここで処理される人はi行目の人
    # その人の両親は
    tmp.father<-M.cousins[i,1]
    tmp.mother<-M.cousins[i,2]
    # tmp.father,tmp.motherのジェノタイプは
    if(tmp.father==0){
      father.genotype<-MakeHaplotype(NaS,PaS,2)
    }else{
      # 親のジェノタイプはこのループ処理ですでに決まっている
      father.genotype<-sim.G[[tmp.father]]
    }
    if(tmp.mother==0){
      mother.genotype<-MakeHaplotype(NaS,PaS,2)
    }else{
      mother.genotype<-sim.G[[tmp.mother]]
    }
    # 両親のジェノタイプが決まったので、この人のジェノタイプを決める
    sim.G[[i]]<-MakeOffspring(father.genotype,mother.genotype,Dg)
  }
  sim.G
}
# 使ってみる
sim.G<-Simulate.Genotype(M.cousins,NaS,PaS,Dg)
sim.G

###############
### 12 
###############

# 染色体ごとにNaS,PaS,Dgを作る
# ここでは簡単のために、すべて同じものとする
num.Chromosomes<-22
NaS.list<-list()
PaS.list<-list()
Dg.list<-list()
for(i in 1:num.Chromosomes){
  NaS.list[[i]]<-NaS
  PaS.list[[i]]<-PaS
  Dg.list[[i]]<-Dg
}
sim.Gs<-list()
for(i in 1:num.Chromosomes){
  sim.Gs[[i]]<-Simulate.Genotype(M.cousins,NaS.list[[i]],PaS.list[[i]],Dg.list[[i]])
}

###############
### 13 
###############

# 入力値はM.cousins,NaS.list,PaS.list,Dg.list
Simulate.Genotype.Multi<-function(M,NaS.list,PaS.list,Dg.list){
  num.Chromosomes<-length(NaS.list)
  sim.Gs<-list()
  for(i in 1:num.Chromosomes){
    sim.Gs[[i]]<-Simulate.Genotype(M.cousins,NaS.list[[i]],PaS.list[[i]],Dg.list[[i]])
  }
  sim.Gs
}
# 使ってみる
sim.Gs<-Simulate.Genotype.Multi(M.cousins,NaS.list,PaS.list,Dg.list)


###############
### 14 
###############

# 繰り返し回数
Niter<-10
sim.Gs.list<-list()
for(i in 1:Niter){
  sim.Gs.list[[i]]<-Simulate.Genotype.Multi(M.cousins,NaS.list,PaS.list,Dg.list)
}

まとめ

●必要な変数
No. markers→genetic distance btw markerも必要
No. allels
Allele.freq
familytree

●simulationする
●あとはfamilyの数を増やしたり、減らしたり