１倍体で場合分け ソース ベータ版

• 依存パッケージ

# packages
# 依存パッケージ
library(kinship) # 家系図を描く
library(MCMCpack) # １倍体場合分けのみなら不要かも
library(gtools) # 場合分け　１倍体場合分けのみなら不要かも
library(sets) # 集合
library(paramlink) # 連鎖解析 核家族 # １倍体場合分けのみなら不要かも
library(rSymPy) # 数式演算

• ハプロタイプグラフを作る
  • ハプロタイプの親子関係を作る(複数のグラフが混在)

MakeHaplotypeGraph2<-function(p){
	ns<-length(p[,1])
	ret<-matrix(0,ns*2,2)
	# 2*i-1,2*iはi番個人の母方と父方の染色体
	for(i in 1:ns){
		tmpmom<-p[i,2]
		tmpdad<-p[i,3]
		ret[2*i-1,1]<-2*tmpmom-1
		ret[2*i-1,2]<-2*tmpmom
		ret[2*i,1]<-2*tmpdad-1
		ret[2*i,2]<-2*tmpdad
	}
	ret[which(ret<0)]<-0
	ret
}

• • ハプロタイプグラフを連結グラフ別に名前をつける

SeparateGraphs<-function(hG){
	ns<-length(hG[,1])
	assigned<-rep(0,ns)
	cnt<-1
	for(i in 1:length(hG[,1])){
		if(hG[i,1]!=0){
			tmpself<-assigned[i]
			tmpp1<-assigned[hG[i,1]]
			tmpp2<-assigned[hG[i,2]]
			M<-max(tmpself,tmpp1,tmpp2)
			if(M==0){
				M<-cnt
				cnt<-cnt+1
				assigned[c(i,hG[i,1],hG[i,2])]<-M
			}else{
				#M<-max(tmpself,tmpp1,tmpp2)
				trioID<-c(i,hG[i,1],hG[i,2])
				trio<-c(tmpself,tmpp1,tmpp2)
				for(j in 1:length(trio)){
					if(trio[j]!=0){
						assigned[which(assigned==trio[j])]<-M
					}else{
						assigned[trioID[j]]<-M
					}
				}
			}
			#print(assigned)
			
		}
	}
	#assigned
	G<-list()
	cnt<-1
	hG2<-cbind(1:ns,hG)
	for(i in 1:max(assigned)){
		if(length(which(assigned==i))!=0){
			G[[cnt]]<-hG2[assigned==i,]
			cnt<-cnt+1
		}
		
	}
	G
}

• 伝達関係の場合分けを列挙

MakeDecsendPattern<-function(hG){
	Alt<-which(hG[,1]!=0)
	#Alt<-1:length(hG[,1])
	# その数
	nAlt<-length(Alt)
	# 選択の場合を列挙
	s<-expand.grid(rep(list(c(1,2)),nAlt))
	for(i in 1:(2^nAlt)){
		tmp<-matrix(0,nAlt,2)
		for(j in 1:nAlt){
			tmp[j,1]<-Alt[j]
			tmp[j,2]<-hG[Alt[j],s[i,j]]
		}
	}
	Ns<-length(hG[,1])/2
	Ls<-length(s[,1])
	s2<-matrix(0,Ls,2*Ns)
	cnt<-1
	for(i in 1:Ns){
		if(p[i,2]==0){
			s2[,i*2-1]<-1
			s2[,i*2]<-2
		}else{
			s2[,i*2-1]<-s[,cnt]
			cnt<-cnt+1
			s2[,i*2]<-s[,cnt]
			cnt<-cnt+1
		}
	}

	s2

}

• 二つのアレルを父方・母方に振り分ける場合分け

# ジェノタイプknownの個人の２アレルを
# paternal/maternalに割り付ける場合分けを列強する
# ホモの場合も２パターンを作る方が
# 計算は多くなるが、間違いが少なくて済みそうだ…
MakePatMatPattern<-function(g){
	GenotypePlus<-which(g[,1]!=0)

	heteros<-rep(2,length(GenotypePlus))
	heteros[which(g[GenotypePlus,1]==g[GenotypePlus,2])]<-1
	heterolist<-list()

	for(i in 1:length(heteros)){
		heterolist[[i]]<-1:heteros[i]
	}
	s3<-expand.grid(heterolist)
	s3
}

• 母方・父方の場合分けに応じてアレルを割り当て

AssignHaplotype<-function(g,s){
	Ns<-length(g[,1])
	hapg<-rep(0,Ns*2)
	GenotypePlus<-which(g[,1]!=0)
	for(j in 1:length(GenotypePlus)){
		#print(g[GenotypePlus[j],unlist(s[j]]))
		hapg[GenotypePlus[j]*2-1]<-g[GenotypePlus[j],unlist(s[j])]
		hapg[GenotypePlus[j]*2]<-g[GenotypePlus[j],(unlist(s[j])-1.5)*(-1)+1.5]
	}
	hapg

}

• 伝達木(林)の作成
  • ハプロタイプグラフ情報と伝達パターンを決めれば、複数の木(林)になる

SelectTrees<-function(sepG,v){
	ret<-list()
	for(i in 1:length(sepG)){
		tmp<-sepG[[i]][,1:2]
		for(j in 1:length(tmp[,1])){
			tmp[j,2]<-sepG[[i]][j,v[tmp[j,1]]+1]
		}
		#print(tmp)
		tmp2<-rep(0,max(tmp))
		cnt<-1
		for(j in 1:length(tmp[,1])){
			if(tmp[j,1]*tmp[j,2]!=0){
				M<-max(tmp2[tmp[j,1]],tmp2[tmp[j,2]])
				#print(M)
				if(M==0){
					M<-cnt
					cnt<-cnt+1
					tmp2[tmp[j,1]]<-tmp2[tmp[j,2]]<-M
				}else{
					tmp2[tmp[j,1]]<-tmp2[tmp[j,2]]<-M
				}
			}
		}
		ret[[i]]<-tmp2[tmp[,1]]
	}
	ret
}

• 林のメンデルチェックと木の上のアレルを抽出
  • 林の木の上のジェノタイプがすべて同じかどうかのチェックをし、木のアレル、木の「集団からの直接の子」の数を数える

# 木全体のメンデルチェック
MendelCheck<-function(sTrees,sepGraphs,hapg){
	ret<-TRUE
	ret2<-c()
	ret3<-c()
	for(i in 1:length(sTrees)){ # グラフごとのループ
	#print("---")
	#print(max(sTrees[[i]]))
		for(j in 1:max(sTrees[[i]])){ # グラフにある木をループ
			# 木にある、"0"でないアレルはすべて同一である必要がある
			nodes<-sepGraphs[[i]][,1][which(sTrees[[i]]==j)]
			#print(hapg)
			#print(nodes)
			numanc<-length(which(sepGraphs[[i]][,2][which(sTrees[[i]]==j)]==0))
			tmp<-hapg[nodes][which(hapg[nodes]!="0")]
			#print(tmp)
			#print("length set")
			#print(length(as.set(tmp)))
			if(length(as.set(tmp))!=1){
				if(length(as.set(tmp))>1){
					ret<-FALSE
				}
				#ret<-FALSE
				ret2<-c(ret2,"0")
				ret3<-c(ret3,1)

				#return(list(mendel=ret,alleles=ret2,nAncestor=ret3))
			}else{
				#ret<-TRUE
				ret2<-c(ret2,hapg[nodes][which(hapg[nodes]!="0")][1])
				ret3<-c(ret3,numanc)
			}
		}
	}
	return(list(mendel=ret,alleles=ret2,nAncestor=ret3))
}

• マーカーごとに確率計算
  • メンデルチェックを通った場合というのは、「あり得る」場合
  • その確率は、「集団からの直接の子」が木のアレルを持つ確率と、場合分けの取り分による
  • すべての場合を加算する
  • 計算式も作る

ProbPerMarker<-function(p,g,As,Ps,hG,sepGraphs,s2,Expression=TRUE){
	cnt<-0
	s3<-MakePatMatPattern(g)
	if(Expression){
		Vars<-list()
		for(i in 1:length(As)){
			Vars[[i]]<-Var(letters[i])
			#print(Vars[[i]])
		}
		retExpression<-Var("0")
	}

	retProb<-0
	for(i in 1:length(s3[,1])){
		# Pat/Mat割り付けパターンでハプロタイプアレルを確定する
		hapg<-AssignHaplotype(g,s3[i,])
		for(j in 1:length(s2[,1])){
			# sepGraphsと伝達パターンから、木を作って取り出す
			sTrees<-SelectTrees(sepGraphs,s2[j,])
			# 木全体のメンデルチェック
			MendelOK<-MendelCheck(sTrees,sepGraphs,hapg)
			
			# MendelOK$mendelがTRUEなら確率加算
			if(MendelOK$mendel){
				#print(MendelOK)
				cnt<-cnt+1
				tmpexp<-Var(1/length(s2[,1]))
				
				tmpProb<-rep(0,length(MendelOK$alleles))
				for(k in 1:length(tmpProb)){
					if(!MendelOK$alleles[k]=="0"){
						aid<-which(As==MendelOK$alleles[k])
						#print(MendelOK$alleles[k])
						#print(aid)
						if(Expression){
							tmptmpexp<-Vars[[aid]]
							if(MendelOK$nAncestor[k]>1){
								for(l in 2:MendelOK$nAncestor[k]){
									tmptmpexp<-tmptmpexp*Vars[[aid]]
								}
							}
							tmpexp<-tmpexp*tmptmpexp
						}
						
						tmpProb[k]<-Ps[which(As==MendelOK$alleles[k])]
						
					}else{
						tmpProb[k]<-1
					}
					
				}
				if(Expression){
					retExpression<-retExpression+tmpexp
				}
				
				#print(tmpProb)
				#print(prod(tmpProb^MendelOK$nAncestor)/length(s2[,1]))
				retProb<-retProb+prod(tmpProb^MendelOK$nAncestor)/length(s2[,1])
			}
			if(!MendelOK$mendel){
				#print(sTrees)
				#print(MendelOK)
			}
				
		}
	}
	#print(cnt)
	if(!Expression) retExpression=Var("0")
	list(prob=retProb,express=retExpression)
}

• マーカーごとに尤度比を計算する
  • ２つの仮説の尤度を計算してその比を取る

CalcLikeRatioPerMarker<-function(p,g1,g2,searched,As,Ps,hG,sepGraphs,s2,Expression){
	retProb1<-ProbPerMarker(p,g1,As,Ps,hG,sepGraphs,s2,Expression)
	#print("-----")
	retProb2<-ProbPerMarker(p,g2,As,Ps,hG,sepGraphs,s2,Expression)
	expres1<-retProb1$express
	expres2<-retProb2$express
	#print(expres1)
	#print(expres2)
	ProbGenPop<-1
	genExp<-Var(1)
	for(i in 1:length(searched)){
		tmpallele1<-g2[searched[i],1]
		tmpallele2<-g2[searched[i],2]
		P1<-Ps[which(As==tmpallele1)]
		P2<-Ps[which(As==tmpallele2)]
		#print(P1)
		#print(P2)
		tmp<-P1*P2
		tmpVar1<-Var(letters[which(As==tmpallele1)])
		tmpVar2<-Var(letters[which(As==tmpallele2)])
		tmpgenExp<-tmpVar1*tmpVar2
		if(P1!=P2){
			tmp<-2*tmp
			tmpgenExp<-2*tmpgenExp
		}
		ProbGenPop<-ProbGenPop*tmp
		genExp<-genExp*tmpgenExp
	}
	#print(genExp)
	#print(ProbGenPop)
	retProb1x<-retProb1$prob*ProbGenPop
	#print(retProb1$prob)
	#print(retProb2$prob)
	list(P1=retProb1x,P2=retProb2$prob,Pratio=retProb2$prob/retProb1x,expression=expres2/(expres1*genExp))
}

• すべてのマーカーについて繰り返して結果を格納する
  • pは家系情報
  • gsはジェノタイプ情報
  • searchedは被捜索者リスト
  • poolidはsearchedがGpoolのどのジェノタイプに対応させるか
  • Allelesはアレル名
  • Probsはアレル頻度
  • Express　は計算式を出すか出さないかのオプション

CalcLikeRatioAllMarker<-function(p,gs,searched,poolid,Gpool,Alleles,Probs,Express=TRUE){
	CLRAM<-list()
	ped<-MakePedigreeFromFamilyInfo(p)
	plot(ped)
	# ハプロタイプグラフを作る
	hG<-MakeHaplotypeGraph2(p)

	# ハプロタイプグラフを分ける
	sepGraphs<-SeparateGraphs(hG)

	# 伝達の場合分けを列挙する

	s2<-MakeDecsendPattern(hG)
	
	for(na in 1:length(Alleles)){
		As<-Alleles[[na]]
		Ps<-Probs[[na]]
		g<-gs[,,na]
		g2<-g
		#print(Gpool)
		g2[searched,]<-Gpool[poolid,,na]
		CLRAM[[na]]<-CalcLikeRatioPerMarker(p,g,g2,searched,As,Ps,hG,sepGraphs,s2,Express)

	}
	CLRAM
}

• 出力関数

PrintStringCLRPM<-function(CLR){
	LocusName<-c("D8S1179","D21S11","D7S820","CSF1PO","D3S1358","TH01",
	"D13S317","D16S539","D2S1338","D19S433","VWA","TPOX","D18S51","D5S818","FGA")
	st<-paste("Locus","LR","Expression",sep="\t")
	st2<-paste("Locus","LR","Expression",sep=" ")
	print(st2)
	cprod<-1
	for(i in 1:length(CLR)){
		tmp<-paste(LocusName[i],CLR[[i]]$Pratio,sympy(CLR[[i]]$expression),sep="\t")
		tmp2<-paste(LocusName[i],CLR[[i]]$Pratio,sympy(CLR[[i]]$expression),sep=" ")
		st<-paste(st,tmp,sep="\n")
		#print(st)
		#st<-st+tmp+"\n"
		print(tmp2)
		#st<-paste(st,tmp,collapse="\n")
		cprod<-cprod*CLR[[i]]$Pratio
	}
	st2<-paste("cumulative LR",cprod,sep="\t")
	st<-paste(st,st2,sep="\n")
	#print(st)
	print(st2)
	st
}

• 実行コマンド
  • 式情報とかをため込むので、たくさんの家系で回すと重い
  • 適宜、回し方は調整を要す

searched<-3
st<-""
for(fi in 1:6){
	p<-pedigrees[[fi]]

	gs<-genotypesFamily[[fi]]


	CLout<-CalcLikeRatioAllMarker(p,gs,searched,fi,Gpool[,,],Alleles,Probs,Express=TRUE)
	st<-paste(st,PrintStringCLRPM(CLout),sep="\n")
}

st