数学・コンピュータ関連の姉妹ブログ『ryamadaのコンピュータ・数学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

核家族のルール〜２倍体生殖は面倒くさい〜

R メンデルの法則
  • ここ数日、Rのsetsパッケージの勉強をしている
  • それを使って核家族のルールを集合演算化できないかといじっている
  • いわゆる集合だけでなく、ファジィマルチ集合も投入してみた
  • どれをどう使ってもかなり面倒くさい(要するに、集合演算としてさらさらとコードを書くことを断念したということ)
  • どこが面倒くさいのかについて考えてみる
  • できたコードは、集合的な考え方をしつつ、ベクトルで回すコードになったわけだが、考えてみるに当たり、その『でき』から逆算的に面倒くささについて書いてみる
    • 家系図にある要素について
      • 個人(を表すノード)がある
      • 個人間のつながり(エッジ)がある
    • 個人ノードは「ディプロタイプ」という情報を持つ(持たずに推定されることもある)
      • この情報は「要素数１もしくは２の集合」である。実験観察のみに基づけば、「順序なし」である
    • 個人間エッジは「親子関係」を表していると同時に、そのエッジには「伝達アレル」が紐づいている。
      • 「親子関係」は個人が持つディプロタイプという集合をつなぐ
      • 「伝達アレル」に着目すれば、エッジには、要素数１の集合が属している
      • この「伝達アレル」を個人の属性として考えることもできて、その場合は、個人の持つ２つのアレルに由来親の区別をつけることに相当する
    • 親子トリオ
      • 親子トリオのそれぞれの「ディプロタイプ集合」は、必ず交わりのある集合である
      • トリオの「ディプロタイプ集合」はすべて重なることもあれば、そうならないときもあるが、子の集合はかならず親の集合と１個以上の共通要素を持つこと、３集合の和集合の要素数の最大値は４(最小値は１)である
      • 「伝達アレル」という要素数１の集合は、この「ディプロタイプ３集合」のベン図上に割り当てることができる
      • 要素数が限定されている(1，2)ことから、「相対的補集合」の有無が限定的であるために、３人のうちの２人のディプロタイプ集合情報、伝達情報から、残りの１人のディプロタイプ集合情報、その１人を含む伝達アレル集合情報が、確定的に決まりうるなどの特殊性が生じる
    • 不明ディプロタイプ、不明伝達アレルがある場合には、ある程度、限定できるのは、要素数が1,2に限定された集合が狭い空間(総要素数はアレル数の上限がある)にひしめいているから
    • このように「狭い空間」であること、「相対的補集合」がとりやすいことから、推定作業にあたっては、「各集合に絶対に属する要素」「各集合に属することがありえない要素」「各集合に属するか属さないか確定できない要素」の３状態に分けてわりふりすることが便利になる
    • 「３状態」での割り振りは、「普通の集合」が２状態に関する処理のためにできていることから、「普通の集合」として取扱いにくくしているようだ
    • 「３状態」であるから、「ファジィ集合」を使うとうまくいきそうだが、要素数が少ないことを用いた「確定作業」をするときには「０，１」っぽい扱いが便利であるために、やりにくさが生じ、また、要素数が２の集合と要素数が１の集合とが入り乱れることも、「ファジィ集合」を導入しさえすれば楽になる、というわけではないことの原因のようだ
  • このように苦労してみると、２倍体というシステムが、本当によくできていると思える。非常に少ない数で非常に単純なルールしかないけれど、データのハンドリング自体は面倒くさい、というのが、「よくできている」という意味である
  • 「よくできている」ことに感心すればするほど、その「よくできた仕組み」を逸脱する介入にはネガティブな立場になってしまう
# 家系を作る
# ジェノタイプ不明者(第５カラムが２もしくは３)も入れる
p<-matrix(
c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,15,16,17,18,19,
  0, 0, 0, 0, 2, 2, 4, 4, 6,  6,  0,  0, 12, 13,13,13,13,13,13,
  0, 0, 0, 0, 1, 1, 3, 3, 7,  7,  0,  0, 11, 10,10,10,10,10,10,
  0, 1, 0, 1, 0, 1, 0, 1, 0,  0,  0,  1,  1,  1,1,1,1,1,1,
  3,1,1,1,3,1,1,1,1,2,1,1,2,1,1,1,1,1,1),
  ncol=5)
# アレル名と頻度
A<-c( "20", "21", "22", "23")
P<-runif(length(A))
P<-P/sum(P) 
# メンデルの法則を満足するディプロタイプを代入する
g<-matrix(c("0", "0",
 "20", "21",
 "22", "21",
 "20", "21",
 "0" , "0" ,
 "21", "21",
 "20", "21",
 "20", "21",
 "21", "21",
 "0" , "0" ,
 "21", "20",
 "23", "23",
 "0" , "0" ,
 "23", "21",
 "21", "21",
 "23", "21",
 "21", "21",
 "23", "21",
 "23", "21"),
 ncol=2,byrow=TRUE)
  • これから、ディプロタイプと伝達アレルに関して「あり得る場合」を抜き出してみる
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]])
	}
	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)
print("diplotype")
print(LDZout$D)
print("maternal allele")
print(LDZout$H[,1,])
print(LDZout$D)
print("paternal allele")
print(LDZout$H[,2,])

> print(LDZout$D)
      [,1] [,2] [,3] [,4] [,5]
 [1,]    0    0    0    0    0
 [2,]    1   -1   -1    1   -1
 [3,]    1    1   -1   -1   -1
 [4,]   -1   -1   -1    1   -1
 [5,]    0    0    0    0    0
 [6,]    1   -1   -1    1   -1
 [7,]    1   -1   -1    1   -1
 [8,]   -1    1   -1    1   -1
 [9,]    1   -1   -1    1   -1
[10,]    1   -1   -1    1   -1
[11,]   -1    1   -1   -1   -1
[12,]   -1    1   -1    1   -1
[13,]   -1    1   -1    1   -1
[14,]   -1    1   -1    1   -1
[15,]   -1    1   -1    1   -1
[16,]   -1    1   -1    1   -1
[17,]   -1    1   -1    1   -1
[18,]   -1    1   -1    1   -1
[19,]   -1    1   -1    1   -1
> print("maternal allele")
[1] "maternal allele"
> print(LDZout$H[,1,])
      [,1] [,2] [,3] [,4] [,5]
 [1,]    0    0    0    0    0
 [2,]    0   -1   -1    0   -1
 [3,]    0    0   -1   -1   -1
 [4,]   -1   -1   -1    0   -1
 [5,]    0   -1   -1    0   -1
 [6,]    0   -1   -1    0   -1
 [7,]   -1   -1   -1    1   -1
 [8,]   -1   -1   -1    1   -1
 [9,]    0   -1   -1    0   -1
[10,]    0   -1   -1    0   -1
[11,]   -1    0   -1   -1   -1
[12,]   -1    0   -1    0   -1
[13,]   -1   -1   -1    1   -1
[14,]   -1    1   -1   -1   -1
[15,]   -1    1   -1   -1   -1
[16,]   -1    1   -1   -1   -1
[17,]   -1    1   -1   -1   -1
[18,]   -1    1   -1   -1   -1
[19,]   -1    1   -1   -1   -1
> print(LDZout$D)
      [,1] [,2] [,3] [,4] [,5]
 [1,]    0    0    0    0    0
 [2,]    1   -1   -1    1   -1
 [3,]    1    1   -1   -1   -1
 [4,]   -1   -1   -1    1   -1
 [5,]    0    0    0    0    0
 [6,]    1   -1   -1    1   -1
 [7,]    1   -1   -1    1   -1
 [8,]   -1    1   -1    1   -1
 [9,]    1   -1   -1    1   -1
[10,]    1   -1   -1    1   -1
[11,]   -1    1   -1   -1   -1
[12,]   -1    1   -1    1   -1
[13,]   -1    1   -1    1   -1
[14,]   -1    1   -1    1   -1
[15,]   -1    1   -1    1   -1
[16,]   -1    1   -1    1   -1
[17,]   -1    1   -1    1   -1
[18,]   -1    1   -1    1   -1
[19,]   -1    1   -1    1   -1
> print("paternal allele")
[1] "paternal allele"
> print(LDZout$H[,2,])
      [,1] [,2] [,3] [,4] [,5]
 [1,]    0    0    0    0    0
 [2,]    0   -1   -1    0   -1
 [3,]    0    0   -1   -1   -1
 [4,]   -1   -1   -1    0   -1
 [5,]    0    0    0    0    0
 [6,]    0   -1   -1    0   -1
 [7,]    1   -1   -1   -1   -1
 [8,]   -1    1   -1   -1   -1
 [9,]    0   -1   -1    0   -1
[10,]    0   -1   -1    0   -1
[11,]   -1    0   -1   -1   -1
[12,]   -1    0   -1    0   -1
[13,]   -1    1   -1   -1   -1
[14,]   -1   -1   -1    1   -1
[15,]   -1   -1   -1    1   -1
[16,]   -1   -1   -1    1   -1
[17,]   -1   -1   -1    1   -1
[18,]   -1   -1   -1    1   -1
[19,]   -1   -1   -1    1   -1
  • あり得るディプロタイプと伝達アレルを書き出してみよう
LDZout2<-LimitDiplotypeZ(p,g2,A)
printLDZout<-function(ldz,A){
	for(i in 1:length(ldz$D[,1])){
		dset<-as.set(which(ldz$D[i,]>=0))
		hset1<-as.set(which(ldz$H[i,1,]>=0))
		hset2<-as.set(which(ldz$H[i,2,]>=0))
		print("---")
		print(i)
		#print(dset)
		print(A[unlist(dset)])
		print(as.set(hset1*hset2))
		print(A[unlist(as.set(hset1*hset2))])

		print(hset1)
		print(A[unlist(hset1)])
		print(hset2)
		print(A[unlist(hset2)])

	}
}
printLDZout(LDZout,A)

> printLDZout(LDZout,A)
[1] "---"
[1] 1
[1] "14" "15" "16" "17" "18"
{(1L, 1L), (1L, 2L), (1L, 3L), (1L, 4L), (1L, 5L), (2L, 1L), (2L, 2L), (2L,
 3L), (2L, 4L), (2L, 5L), (3L, 1L), (3L, 2L), (3L, 3L), (3L, 4L), (3L, 5L),
 (4L, 1L), (4L, 2L), (4L, 3L), (4L, 4L), (4L, 5L), (5L, 1L), (5L, 2L), (5L,
 3L), (5L, 4L), (5L, 5L)}
 [1] "14" "14" "14" "15" "14" "16" "14" "17" "14" "18" "15" "14" "15" "15" "15" "16"
[17] "15" "17" "15" "18" "16" "14" "16" "15" "16" "16" "16" "17" "16" "18" "17" "14"
[33] "17" "15" "17" "16" "17" "17" "17" "18" "18" "14" "18" "15" "18" "16" "18" "17"
[49] "18" "18"
{1L, 2L, 3L, 4L, 5L}
[1] "14" "15" "16" "17" "18"
{1L, 2L, 3L, 4L, 5L}
[1] "14" "15" "16" "17" "18"
[1] "---"
[1] 2
[1] "14" "17"
{(1L, 1L), (1L, 4L), (4L, 1L), (4L, 4L)}
[1] "14" "14" "14" "17" "17" "14" "17" "17"
{1L, 4L}
[1] "14" "17"
{1L, 4L}
[1] "14" "17"
[1] "---"
[1] 3
[1] "14" "15"
{(1L, 1L), (1L, 2L), (2L, 1L), (2L, 2L)}
[1] "14" "14" "14" "15" "15" "14" "15" "15"
{1L, 2L}
[1] "14" "15"
{1L, 2L}
[1] "14" "15"
[1] "---"
[1] 4
[1] "17"
{(4L, 4L)}
[1] "17" "17"
{4L}
[1] "17"
{4L}
[1] "17"
[1] "---"
[1] 5
[1] "14" "15" "16" "17" "18"
{(1L, 1L), (1L, 2L), (1L, 3L), (1L, 4L), (1L, 5L), (4L, 1L), (4L, 2L), (4L,
 3L), (4L, 4L), (4L, 5L)}
 [1] "14" "14" "14" "15" "14" "16" "14" "17" "14" "18" "17" "14" "17" "15" "17" "16"
[17] "17" "17" "17" "18"
{1L, 4L}
[1] "14" "17"
{1L, 2L, 3L, 4L, 5L}
[1] "14" "15" "16" "17" "18"
[1] "---"
[1] 6
[1] "14" "17"
{(1L, 1L), (1L, 4L), (4L, 1L), (4L, 4L)}
[1] "14" "14" "14" "17" "17" "14" "17" "17"
{1L, 4L}
[1] "14" "17"
{1L, 4L}
[1] "14" "17"
[1] "---"
[1] 7
[1] "14" "15" "17"
{(4L, 1L), (4L, 2L)}
[1] "17" "14" "17" "15"
{4L}
[1] "17"
{1L, 2L}
[1] "14" "15"
[1] "---"
[1] 8
[1] "15" "17"
{(4L, 2L)}
[1] "17" "15"
{4L}
[1] "17"
{2L}
[1] "15"
[1] "---"
[1] 9
[1] "14" "15" "17"
{(1L, 1L), (1L, 2L), (1L, 4L), (4L, 1L), (4L, 2L), (4L, 4L)}
 [1] "14" "14" "14" "15" "14" "17" "17" "14" "17" "15" "17" "17"
{1L, 4L}
[1] "14" "17"
{1L, 2L, 4L}
[1] "14" "15" "17"
[1] "---"
[1] 10
[1] "14" "15" "17"
{(1L, 1L), (1L, 2L), (1L, 4L), (4L, 1L), (4L, 2L), (4L, 4L)}
 [1] "14" "14" "14" "15" "14" "17" "17" "14" "17" "15" "17" "17"
{1L, 4L}
[1] "14" "17"
{1L, 2L, 4L}
[1] "14" "15" "17"
[1] "---"
[1] 11
[1] "15"
{(2L, 2L)}
[1] "15" "15"
{2L}
[1] "15"
{2L}
[1] "15"
[1] "---"
[1] 12
[1] "15" "17"
{(2L, 2L), (2L, 4L), (4L, 2L), (4L, 4L)}
[1] "15" "15" "15" "17" "17" "15" "17" "17"
{2L, 4L}
[1] "15" "17"
{2L, 4L}
[1] "15" "17"
[1] "---"
[1] 13
[1] "15" "17"
{(2L, 2L), (4L, 2L)}
[1] "15" "15" "17" "15"
{2L, 4L}
[1] "15" "17"
{2L}
[1] "15"
[1] "---"
[1] 14
[1] "15" "17"
{(2L, 2L), (2L, 4L), (4L, 2L), (4L, 4L)}
[1] "15" "15" "15" "17" "17" "15" "17" "17"
{2L, 4L}
[1] "15" "17"
{2L, 4L}
[1] "15" "17"
[1] "---"
[1] 15
[1] "15" "17"
{(2L, 2L), (2L, 4L), (4L, 2L), (4L, 4L)}
[1] "15" "15" "15" "17" "17" "15" "17" "17"
{2L, 4L}
[1] "15" "17"
{2L, 4L}
[1] "15" "17"
[1] "---"
[1] 16
[1] "15" "17"
{(2L, 2L), (2L, 4L), (4L, 2L), (4L, 4L)}
[1] "15" "15" "15" "17" "17" "15" "17" "17"
{2L, 4L}
[1] "15" "17"
{2L, 4L}
[1] "15" "17"
[1] "---"
[1] 17
[1] "15" "17"
{(2L, 2L), (2L, 4L), (4L, 2L), (4L, 4L)}
[1] "15" "15" "15" "17" "17" "15" "17" "17"
{2L, 4L}
[1] "15" "17"
{2L, 4L}
[1] "15" "17"
[1] "---"
[1] 18
[1] "15" "17"
{(2L, 2L), (2L, 4L), (4L, 2L), (4L, 4L)}
[1] "15" "15" "15" "17" "17" "15" "17" "17"
{2L, 4L}
[1] "15" "17"
{2L, 4L}
[1] "15" "17"
[1] "---"
[1] 19
[1] "15" "17"
{(2L, 2L), (2L, 4L), (4L, 2L), (4L, 4L)}
[1] "15" "15" "15" "17" "17" "15" "17" "17"
{2L, 4L}
[1] "15" "17"
{2L, 4L}
[1] "15" "17"

logo

図書室

日記全体の目次
(記事一覧)

疾患遺伝子解析を学ぶためのカタログ(H18 3月改訂)の目次

 教科書

京都大学大学院医学研究科集中講義 全８限Part I(2007)
京都大学大学院医学研究科集中講義 全８限Part II(2008)
京都大学大学院医学研究科統計遺伝学講義 全８限(2009)

以下の関連サイトへのリンク道案内ページ
遺伝子多型解析Wiki
東大医科研HGCゲノム機能解析分野(旧)
理研リウマチチーム解析資料ページ


疾患遺伝子解析を学ぶためのカタログ(H18 3月改訂)
遺伝統計学のための統計学基礎 全９限
Haploviewを用いた連鎖不平衡マッピング解析実習 全９限
Arlequinを用いた集団遺伝学実習 全７限
一から始める遺伝子多型・形質関連解析
多数のGWAS/GWASメタアナリシスとその周辺解析用のアプリケーションへのリンク

複合遺伝性疾患とは

----------News----------
20万アクセス
2008/07/20
40万アクセス
2010/09/12


-------------------------