ジェノタイプ推定で学ぶベイジアンネットワーク - ryamadaの遺伝学・遺伝統計学メモ

ベイジアンネットワークは確率事象の関係をグラフにして、情報が与えられたときの事後確率を計算してくれる
家系図があるとき、一部個人のジェノタイプがわかっているとき、ほかの個人のジェノタイプを推定する作業をベイジアンネットワーク化することができる
家系図を用いたジェノタイプ推定は、ノード間の関係が単純なので、ベイジアンネットワーク操作を考える上で、わかりやすい点も多いことから、これを題材にRのgRainパッケージの使用にトライする
kindle版
作者: ryamada
発売日: 2014/11/13
メディア: Kindle版
この商品を含むブログ (1件) を見る
以下のRmdファイル(Rmdファイルからhtml形式ファイルを作れば無料です。作り方はこちら)
注意。gRainパッケージのsetEvidence()関数の引数 nslistが廃止されたので、それに対応して、Rmdファイルは修正してあります。epubファイルは修正していません。
---
title: "STR型推定で学ぶベイジアン・ネットワーク"
author: "ryamada"
date: "Friday, November 14, 2014"
output: html_document
---

## 家系データの基本
DNA型推定における家系データの最小限の情報は、

* 構成員全員のIDベクトル
* 構成員の父親IDベクトル
* 構成員の母親IDベクトル
* 性別ベクトル
* 試料提供者を表すベクトル
* DNA型推定対象者を表すベクトル

の４ベクトルです。

構成員IDに通し番号をつけ、父母が構成員に含まれないときのIDを0とし、男女を1,2(不明を0)すれば、次のように家系図が描ける。

今、DNA試料が得られる個人を黒塗、遺伝子型を推定したい対象を赤塗にすれば、

```{r}
library(kinship2)
id <- 1:15
sex <- c(2,1,1,2,1,2,1,2,2,1,1,1,1,2,1)
dadid <- c(0,0,0,2,2,0,0,3,3,5,5,7,7,10,10)
momid <- c(0,0,0,1,1,0,0,4,4,6,6,8,8,9,9)
# 試料提供者と推定対象者を1にする
offer <- rep(0,length(id))
offer[c(5,6,10,11,12,13,14,15)] <- 1
# 推定対象者を赤指定する
target <- rep(1,length(id))
target[c(5,6,10,11)] <- 2
my.ped <- pedigree(id,dadid,momid,sex,offer)
# 家系図をプロットする
plot(my.ped,col=target)
```

## 遺伝子型推定のためのベイジアン・ネットワークの基礎

ベイジアン・ネットワークでは、有向グラフを作る。特に、ぐるぐる回り続けることが無いような有向グラフを作るが、それをDirected acyclic graph (DAG)と呼ぶ。
グラフなので、ノードとエッジとから成る。

遺伝子型推定のためのベイジアンネットワークでは、型を持つハプロタイプと、実験観察されるジェノタイプとをノードとする。

個人は、一つのジェノタイプを持ち、また、二つのハプロタイプを持つ。
個人の二つのハプロタイプは、片方が母親由来、もう片方が父親由来である。

ジェノタイプは、常染色体の場合、２つのハプロタイプのペアとして決まる。

ある個人の母由来ハプロタイプは、母の二つのハプロタイプのどちらか片方を確率0.5で受け継ぐ。

したがって、ハプロタイプを表すノードには、二つのハプロタイプのノードからエッジが入り、ジェノタイプのノードには、二つのハプロタイプのノードからエッジが入る。

ベイジアン・ネットワークでは、個々のノードと、それにエッジを入れるノードとのセットを１単位として考える。
ノードが入ってくるノードを子ノードと呼ぶことにし、エッジを入れるノードを親ノードと呼ぶことにすれば、このセットには、必ず１個の子ノードがあり、１個以上の親ノードがある。
子ノードがとる値(ハプロタイプノードの場合はアレル、ジェノタイプノードの場合は、ディプロイドジェノタイプ)は、親ノードがとる値によって確率的に決まるが、それはテーブルで表すことができる。

## Rの準備

Rでベイジアン・ネットワークを使うために、gRainパッケージを使う。
この"gRain"の使用には、あらかじめ、"graph","RBGL","Rgraphviz"の３パッケージをバイオコンダクターサイトから取得し読み込んでおく必要がある。
```{}
source("http://bioconductor.org/biocLite.R")
biocLite(c("graph","RBGL","Rgraphviz"))
```
```{r}
library(RBGL)
```

そのうえで、gRainパッケージとgRbaseパッケージとをCRANから取得して読み込む

```{}
install.packages(c("gRain","gRbase"))
```
```{r}
library(gRbase)
library(gRain)
```
## 一人のジェノタイプとハプロタイプのみのネットワーク

ごく簡単な例でベイジアン・ネットワークの使い方を学ぶこととする。

ジェノタイプノード g と父母由来のハプロタイプノード hp,hmとから成る。

アレル数naを4とし、アレル名を、A=("a1",...,"ana")とする。
またアレル頻度をP=(p1,p2,p3,p4)とする。

二つのアレルが作るディプロタイプGは"a1=a1","a2=a2",...,"ana=ana"というホモタイプと、
"a1=a2","a1=a3", ... "a(na-1)=ana"というヘテロタイプとがある。
ホモタイプの種類数はna、ヘテロタイプの種類数はna(na-1)/2を併せた、ng = na+na(na-1)/2=(na+1)na/2である。

```{r}
na <- 4
P <- c(0.1,0.2,0.3,0.4)
A <- c("a1","a2","a3","a4")
G <- c("a1=a1","a2=a2","a3=a3","a4=a4","a1=a2","a1=a3","a1=a4","a2=a3","a2=a4","a3=a4")
```
次にgRainパッケージを用いて、ネットワークを作る。

二つのハプロタイプとその子にあたるジェノタイプとの関係は、hpがaiでありhmがajであるときには、
aiとajが同じならば、gは"ai=ai"になり、そうでなければ"ai=aj" (ただしi < j)になるわけだが、それを

$ng \times na \times na$のアレイに収める。

```{r}
hh2g <- function(n){
  N <- (n)*(n+1)/2
  ret <- array(0,c(n,n,n,n))
	for(i in 1:n){
		for(j in 1:n){
			ret[i,j,i,j] <- ret[i,j,i,j] + 0.5
			ret[i,j,j,i] <- ret[i,j,j,i] + 0.5
		}
	}
	ret2 <- array(0,c(N,n,n))
	for(i in 1:n){
		for(j in 1:n){
			tmp <- ret[i,j,,]
			tmp2 <- tmp + t(tmp)
			tmp3 <- diag(tmp)
			for(k in 1:(n-1)){
				tmp3 <- c(tmp3,tmp2[k,(k+1):n])
			}
			ret2[,i,j] <- tmp3
		}
	}
	ret2
}
HH2G <- hh2g(na)
```

実際、a1とa1からは第１のディプロタイプである"a1=a1"が確率1で生じ、
a1とa2からは、"a1=a2"なるヘテロ型が確率1で生じるので、この３次元アレイは、要素が0か1である。
```{r}
HH2G[,1,1]
HH2G[,1,2]
HH2G[,3,3]
```
二つのハプロタイプのノードは、名称"hp","hm"を指定しつつ、とりうる値(levels)はAとして与え、その確率(values)は、頻度ベクトルPで与えている。
```{r}
h.p <- cptable(~hp,values=P,levels=A)
h.m <- cptable(~hm,values=P,levels=A)
h.p
h.m
```
ジェノタイプのノードは、確率が、親ノードに依存しているので、親ノードを指定しつつ、うえで作った、アレイを与えて作成する。
```{r}
g <- cptable(~g : hp + hm,values=HH2G,levels=G)
g
```

作成したベイジアン・ネットワークノードを表示させると

{v,pa(v)}というものが現れるが、これは、自身の名称と、親ノード(pa(v))の名称である。ハプロタイプノードは親ノードがないので、自身の名称のみが記載され、ジェノタイプノードの場合は、親ノード名が現れている。

次項のlevels of vはとりうる値であり、
valuesは、親のないノードの場合は、その確率が、親ノードがある場合には、アレイのセルの値が記載される。

### ノードをネットワークにする

ノードの作成が終わったら、それをネットワークにまとめる。
```{r}
# ノードをリストオブジェクトにして
tmplist <- list(h.p,h.m,g)
# それをネットワークにするには２つの関数compileCPT(),grain()を使うので、それを関数にまとめて
my.grain <- function(L){
  tmp <- compileCPT(L)
  grain(tmp)
}
# ベイジアン・ネットワークオブジェクトを作る
bn <- my.grain(tmplist)
summary(bn)
# プロット
plot(bn)
# ノードの値の確率を問い合わせる
querygrain(bn)
```

描かれるネットワーク図は、期待したものになっており、各ノードの値の確率も、期待した通りである。

ここで、ジェノタイプを観察して、それが"a1=a2"だったとする。
```{r}
bn.1 <- setEvidence(bn, evidence=list(g="a1=a2"))
plot(bn.1)
summary(bn.1)
```

ベイジアン・ネットワークでの計算にあたって、グラフを閉じて(部分的に完全グラフ(クリーク)を作る)確率伝播法を適用する方が効率がよいので、それを適用し、グラフが三角形に閉じている。
このクリークとは
グラフの中の一部であって、あるノードのセットを見ると、すべてのノードペアにエッジがあるもののことである。そのクリーク(clique)数がサマリーに書き出されている。また、クリークの中で最大のもの(ここでは三角形)のノード数も記されている。

さらにMaximum state space in cliques が160というのは、hp,hmのそれぞれの場合が4,4、gの場合が10あるわけだが、それの組み合わせの場合の数$4 \times 4 \times 10 = 160$ のことである。

ジェノタイプの観察結果を与えた後で、各ノードの値の確率を問い合わせると(gはわかっているので、出ないが)、返ってくる。
```{r}
querygrain(bn.1,result="data.frame")
```

また、hp,hm,gの3つの組み合わせではどうなっているかを問い合わせることもできる。

```{r}
querygrain(bn.1,type="joint",result="data.frame")
```
わかりきっているgの情報を出さないようにしてやると、
```{r}
querygrain(bn.1,type="joint",nodes = c("hp","hm"),result="data.frame")
```
少し簡単になる。

これを見ると、(hp,hm)=("a1","a2")の確率が0.5で、(hp,hm)=("a2","a1")の確率が0.5であり、それですべてであることがわかる。

組み合わせで問い合わせないと、hpとhmとが、それぞれ、"a1","a2"の確率が0.5であることがわかったが、それだけでは、hp,hmの両方が"a1"である可能性もあることはわからない。
このように、細かい組み合わせの情報は出さずに、合算して出すことを「Marginal(合算)」で出す、と言う。
練習として、異なるジェノタイプを指定したり、ハプロタイプの値を指定したりして、問い合わせの結果がどう変わるかを見ておく。

```{r}
#querygrain(setEvidence(bn,nslist=list(hp="a3")),result="data.frame")
#querygrain(setEvidence(bn,nslist=list(hp="a3")),type="joint",nodes = c("hm","g"),result="data.frame")
#querygrain(setEvidence(bn,nslist=list(hp="a3",hm="a2")),result="data.frame")
querygrain(setEvidence(bn,evidence=list(hp="a3")),result="data.frame")
querygrain(setEvidence(bn,evidence=list(hp="a3")),type="joint",nodes = c("hm","g"),result="data.frame")
querygrain(setEvidence(bn,evidence=list(hp="a3",hm="a2")),result="data.frame")
```

## もう一つの関係、ハプロタイプの親から子への伝搬

ある個人が２つのハプロタイプhp.1,hm.1を持っているときに、それを子供に伝えるとき、父親として伝えるなら、子供の父由来ハプロタイプhp.2へとエッジが入る(母親として伝えるときは子供の母由来ハプロタイプhm.2へとエッジが入る)。

どのアレルが伝わるかは0.5の確率なので、以下のようにしてネットワークを作成する。

まずは、伝達具合のアレイを作成する。このアレイは、親の二つのハプロタイプと、子供の一つのハプロタイプとの関係なので、na x na x na のアレイとなる。
```{r}
hh2h <- function(n){
  ret <- array(0,rep(n,3))
	for(i in 1:n){
		for(j in 1:n){
			ret[i,i,j] <- ret[i,i,j] + 0.5
			ret[j,i,j] <- ret[j,i,j] + 0.5
		}
	}
	ret
}
HH2H <- hh2h(na)
```
このアレイの中身を少し見てみる。
```{r}
HH2H[,1,1] # 親の二つのハプロタイプのアレルが、ともに"a1"の場合
HH2H[,2,3] # 親の二つのハプロタイプのアレルが"a2","a3"の場合
```
```{r}
h.p1 <- cptable(~hp1,values=P,levels=A)
h.m1 <- cptable(~hm1,values=P,levels=A)
h.p1
h.m1
```

このアレイを使って、子供のハプロタイプノードを作成する。
```{r}
h.p2 <- cptable(~hp2 : hp1 + hm1,values=HH2H,levels=A)
h.p2
```

この先の手順はジェノタイプの場合と同じで、ネットワーク化して、情報を与えて、問い合わせる、となる。
```{r}
# ベイジアン・ネットワークオブジェクトを作る
tmplist <- list(h.p1,h.m1,h.p2)
bn2 <- my.grain(tmplist)
summary(bn2)
# プロット
plot(bn2)
# ノードの値の確率を問い合わせる
querygrain(bn2)
```

描かれるネットワーク図は、期待したものになっており、各ノードの値の確率も、期待した通りである。

ここで、子の父由来ハプロタイプを観察して、それが"a1"だったとする。
```{r}
bn2.1 <- setEvidence(bn2, evidence=list(hp2="a1"))
plot(bn2.1)
summary(bn2.1)
querygrain(bn2.1,result="data.frame",type="marginal")
querygrain(bn2.1,result="data.frame",type="joint",nodes =c("hp1","hm1"))
```

## 親子関係を作る

ここまでで、二つのことをネットワーク化した。

*個人のハプロタイプノード二つとジェノタイプノード一つの組み合わせ
*親から子へのハプロタイプの伝達

の二つである。

今度はこれを組み合わせ、両親と一人の子供の三人をネットワーク化する。

各個人は、ハプロタイプノードを二つ、ジェノタイプノードを一つ持つので、全部で3 x 3 = 9ノードで構成される。

子のジェノタイプを観察して、父親と母親のジェノタイプを推定することにする。

家系図を描いておく。
```{r}
id <- 1:3
sex <- c(1,2,1)
dadid <- c(0,0,1)
momid <- c(0,0,2)
# 試料提供者と推定対象者を1にする
offer <- rep(0,length(id))
offer[c(1,2,3)] <- 1
# 推定対象者を赤指定する
target <- rep(1,length(id))
target[c(1,2)] <- 2
my.ped <- pedigree(id,dadid,momid,sex,offer)
# 家系図をプロットする
plot(my.ped,col=target)
```

ノードの数が増えてくると、手作業は大変なので、家系図に与える情報(個人IDリスト、父親IDリスト、母親IDリスト)と、アレル名、その頻度を与えて、ネットワーク化する関数を作っておく。


```{r}
# id 個人ID ベクトル
# dadid 父親ID ベクトル
# momid 母親ID ベクトル
# haplo.lab アレル名ベクトル
# freq.hap アレル頻度ベクトル
make.genotype.bn <- function(id,dadid,momid,haplo.lab,freq.hap){
  # 人数
  n.ind <- length(id)
  # アレル数
	n.allele <- length(haplo.lab)
  # ジェノタイプ名を作成
	gen.lab <- c()
	for(i in 1:n.allele){
		gen.lab <- c(gen.lab,paste(haplo.lab[i],"=",haplo.lab[i],sep=""))
	}
	for(i in 1:(n.allele-1)){
		for(j in (i+1):n.allele){
			gen.lab <- c(gen.lab,paste(haplo.lab[i],"=",haplo.lab[j],sep=""))
		}
	}
  # ノード連結アレイの作成
	outh <- hh2h(n.allele)
	outg <- hh2g(n.allele)
	
	# freq.gen　ジェノタイプの頻度ベクトル(HWE仮定)
	freq.gen <- freq.hap^2
	for(i in 1:(n.allele-1)){
		for(j in (i+1):n.allele){
			freq.gen <- c(freq.gen,freq.hap[i]*freq.hap[j]*2)
		}
	}
	
	# 父型ハプロタイプ
	h.p <- paste("h.p",1:n.ind,sep="")
	# 母型ハプロタイプ
	h.m <- paste("h.m",1:n.ind,sep="")
	# ジェノタイプ
	g <- paste("g",1:n.ind,sep="")
	
	# ハプロタイプ・ジェノタイプのtablesを仮に作る
	hps <- hms <- gs <- list()
	
	for(i in 1:n.ind){
		hps[[i]] <- cptable(~hp,values=freq.hap,levels=haplo.lab)
		hms[[i]] <- cptable(~hm,values=freq.hap,levels=haplo.lab)
		gs[[i]] <- cptable(~g,values=freq.gen,levels=gen.lab)
	}
  # 伝達関係、個人内ハプロタイプ・ジェノタイプ関係を登録する
	for(i in 1:n.ind){
		if(dadid[i]!=0){
			hps[[i]]$vpa <- c(paste("hp",i,sep=""),paste("hp",dadid[i],sep=""),paste("hm",dadid[i],sep=""))
			hps[[i]]$values=outh
		}else{
			hps[[i]]$vpa[1] <- paste("hp",i,sep="")
		}
		if(momid[i]!=0){
			hms[[i]]$vpa <- c(paste("hm",i,sep=""), paste("hp",momid[i],sep=""), paste("hm",momid[i],sep=""))
			hms[[i]]$values <- outh
		}else{
			hms[[i]]$vpa[1] <- paste("hm",i,sep="")
		}
		gs[[i]]$vpa <- c(paste("g",i,sep=""), paste("hp",i,sep=""), paste("hm",i,sep=""))
		gs[[i]]$values <- outg
	}
  # すべてをリストにまとめる
	tmplist <- list()
	cnt <- 1
	for(i in 1:n.ind){
		tmplist[[cnt]] <- hps[[i]]
		cnt <- cnt+1
		tmplist[[cnt]] <- hms[[i]]
		cnt <- cnt+1
		tmplist[[cnt]] <- gs[[i]]
		cnt <- cnt+1
	}
  # ネットワークオブジェクトにして返す
	my.grain(tmplist)
}
```

```{r}
bn.trio <- make.genotype.bn(id,dadid,momid,A,P)
plot(bn.trio)
```
確かに、九つのノードが作成され、その関係も適切であることが見て取れる。

```{r}
querygrain(bn.trio,result="data.frame",type="marginal")
```

子のジェノタイプが"a1=a1"とする。

```{r}
bn.trio1 <- setEvidence(bn.trio,evidence=list(g3="a1=a1"))
querygrain(bn.trio1,result="data.frame",type="marginal")
```

父親と母親のジェノタイプの組み合わせで問い合わせることにする。

場合の数が大きくなって計算が重すぎるので、少し工夫をする。
まず、ネットワークに問い合わせる前に、あらかじめ、「g1,g2の組み合わせに興味がある」ということを、以下のように指定して、新たなネットワークオブジェクトを作成する。

```{r}
bn.trio2 <- compile(bn.trio,root=c("g1","g2"),propagate=TRUE)
```

その上で、情報を与え、興味のある組み合わせについて、Joint オプションで問い合わせる。
```{r}
bn.trio3 <- setEvidence(bn.trio2,evidence=list(g3="a1=a1"))
querygrain(bn.trio3,result="data.frame",type="marginal")
querygrain(bn.trio3,nodes=c("g1","g2"),result="data.frame",type="joint")
```

## 本番

冒頭で考えた家系図について実行する。

```{r}
id <- 1:15
sex <- c(2,1,1,2,1,2,1,2,2,1,1,1,1,2,1)
dadid <- c(0,0,0,2,2,0,0,3,3,5,5,7,7,10,10)
momid <- c(0,0,0,1,1,0,0,4,4,6,6,8,8,9,9)
# 試料提供者と推定対象者を1にする
offer <- rep(0,length(id))
offer[c(5,6,10,11,12,13,14,15)] <- 1
# 推定対象者を赤指定する
target <- rep(1,length(id))
target[c(5,6,10,11)] <- 2
my.ped <- pedigree(id,dadid,momid,sex,offer)
# 家系図をプロットする
plot(my.ped,col=target)
```


```{r}
bn.pedigree <- make.genotype.bn(id,dadid,momid,A,P)
plot(bn.pedigree)
```

推定対象者のジェノタイプを指定してネットワークを作りかえる。

```{r}
targets <- paste("g",c(5,6,10,11),sep="")
bn.pedigree2 <- compile(bn.pedigree,root=targets,propagate=TRUE)
```

```{r}
out2m <- querygrain(bn.pedigree2,nodes = targets,result="data.frame",type="marginal")
out2m
out2j <- querygrain(bn.pedigree2,nodes = targets,result="data.frame",type="joint")
# 場合の数が多すぎるので、場合の数を表示する
print(length(out2j$Freq))
# 確率が0ではない場合の数は
print(length(out2j$Freq[which(out2j$Freq!=0)]))
```

12,13,14,15のジェノタイプが、すべて"a1=a1"だとする。

```{r}
bn.pedigree3 <- setEvidence(bn.pedigree2,evidence=list(g12="a1=a1",g13="a1=a1",g14="a1=a1",g15="a1=a1"))
out3m <- querygrain(bn.pedigree3,nodes = targets,result="data.frame",type="marginal")
out3m
out3j <- querygrain(bn.pedigree3,nodes = targets,result="data.frame",type="joint")
# 確率が0ではない場合の数は
print(length(out3j$Freq[which(out3j$Freq!=0)]))
out3j[which(out3j$Freq!=0),]
```