2011-05-07

マーカーごとに可能なディプロタイプ・ハプロタイプを選ぶ

R DNA鑑定

マーカーが相互に独立な場合には、個々のマーカーに関して確率・尤度を計算して、積を取ることができる
ゲノム全体にぱらぱらと置いたIdentifilerの場合には、この仮定でよい
マーカーごとに次のように考える
アレル数が $Na$ のとき、ディプロタイプのタイプ数は $Ng=\frac{Na(Na+1)}{2}$ ある
ディプロタイプの情報がなければ、すべてのディプロタイプをすべての個人が取りうるものとして考えることになります。人数がならば通りが考えられます
- これは、計算量が膨大です
ディプロタイプが既知のメンバー(家系情報行列の第５カラムが1)がいるときには、その分を減らせるが、家系が大きかったり、タイプが2,3のメンバーの数が複数のときには、計算量が十分大きく、現実的ではない
家系情報とディプロタイプ既知メンバーのディプロタイプ情報を用いて、ディプロタイプ不明メンバーの取りうるディプロタイプや、伝達(継承)情報が限定できる
そのルールは『子は両親の持つ２アレルの片方を受け継ぐ』こと
- これが以下のようにいろいろな制約をもたらす
  - 子は由来親の持っていないアレルは受け取れない
  - 子は両親から伝達される可能性のないアレルはディプロタイプに持ちえない
  - 子の由来親から受け取れないアレルの数がアレル総数-1に等しいとき、そのアレルは受け取るアレルとして確定する
  - 子が必ず持つアレルことがわかっているアレルのうち、父(母)由来でないアレルは、必ず母(父)からの伝達アレルである
  - 子が由来親から伝達することが確実なアレルが１つあるとき、それ以外のアレルが伝達されることはない
  - 子が由来親から伝達されることがあり得ないアレルの数がアレル総数-1に等しいとき、伝達される可能性のあるアレルは確実に伝達アレルであるとわかる
- この処理に必要なのは、家系情報とディプロタイプとアレル総数である
- 以下のソースは、長く、また汚く、返り値もごちゃごちゃした段階であるが、「ディプロタイプ」と「伝達ハプロタイプ」をそれぞれ出力すること、それを集合として扱うことと、ベクトルとして扱うこと、また、アレル名称で扱うことと、アレルにつけたシークエンスIDで扱うために、重複していることが主な理由で、その視点でみると、それほど複雑ではない

LimitDiplotypeZ<-function(p,g,A){
	unknown="0"
	ns<-length(p[,1])
	na<-length(A)
	#D ディプロタイプ
	#H ハプロタイプ
	#0 不明 1必ず持つ -1決して持たない
	D<-matrix(0,ns,na)
	H<-array(0,c(ns,2,na))
		# 確定ディプロタイプで代入
		# 確定ディプロタイプが示すアレルを持つことは確定し
		# 持たないアレルを持たないことも確定する
		# 持たないアレルが伝達することは決してない
		for(i in 1:ns){
			if(g[i,1]!=unknown){
				tmp<-rep(0,na)
				tmp[which(A==g[i,1])]<-1
				tmp[which(A==g[i,2])]<-1
				D[i,which(tmp==1)]<-1
				D[i,which(tmp==0)]<-(-1)
				H[i,1,which(tmp==0)]<-(-1)
				H[i,2,which(tmp==0)]<-(-1)
			}
		}

	# ループ
	loop<-TRUE
	while(loop){
		tmpH<-H
		tmpD<-D
		# 親子トリオ
		for(i in 1:ns){
			# 由来親の持っていないアレルは受け取れない
			# 由来親から受け取れないアレルの数がアレル総数-1に等しいとき、
			# そのアレルは受け取るアレルとして確定する
			if(p[i,2]!=0){# mother known
				H[i,1,which(D[p[i,2],]==-1)]<-(-1)
				if(length(which(D[p[i,2],]>=0))==1){
					#print(i)
					#print(H[i,1,])
					H[i,1,which(D[p[i,2],]==1)]<-1
					H[i,1,which(D[p[i,2],]!=1)]<-(-1)
					#print(H[i,1,])
				}
			}
			if(length(which(H[i,1,]==-1))==(na-1))H[i,1,which(H[i,1,]>(-1))]<-1
			if(length(which(H[i,1,]==1))==1)H[i,1,which(H[i,1,]!=1)]<--1
			D[i,which(H[i,1,]==1)]<-1
			if(p[i,3]!=0){ # father known
				H[i,2,which(D[p[i,3],]==-1)]<-(-1)
				if(length(which(D[p[i,3],]>=0))==1){
					H[i,2,which(D[p[i,3],]==1)]<-1
					H[i,2,which(D[p[i,3],]!=1)]<-(-1)
				}
			}
			if(length(which(H[i,2,]==-1))==(na-1))H[i,2,which(H[i,2,]>(-1))]<-1
			if(length(which(H[i,2,]==1))==1)H[i,2,which(H[i,2,]!=1)]<--1
			D[i,which(H[i,2,]==1)]<-1

			# 子の必ず持つアレルのうち、父(母)由来でないアレルは、必ず母(父)からの伝達である
			H[i,1,which(((D[i,]==1) * (H[i,2,]==-1))==1)]<-1
			H[i,2,which(((D[i,]==1) * (H[i,1,]==-1))==1)]<-1

			D[p[i,2],which(H[i,1,]==1)]<-1
			D[p[i,3],which(H[i,2,]==1)]<-1
			if(length(which(D[i,]==1))==2)D[i,which(D[i,]!=1)]<-(-1)
			
			
			# 両親から伝達される可能性のないアレルはディプロタイプに持ちえない
			for(i in 1:ns){
				D[i,which((H[i,1,]==-1) * (H[i,2,]==-1)==1)]<-(-1)
				if(length(which(D[i,]==-1))>=(na-1))D[i,which(D[i,]>(-1))]<-1
			}
		}
		if(isTRUE(all.equal(D,tmpD)) & isTRUE(all.equal(H,tmpH)))loop<-FALSE
		
	}
	dset<-dsetA<-list()
	dtypeset<-dtypesetA<-list()
	hset1<-hset2<-hset1A<-hset2A<-list()
	for(i in 1:length(D[,1])){
		tmpd1<-as.set(which(D[i,]==1))
		tmpd0<-as.set(which(D[i,]==0))
		tmpd1A<-as.set(A[which(D[i,]==1)])
		tmpd0A<-as.set(A[which(D[i,]==0)])
		if(set_is_empty(tmpd1)){
			#dset[[i]]<-tmpd0*tmpd0
			dset[[i]]<-set_combn(tmpd0,1) | set_combn(tmpd0,2)
			#dsetA[[i]]<-tmpd0A*tmpd0A
			dsetA[[i]]<-set_combn(tmpd0A,1) | set_combn(tmpd0A,2)

		}else{
			if(set_is_empty(tmpd0)){
				dset[[i]]<-set(tmpd1)
				dsetA[[i]]<-set(tmpd1A)
			}else{
				#dset[[i]]<-tmpd1*tmpd0
				dset[[i]]<-dsetA[[i]]<-set()
				tmpd01<-tmpd0 | tmpd1
				for(j in tmpd01){
					dset[[i]]<-dset[[i]] | set(tmpd1 | j)

				}
				tmpd01A<-tmpd0A | tmpd1A
				for(j in tmpd01A){
					dsetA[[i]]<-dsetA[[i]] | set(tmpd1A | j)
				}
				#dset[[i]]<-as.set(sapply(set_combn(tmpd0,1),"|",tmpd1))
				#dsetA[[i]]<-tmpd1A*tmpd0A
				#dsetA[[i]]<-as.set(sapply(set_combn(tmpd0A,1),"|",tmpd1A))
			}
			
		}
		#dset[[i]]<-as.set(which(D[i,]>=0))
		hset1[[i]]<-as.set(which(H[i,1,]>=0))
		hset2[[i]]<-as.set(which(H[i,2,]>=0))
		dtypeset[[i]]<-as.set(hset1[[i]]*hset2[[i]])
		#dsetA[[i]]<-as.set(A[which(D[i,]>=0)])
		hset1A[[i]]<-as.set(A[which(H[i,1,]>=0)])
		hset2A[[i]]<-as.set(A[which(H[i,2,]>=0)])
		dtypesetA[[i]]<-as.set(hset1A[[i]]*hset2A[[i]])
	}
	# 子が持つアレルを１本も持たないディプロタイプを親は持たない
	for(i in 1:length(D[,1])){
		tmpdset<-set()

		for(j in dset[[i]]){
			offsOK<-TRUE
			oyajudge<-FALSE

			oyanashi<-TRUE
			for(k in 1:length(p[,1])){
				if(p[k,2] == i | p[k,3] == i){
				offsOKperchild<-FALSE

					oyajudge<-TRUE
					oyanashi<-FALSE
						for(l in dset[[k]]){
							if(!set_is_empty(j & l))offsOKperchild<-TRUE
							#print(j)
							#print(l)
							#print(offsOKperchild)
							#print("---")
						}
						if(!offsOKperchild)offsOK<-FALSE

				}

			}
			hapsOK<-FALSE
			if(!(set_is_empty(j & hset1[[i]])) & !(set_is_empty(j & hset2[[i]])) &
			length(j & (hset1[[i]] | hset2[[i]]) )==length(j) )hapsOK<-TRUE
			if(offsOK & hapsOK){
				#oyajudge=TRUE
				tmpdset<-tmpdset | set(j)
			}
			if(oyanashi & hapsOK){
				oyajudge=TRUE
				tmpdset<-tmpdset | set(j)
			}
		}
		#if(oyajudge){
			dset[[i]]<-tmpdset
		#}
		
	}
	for(i in 1:length(D[,1])){
		tmpdset<-set()

		for(j in dsetA[[i]]){
			offsOK<-TRUE
			oyajudge<-FALSE

			oyanashi<-TRUE
			for(k in 1:length(p[,1])){
				if(p[k,2] == i | p[k,3] == i){
				offsOKperchild<-FALSE

					oyajudge<-TRUE
					oyanashi<-FALSE
						for(l in dsetA[[k]]){
							if(!set_is_empty(j & l))offsOKperchild<-TRUE
							#print(j)
							#print(l)
							#print(offsOKperchild)
							#print("---")
						}
						if(!offsOKperchild)offsOK<-FALSE

				}

			}
			hapsOK<-FALSE
			if(!(set_is_empty(j & hset1A[[i]])) & !(set_is_empty(j & hset2A[[i]])) &
			length(j & (hset1A[[i]] | hset2A[[i]]) )==length(j) )hapsOK<-TRUE
			if(offsOK & hapsOK){
				#oyajudge=TRUE
				tmpdset<-tmpdset | set(j)
			}
			if(oyanashi & hapsOK){
				oyajudge=TRUE
				tmpdset<-tmpdset | set(j)
			}
		}
		#if(oyajudge){
			dsetA[[i]]<-tmpdset
		#}
	}
	dlist<-dlistA<-list()
	
	#for(i in 1:length(dtypeset)){
	for(i in 1:length(dset)){
		cnt<-1
		dlist[[i]]<-list()
		#for(j in dtypeset[[i]]){
		for(j in dset[[i]]){
			dlist[[i]][[cnt]]<-as.set(j)
			cnt<-cnt+1
		}
	}
	#for(i in 1:length(dtypesetA)){
	for(i in 1:length(dsetA)){
		cnt<-1
		dlistA[[i]]<-list()
		#for(j in dtypesetA[[i]]){
		for(j in dsetA[[i]]){
			dlistA[[i]][[cnt]]<-as.set(j)
			cnt<-cnt+1
		}
	}
	list(D=D,H=H,dlist=dlist,dlistA=dlistA,dset=dset,dsetA=dsetA,dtypeset=dtypeset,dtypesetA=dtypesetA,
	hset1=hset1,hset2=hset2,hset1A=hset1A,hset2A=hset2A)
}

LDZout<-LimitDiplotypeZ(p,g,A)