2011-01-16

IBDをトラッキング

R IBD

添付のような家系図があるとする（第2、第3世代のみが描かれている。以下のソースでは第1世代の両親を加えてある）
ここで、22常染色体の交叉組換えをランダムにシミュレーションしてやる
発端者と責任座位を与えたうえで、責任座位に関してIBDにあるかどうかを判定してやろう

scale<-1 # 染色体の長さの拡縮スケール
nch<-22 # 常染色体本数
# その長さ（ほぼ実数）
chromLen<-round(c(247249719,242951149,199501827,191273063,180857866,
170899992,158821424,146274826,140273252,135374737,134452384,
132349534,114142980,106368585,100338915,88827254,78774742,
76117153,63811651,62435964,46944323,49691432,154913754,57772954)/1,0)
# 個人を作る
# 個人はID登録、親IDを2つ登録できる
# transmissionは2本のハプロタイプの先頭塩基が親の親（祖父母）のどちら由来かを指定し
# その後、交叉の起きた箇所をリストで持つことでIBD情報を積む
MakeIndividual<-function(id,parents=NULL,nch=22,chrLength=NULL,dip=2,phenotype=NULL){
	if(is.null(chrLength))chrLength<-c(247249719,242951149,199501827,191273063,180857866,170899992,158821424,146274826,140273252,135374737,134452384,132349534,114142980,106368585,100338915,88827254,78774742,76117153,63811651,62435964,46944323,49691432,154913754,57772954)

	transmission<-list()
	for(i in 1:nch){
		transmission[[i]]<-NULL
	}
	list(id=id,parents=parents,transmission=transmission,nch=nch,chrLength=chrLength,dip=dip,phenotype=phenotype)
}
# 家系の構成人数
Nm<-27
# 人数分の個人のリストで家系メンバーを束ねる
Members<-list()
for(i in 1:Nm){
	Members[[i]]<-MakeIndividual(id=i,phenotype=0)
}
# 伝達情報がないときに、ランダムに伝達情報を作成する
# rは１塩基あたりの交叉確率
MakeRandomTransmission<-function(member,r=10^(-6)){
	#rは１塩基あたりの交叉確率
	ret<-list()
	for(i in 1:member$nch){
		ret[[i]]<-list()
		for(j in 1:2){
			initparent<-sample(1:2,1)
			numcross<-member$chrLength[i]*r
			meancross<-member$chrLength[i]/numcross
			rexps<-rexp(rpois(1,numcross),1/meancross)
			rexps<-cumsum(rexps)
			ret[[i]][[j]]<-list(initparent=initparent,cross=rexps)
		}
	}
	ret
}
# 親子関係を3列行列で作る
# 家系図から地道に作ったもの
Oyako<-matrix(c(3,1,2,4,1,2,5,1,2,6,1,2,7,1,2,8,1,2,15,3,9,16,3,9,17,4,10,18,4,10,19,5,11,20,5,11,22,6,12,23,6,12,24,6,12,25,7,13,26,7,13,27,8,14),ncol=3,byrow=TRUE)
# 交叉確率を指定する
# 染色体長のスケールに応じて交叉確率もスケーリングする
r<-5*10^(-9)*10^(scale)
for(i in 1:length(Oyako[,1])){
	ko<-Oyako[i,1]
	oya1<-Oyako[i,2]
	oya2<-Oyako[i,3]
	Members[[ko]]$parents=list(Members[[oya1]],Members[[oya2]])
	Members[[ko]]$transmission<-MakeRandomTransmission(Members[[ko]],r=r)
}
# 第３メンバーについて１番染色体の交叉箇所を確認してみる
length(Members[[3]]$transmission[[1]][[2]])
# IBD関係を確認する
Niter<-10
for(i in 1:Niter){
	# 発端者を決める
	#resInd<-sample(1:Nm,1)
	resInd<-3
	# 責任座位を与える
	resCh<-sample(1:nch,1)
	resLoc<-sample(1:chromLen[resCh],1)
	resdip<-1
	tmp<-matrix(0,Nm,2)
	for(j in 1:Nm){
# 子の1つ目の相同染色体、2つ目のそれがどちらの親からのそれかをシャッフルする
		dipshuffle<-sample(1:2,2) 
		for(k in 1:2){
			loopjudge<-TRUE
			currentID<-j
			currenthap<-dipshuffle[k]
			while(loopjudge){
				if(is.null(Members[[currentID]]$parents)){
					loopjudge<-FALSE
				}else{
					tmptrans<-Members[[currentID]]$transmission[[resCh]][[currenthap]]
					numcross<-length(which(tmptrans[[2]]<resLoc))
					tmp[j,k]<-Members[[currentID]]$parents[[currenthap]]$id
					currenthap<-(tmptrans[[1]]+(numcross%%2))%%2+1
					currentID<-tmp[j,k]
				}
			}
		}
		
		
	}
print(tmp) # IBDを示す
}