RmarkdownのなんちゃってKyotoPRO
- 昨日の記事で、BRACAPROとかその姉妹版がどんなものかを理解するためのイントロを書いた
- 今日の記事はその清書
- Rmarkdownでepub可能な形にしました
- この本の姉妹編(分離比解析)はこちら
- 詳しくは以下の通り
- Rmdファイルです。html化、epub化できます(やり方はこちら。やり方を覚えれば無償です)
- html化、epub化が面倒くさければ、kindleで2米ドルでも(こちらに現れる予定)
発病リスク推定のためのベイジアンネットワーク?オミックス統計学2014_6? 統計遺伝学
- 作者: ryamada
- 発売日: 2014/02/16
- メディア: Kindle版
- この商品を含むブログ (1件) を見る
- 2014年度京大医学研究科修士・博士向け講義「オミックス統計学入門」の前期7月くらいにやる予定(場合によっては後期にずれ込む可能性も大)
- 必要パッケージ関連情報
# ジェノタイプとフェノタイプとでベイジアンネットワーク 家系図があり、そこに、ある遺伝性フェノタイプが観察されているとする。 このフェノタイプは加齢とともに発病率が高くなるようなものであり、その加齢性上昇のパターンがジェノタイプによって決まるモデルを考える。 そのようなモデルにおいて、家系図上の個人について、登録年齢と発病年齢と登録時フェノタイプの情報が得られたり、ジェノタイプ情報が得られたりすると、それらの情報によって、情報が得られていない個人のジェノタイプが推測ができ、その推測フェノタイプにより発病年齢予測ができたりするだろう。 この情報活用をベイジアンネットワークを用いて行うこととする。 このベイジアンネットワークの活用を通じて、家系集積性フェノタイプについての基本も学ぶ。 ## 遺伝的形質のベイジアンネットワークの基礎 ### １人のベイジアンネットワーク １人について考える。 この人はジェノタイプとフェノタイプを持つ。 今、ジェノタイプはMM,Mm,mmの３通りあり、MMなら発病し、Mm,mmならば発病しないとする。 これは劣性遺伝形式である。 ジェノタイプを測定していないとき、MM,Mm,mmのどれなのかは不明ですが、確率的にP(MM),P(Mm),P(mm)と考えてもい良い。これをXgとする。 フェノタイプもわかっていないときには、病気である=1、病気ではない=0のどちらかになりますが、これをXpで表して、それぞれの確率をQ(1),Q(0)とする。 XpはXgによって決まるので、ベイジアンネットワークではXgからXpへ矢印(エッジ)を引く。 そして、Xpの確率はXgの確率によって決まるが、その関係は ||X=MM|X==Mm|X=mm| |:--:|:--:|:--:|:--:| |P=1|1|0|0| |P=0|0|1|1| となる。 今、P(MM),P(Mm),P(mm)がある集団で、0.01,0.18,0.81の割合であるとすれば、この個人のジェノタイプの確率もこの割合そのものであるとみなしてもよく、そうすると、この個人が発病者である条件付き確率は、0.01である。 これをRのgRainパッケージを使ってベイジアンネットワークにしてみる。 まず、ジェノタイプをベイジアンネットワーク的オブジェクトにする。 ```{r} library(gRain) Xg <- cptable(~Xg,values=c(0.01,0.18,0.81),levels=c("MM","Mm","mm")) Xg ``` <pre> ## {v,pa(v)} : Xg </pre> は、ここには、確率変数自身=vと、それに影響を与えうる親変数(pa(v))とを納めるが、この確率変数Xgの場合は、Xgという名前であって、それ以外の変数(pa(v))は一つもない、ということを宣言している。 次に <pre> ## levels of v : MM, Mm, mm </pre> は、このXgにはMM,Mm,mmの３通りの値がありえることを示している。 その次の <pre> ## values : 0.01, 0.18, 0.81 </pre> は、３通りの値に事前確率が与えられていることを表している。 最後の <pre> ## normalize=TRUE, smooth=0.000000 </pre> はオブジェクトを作るときの条件の記述だが、今は無視してよい。 次に、Xgに影響を受けるXpのためのオブジェクトを作る。 ```{r} Xp <- cptable(~Xp:Xg,values=c(1,0,0,1,0,1),levels=c("1","0")) Xp ``` 今度は、 <pre> ## {v,pa(v)} : Xp, Xg </pre> となって、この確率変数は、自身の名前はXpと言い、それには親変数(pa(v))が１つあって、それはXgであることを宣言している。 Xpの取りうる値は、levels of vが1,0である、として示されている。 次のvaluesには6つの数がある。これは、条件付き確率を表した以下の表の6つの数を、左上→左下→中上→中下→右上→右下の順に並べたものである。 ||X=MM|X=Mm|X=mm| |:--:|:--:|:--:|:--:| |P=1|1|0|0| |P=0|0|1|1| この２つしか確率変数がない単純なベイジアンネットワークの構成要素はこれですべてである。 これらを使って、ベイジアンネットワークオブジェクトにする。 ```{r} # リストにする x.list <- list(Xg,Xp) # ネットワークの個々の要素の相互関係をチェックして束ねる cptlist <- compileCPT(x.list) # ベイジアンネットワーク的に今ある条件の下での確率を計算する pn <- grain(cptlist) # プロットする plot(pn) # 中味の確率を問い合わせる querygrain(pn) ``` ジェノタイプはMM,Mm,mmの確率がそれぞれ0.01,0.18,0.81であると計算され、フェノタイプは1,0の確率が0.01,0.99であると計算された。 さて、ここで、ジェノタイプを測定したところ、MMだったとする。 これは情報・証拠であるので、この情報をネットワークに投入する。 ```{r} # すでにあるpnというネットワークのXgという確率変数がMMという値である、という情報を入れて、pn2という新たなネットワークを作る pn2 <- setEvidence(pn,"Xg","MM") ``` その中身を問い合わせると、フェノタイプが1である確率が1であり、フェノタイプが0である確率が0である、という情報が、以下のように示される。 ```{r} querygrain(pn2) ``` では、ジェノタイプを測定する代わりに、フェノタイプがわかったとする。 ```{r} pn3 <- setEvidence(pn,"Xp","0") querygrain(pn3) ``` 今度は、ジェノタイプがMM,Mm,mmのいずれなのか、どれくらいの確率なのかが示される。 フェノタイプが0のとき、ジェノタイプがMMであることはありえないから、0になっている。 Mmなのか、mmなのかは決まらないが、もともとMm vs. mmは0.18 vs. 0.81だったのでそれらの値は、MMの可能性が消滅した分だけ少し増えている。 Rで計算すれば ```{r} c(0.18,0.81)/(0.18+0.81) ``` と計算できる。 では、条件付き確率表を少し変えて、同じことをやってみる。 ||X=MM|X=Mm|X=mm| |:--:|:--:|:--:|:--:| |P=1|0.8|0.6|0.4| |P=0|0.2|0.4|0.6| どのジェノタイプもフェノタイプ1を生じるが、MM > Mm > mmの順になっている。 ベイジアンネットワークの処理自体は、Xpの条件付き確率の指定を変えるだけなので以下のようになる。 ```{r} Xg <- cptable(~Xg,values=c(0.01,0.18,0.81),levels=c("MM","Mm","mm")) Xg Xp <- cptable(~Xp:Xg,values=c(0.8,0.2,0.6,0.4,0.4,0.6),levels=c("1","0")) Xp # リストにする x.list <- list(Xg,Xp) # ネットワークの個々の要素の相互関係をチェックして束ねる cptlist <- compileCPT(x.list) # ベイジアンネットワーク的に今ある条件の下での確率を計算する pn <- grain(cptlist) # プロットする plot(pn) # 中味の確率を問い合わせる querygrain(pn) # すでにあるpnというネットワークのXgという確率変数がMMという値である、という情報を入れて、pn2という新たなネットワークを作る pn2 <- setEvidence(pn,"Xg","MM") querygrain(pn2) # フェノタイプが0であるとわかったすると pn3 <- setEvidence(pn,"Xp","0") querygrain(pn3) ``` それぞれの出力がどういう意味なのかを確認することができるはずである。 ### 親子の場合 親子３人(pa,ma,ch)を考える。 ３人のそれぞれにジェノタイプとフェノタイプがあるので、それぞれののオブジェクトを作る。 ジェノタイプ変数とフェノタイプ変数の関係は１人のときと関係ない(遺伝以外には親が子のフェノタイプに影響を与えたりしない、と仮定している)。 両親のジェノタイプの事前確率も0.01,0.18,0.81であって、他の変数の影響を受けない点は、１人のときと同じである。 子のジェノタイプは親のジェノタイプによって以下のように規定される。 |Xg.pa|Xg.ma|Xg.ch=MM,|Xg.ch=Mm|Xg.ch=mm| |:--:|:--:|:--:|:--:|:--:| |MM|MM|1|0|0| |Mm|MM|0.5|0.5|0| |mm|MM|0|1|0| |MM|Mm|0.5|0.5|0| |Mm|Mm|0.25|0.5|0.25| |mm|Mm|0|0.5|0.5| |MM|mm|0|1|0| |Mm|mm|0|0.5|0.5| |mm|mm|0|0|1| フェノタイプとジェノタイプの関係は、以下の条件付き確率を使って、ベイジアンネットワークを作成する。 ||X=MM|X=Mm|X=mm| |:--:|:--:|:--:|:--:| |P=1|1|0|0| |P=0|0|1|1| ```{r} genotype.freq <- c(0.01,0.18,0.81) genotypes <- c("MM","Mm","mm") Xg.pa <- cptable(~Xg.pa,values=genotype.freq,levels=genotypes) Xg.ma <- cptable(~Xg.ma,values=genotype.freq,levels=genotypes) Xg.ch <- cptable(~Xg.ch:Xg.pa + Xg.ma,values = c(1,0,0,0.5,0.5,0,0,1,0,0.5,0.5,0,0.25,0.5,0.25,0,0.5,0.5,0,1,0,0,0.5,0.5,0,0,1),levels=genotypes) phenotype.values <- c(1,0,0,1,0,1) phenotypes <- c("1","0") Xp.pa <- cptable(~Xp.pa:Xg.pa,values=phenotype.values,levels=phenotypes) Xp.ma <- cptable(~Xp.ma:Xg.ma,values=phenotype.values,levels=phenotypes) Xp.ch <- cptable(~Xp.ch:Xg.ch,values=phenotype.values,levels=phenotypes) ``` 要素は作り終えたので、これをリストに入れて後はネットワークにするだけである。 ```{r} # リストにする x.list <- list(Xg.pa,Xg.ma,Xg.ch,Xp.pa,Xp.ma,Xp.ch) # ネットワークの個々の要素の相互関係をチェックして束ねる cptlist <- compileCPT(x.list) # ベイジアンネットワーク的に今ある条件の下での確率を計算する pn <- grain(cptlist) # プロットする plot(pn) # 中味の確率を問い合わせる querygrain(pn) # すでにあるpnというネットワークのXg.paという確率変数がMMという値である、という情報を入れて、pn2という新たなネットワークを作る pn2 <- setEvidence(pn,"Xg.pa","MM") querygrain(pn2) ``` 父親のジェノタイプが確定したので、それに応じて父親のフェノタイプの確率も変わっている。 父親のジェノタイプの確定は子のジェノタイプの確率に影響が及び、それを通じて子のフェノタイプの確率も変わっている。 ```{r} # Xp.maのフェノタイプが0であるとわかったすると pn3 <- setEvidence(pn,"Xp.ma","0") querygrain(pn3) ``` こちらの場合は、母親のフェノタイプが確定したので、それに応じて母親のジェノタイプの確率が変化し、その結果、子のジェノタイプの確率が変化し、それに連れて子のフェノタイプの確率も変化している。 ここでさらに、子のフェノタイプも確定させてみる。 ```{r} # Xp.maのフェノタイプが0であるとわかったことになっているネットワークpn3に、追加情報として子のフェノタイプが0であるという情報を投入する pn4 <- setEvidence(pn3,"Xp.ch","0") querygrain(pn4) ``` 子のフェノタイプ情報の投入により、子のジェノタイプの取り方が影響を受け、それを受けて、母親と父親との両方のジェノタイプ確率が影響を受けている。 母親のフェノタイプはすでに確定しているので影響されないが、父親のフェノタイプの確率が影響を受けていることが見て取れる。 # 本論 では、本論に入る。 ## 家系データの作成 分離比解析において用いた家系データ作成を踏襲する。 まずはそのときに用いた関数を読み込む。 ```{r} my.pedigree <- function(n.init=3,n.iter=20,mean.kid=2){ trio <- matrix(0,n.init,5) m <- length(trio[1,]) trio[1:n.init,1] <- 1:n.init trio[1:n.init,4] <- sample(c(-1,1),n.init,replace=TRUE) trio[,5] <- 0 cnt <- n.init # 「子がいる場合」の「子の数の平均」がmean.kidなので、ポアソン分布に渡すパラメタはmean.kidより小さい tmp.f <- function(m,k){ (m/(1-exp(-m))-k)^2 } xmin <- optimize(tmp.f, c(0, mean.kid), tol = 0.0001, k = mean.kid) mean.kid. <- xmin[[1]] for(i in 1:n.iter){ s <- sample(trio[,1],1) if(length(which(trio[,2:3]==s))==0){ # 子がいない(相手がいない) trio <- rbind(trio,c(cnt+1,rep(0,m-1))) dp <- dpois(1:100,mean.kid.) dp <- dp/sum(dp) num.kids <- sample(1:100,1,prob=dp) tmp <- matrix(0,num.kids,m) tmp[,1] <- cnt+1+(1:num.kids) if(trio[s,4]==1){ tmp[,2] <- s tmp[,3] <- cnt+1 }else{ tmp[,2] <- cnt+1 tmp[,3] <- s } tmp[,4] <- sample(c(-1,1),num.kids,replace=TRUE) tmp[,5] <- trio[s,5] + 1 trio <- rbind(trio,tmp) trio[cnt+1,4] <- -trio[s,4] trio[cnt+1,5] <- trio[s,5] cnt <- cnt+1+num.kids }else if(trio[s,2] == 0 & trio[s,3] == 0){ # 親登録がない trio <- rbind(trio,c(cnt+1,rep(0,m-1))) trio <- rbind(trio,c(cnt+2,rep(0,m-1))) trio[cnt+1,4] <- 1 trio[cnt+2,4] <- -1 trio[c(cnt+1,cnt+2),5] <- trio[s,5]-1 trio[s,2] <- cnt+1 trio[s,3] <- cnt+2 cnt <- cnt+2 } } trio } make.genotype <- function(trio,prob){ ret <- rep(0,length(trio[,1])) # 家系図の個人を世代の降順に ord <- order(trio[,5]) # 親がいなければ、集団の確率で # 親がいればメンデルの法則で for(i in 1:length(ord)){ if(trio[ord[i],2]==0 & trio[ord[i],3]==0){ ret[ord[i]] <- sample(c(0,1,2),1,prob=prob) }else{ parents <- trio[ord[i],2:3] pat <- ret[trio[ord[i],2]]/2 mat <- ret[trio[ord[i],3]]/2 if(pat==0.5){ pat <- sample(c(0,1),1) } if(mat==0.5){ mat <- sample(c(0,1),1) } ret[ord[i]] <- pat+mat } } ret } make.affected.pedigree <- function(trio,p.g,r.g,h,prev){ # n.lociのジェノタイプを作る genotype <- matrix(0,length(trio[,1]),n.loci) for(i in 1:n.loci){ genotype[,i] <- make.genotype(trio,p.g[i,]) } # ジェノタイプ別リスクで個人リスクの総和を計算する R.mat <- genotype for(i in 1:n.loci){ tmp <- r.g[i,] R.mat[,i] <- tmp[genotype[,i]+1] } R.sum <- apply(R.mat,1,sum) var.R <- var(R.sum) var.E <- var.R*(1-herit)/herit E <- rnorm(length(trio[,1]),0,sqrt(var.E)) R.E <- R.sum + E affected <- rep(0,length(trio[,1])) #affected[which(R.E> quantile(R.E,1-prev))] <- 1 affected[which(R.E>= prev)] <- 1 return(list(pedigree=trio,genotype=genotype,gen.risk.mat=R.mat,gen.risk=R.sum,risk=R.E,affected=affected)) } # 性別が-1,1になっているが、kinsiph2パッケージでは1,2であることを要求するので、それを修正しつつ # kinship2パッケージの家系図描図関数を使う library(kinship2) plot.affected.pedigree <- function(ap){ sex <- ap$pedigree[,4] sex[which(sex==-1)] <- 2 my.ped <- pedigree(ap$pedigree[,1],ap$pedigree[,2],ap$pedigree[,3],sex,ap$affected) plot(my.ped,symbolosize=0.1) } library(igraph) # 家系を抜き出すための関数 affected.pedigree <- function(pedigree,genotype,gen.risk.mat,gen.risk,risk,affected){ return(list(pedigree=pedigree,genotype=genotype,gen.risk.mat=gen.risk.mat,gen.risk=gen.risk,risk=risk,affected=affected)) } # 作成家系図情報行列から、親子関係エッジのリストを作成する。グラフオブジェクトで扱うため。 # そのための関数 el.from.trio <- function(trio){ non.zero <- which(apply((trio[,2:3])^2,1,sum)!=0) el <- rbind(cbind(trio[non.zero,2],trio[non.zero,1]),cbind(trio[non.zero,3],trio[non.zero,1])) el } # 家系図を適当に発生して、その中の一人からグラフ上の距離がkのサンプルを取り出す # g igraphオブジェクト,aff フェノタイプ,n発端ID,kはnからの距離 # 抜き出しにあたり、親IDが登録されていたが、抜き出した中に親がいなくなった場合に親IDを0にするための処理と # 近親者として含まれたらその両親を含むことにするための処理とが、関数の中の主要処理となっている sample.aff.sub <- function(g,aff.pedigree,k,n=NULL){ # 指定発端IDがなければaffected==1の中から１名選ぶ if(is.null(n)){ n <- sample(which(aff.pedigree$affected==1),1) } # 指定IDから距離k以内のノードを取り出す nb <- neighborhood(g,k,n,mode="all")[[1]] # 取り出されたノードのうち、両親がともに取り出したノードに含まれなければ # その両親を0に換える ped <- matrix(aff.pedigree$pedigree[nb,],nrow=length(nb)) tmp <- ped[,2:3] tmp[!is.element(tmp,ped[,1])] <- 0 no.parents <- which(apply(tmp,1,sum)==0) ped[no.parents,2:3] <- 0 tmp2 <- ped[,2:3] # 片親IDのみがある場合は、入っていないもう片親を加える tobe.added <- unique(tmp2[which(!is.element(tmp2,ped[,1]))]) zero <- which(tobe.added==0) tobe.added <- tobe.added[-zero] #print(ped) #print(tobe.added) nb <- c(nb,tobe.added) ped <- aff.pedigree$pedigree[nb,] tmp <- ped[,2:3] tmp[!is.element(tmp,ped[,1])] <- 0 #no.parents <- which(apply(tmp,1,sum)==0) ped[,2:3] <- tmp affected.pedigree(ped,aff.pedigree$genotype[nb,],aff.pedigree$gen.risk.mat[nb,],aff.pedigree$gen.risk[nb],aff.pedigree$risk[nb],aff.pedigree$affected[nb]) } genotype.change <- function(g.mat){ ret <- g.mat ret[which(ret==0)] <- "0/0" ret[which(ret==1)] <- "0/1" ret[which(ret==2)] <- "1/1" ret } make.multi.family.table <- function(F){ family.id <- rep(1,length(F[[1]]$pedigree[,1])) pedigree <- F[[1]]$pedigree genotype <- genotype.change(F[[1]]$genotype) gen.risk.mat <- F[[1]]$gen.risk.mat gen.risk <- F[[1]]$gen.risk risk <- F[[1]]$risk affected <- F[[1]]$affected if(length(F)>1){ for(i in 2:length(F)){ family.id <- c(family.id,rep(i,length(F[[i]]$pedigree[,1]))) pedigree <- rbind(pedigree,F[[i]]$pedigree) genotype <- rbind(genotype,genotype.change(F[[i]]$genotype)) gen.risk.mat <- rbind(gen.risk.mat,F[[i]]$gen.risk.mat) gen.risk <- c(gen.risk,F[[i]]$gen.risk) risk <- c(risk,F[[i]]$risk) affected <- c(affected,F[[i]]$affected) } } tmp <- cbind(family.id,pedigree,affected,risk,gen.risk,gen.risk.mat,genotype) dimnames(tmp)[[2]] <- c("family","ind","father","mother","sex","generation","affection","risk","genetic_risk",paste("locus_risk",1:length(genotype[1,]),sep="_"),paste("genotype",1:length(genotype[1,]),sep="_")) tmp } ``` ## ジェノタイプが発病年齢に影響を与えるモデル このプログラムは学習用の簡易版なので、座位数は2個までが現実的な計算対象である。 以下は座位数2で実施する。 ジェノタイプが発病年齢に影響を与えるモデルでは、無限の時間のうちには、かならず発病する、という仮定をし、その発病年齢が生存期間に比べて若いとき、実際に発病が観察される、というように考える。 複数のイベントが蓄積して一定以上蓄積すると発病することをモデルにするとき、ガンマ分布やワイブル分布を用いることは現実的なので、以下では、ガンマ分布を利用して、仮のデータを作成するが、実観測データがあり、ジェノタイプ別の年齢別発病率が推定できる場合にはそれを用いることも悪くない。 ```{r} # 座位数 n.l <- 2 # 複合遺伝子型 g <- expand.grid(rep(list(0:2),n.l)) # 各座位のアレル１コピーあたりのリスク r.a <- c(8,2) # 相加モデルにて複合遺伝子型のリスクを定める R <- apply(t(g) * r.a,2,sum) # このリスクが発病率に影響を与えるものとする # その１つのモデルとしてガンマ関数をここでは採用する # もし、年齢別の発病率がジェノタイプごとに推定できるのであれば、それを用いてもよい # 以下では、0才から100才までの累積発病率をPsとして、リスクをパラメタとするガンマ関数で与えている age <- 0:100 Ps <- matrix(0,length(R),length(age)) for(i in 1:length(R)){ Ps[i,] <- pgamma(age,(max(R)-R[i]+7)*5,1) } # ジェノタイプ別の matplot(age,t(Ps),type="l",xlab="age",ylab="cummulative prevalence") # 1,...,100才での新規発病率をpsとしている # 上記のPsの差分をとっている ps <- t(apply(Ps,1,diff)) matplot(age[-1],t(ps),type="l",xlab="age",ylab="incidence") ``` ## 家系情報を作る 指定した数の遺伝子座のジェノタイプは、各座位のアレル頻度とHWEを仮定したジェノタイプ頻度を設定し、家系図上で親がいない個人のジェノタイプ頻度にはそれを用いる。 このようにしてできたジェノタイプの定まった個人についてフェノタイプ確認時年齢をランダムに与える。 そのうえで、各個人の複合ジェノタイプを確認し、フェノタイプ確認時年齢まえに発病しているかどうかをリスクに応じてランダムに決め、発病者については、ジェノタイプの年齢別発病率に応じて発病年齢をランダムに定める。 以下、そのソースである。 ```{r} # 0世代の独立個人数 n.init <- 1 # 家系図を膨らませる処理ステップ数 n.iter <- 5 # 平均子供数 mean.kid <-2 # n.loci個のリスク座位 n.loci <- n.l # それらのアレル頻度とHWEを仮定したジェノタイプ頻度 ## 適当にアレル頻度を発生させて p.a <- runif(n.loci) ## HWEに従うジェノタイプ頻度を作成 p.g <- cbind(p.a^2,2*p.a*(1-p.a),(1-p.a)^2) # 各座位のジェノタイプ別リスク r.g <- cbind(r.a*2,r.a,rep(0,n.loci)) # 遺伝率(全分散に占める遺伝要因分散の割合) herit <- 1 # リスク閾値 prev <- 0 # 家系を作り my.ped <- my.pedigree(n.init=n.init,n.iter=n.iter,mean.kid=mean.kid) # 作った家系を用いてジェノタイプ・フェノタイプを割りつける aff.pedigree <- make.affected.pedigree(my.ped,p.g,r.g,herit,prev) # 家系図上の個人に登録時年齢を割り振る n.s <- length(my.ped[,1]) s.age <- sample(10:90,n.s,replace=TRUE) # 次に、各個人のジェノタイプに照らして、登録時年齢までに発病している確率に応じ # 個人のフェノタイプを定める aff <- rep(0,n.s) # 複合ジェノタイプごとにリスクが決まっているので、それを確認する R.type <- rep(0,n.s) for(i in 1:length(R.type)){ R.type[i] <- which(R==aff.pedigree$gen.risk[i])[1] } # リスクごとに累積発病率に照らして、乱数を用いてフェノタイプを割り当てる for(i in 1:n.s){ if(Ps[R.type[i],s.age[i]] > runif(1)){ aff[i] <- 1 } } # 発病している個人については、発病年齢を与える # 登録年齢までのうちのいつ発病しているかは、その個人の年齢別発病率に応じてランダムに決める onset.age <- rep(-1,n.s) for(i in 1:n.s){ if(aff[i]==1){ onset.age[i] <- sample(0:s.age[i],1,prob=ps[R.type[i],1:(s.age[i]+1)]) } } # フェノタイプを登録する aff.pedigree$affected <- aff # プロットする plot.affected.pedigree(aff.pedigree) ``` ## ベイジアンネットワークにする ### ベイジアンネットワークの基礎 ベイジアンネットワークは、確率変数をノードとしたグラフである。 ある確率変数にほかの確率変数の値が影響を与えるとき、影響を与える方の変数のノードから、影響を受ける変数のノードに向けてエッジを引く。 入ってくるエッジがない確率変数（他から影響を受けない確率変数）がどのような値をとるかを事前確率で与える。 入ってくるエッジがある確率変数（他から影響を受ける確率変数）の場合には、影響を与える確率変数の値がいくらの場合に、どのような値をどのような確率でとるか、という情報（条件付き確率）を与える。 このようにしてグラフ（ネットワーク）を作成すると、一部の確率変数の値が確定したときに、他の確率変数がとる値の確率が変化する。 ### 家系図のベイジアンネットワーク 家系図には、個人と個人の間の関係がある。 個人はジェノタイプとフェノタイプとを持つ。 ジェノタイプもフェノタイプも確率変数であり、単純に考えるとき、これ以外を考える必要はない。 したがって、家系図のベイジアンネットワークでは個人の数の２倍のノードのグラフとなる。 子のジェノタイプは両親のジェノタイプに影響を受けるから、家系図のベイジアンネットワークには家系図の形をなぞったようなグラフが含まれる。 また、各個人のフェノタイプは各個人のジェノタイプに影響を受けるから、すべての個人において、ジェノタイプからフェノタイプへとエッジが引かれる。 家系図上で親のいない個人のジェノタイプの確率は事前確率として与える必要がある。 親がいる子のジェノタイプは両親のジェノタイプの条件付き確率であって、これはすべての親子関係に共通である。 各個人のジェノタイプがフェノタイプにどのような条件付き確率を与えるかは、ジェノタイプとジェノタイプ別・年齢別累積発病率情報によって定まる。 以下ではRのgRainパッケージを用いて、この条件を満足するベイジアンネットワークを作成するプロセスである。 コードを読みながら、確認してほしい。 ```{r} # gRainパッケージを使う library(gRain) # 複合ジェノタイプを"g1,g2,..."という文字列にする # 発病フェノタイプを"0,1"という文字列にする g.type <- paste("g",1:(length(R)),sep="") f.type <- c("0","1") g.type f.type ``` ```{r} # 複合ジェノタイプ数 n.g <- length(g[,1]) # ジェノタイプとフェノタイプという確率変数に関するベイジアンネットワーク用の情報を一つずつ作ってリストに格納する # 空の格納リストを作る g.table <- g0.table <- f.table <- list() for(i in 1:n.s){ g.table[[i]] <- f.table[[i]] <- list() } # 一人ずつ作る for(i in 1:n.s){ # 発病者の場合は発病時年齢を、見発病者の場合は登録時年齢を確認する if(aff[i]==0){ tmp.age <- s.age[i] }else if(aff[i]==1){ tmp.age <- onset.age[i] } # 発病者の場合には、すべての複合ジェノタイプごとに、「発病年齢において発病する確率」と「発病しない確率」とをジェノタイプ条件付きのフェノタイプ確率として登録する # 未発病の場合には、すべての複合ジェノタイプごとに、「その年齢までに発病する確率」「その年齢までに発病しない確率」を用いて、ジェノタイプという条件付きのフェノタイプ確率として登録する tmp.p <- matrix(0,n.g,2) for(j in 1:n.g){ if(aff[i]==0){ tmp.p[j,2] <- Ps[j,tmp.age] tmp.p[j,1] <- 1-tmp.p[j,2] }else if(aff[i]==1){ tmp.p[j,2] <- ps[j,tmp.age] tmp.p[j,1] <- 1-tmp.p[j,2] } } # gRainパッケージにおいてノード情報を作成する処理 # 仮に、pという変数がgという変数に影響を受けるとして作成 # pの取りうる値はf.typeとして作成した"0,1" # gは複合ジェノタイプの数だけ場合があるので、gがg1だったときにフェノタイプが0の確率、1の確率、gがg2だったときにフェノタイプが0の確率、1の確率…、の順で数値ベクトルを与えるのがvalues f.table[[i]] <- cptable(~p:g,values=c(t(tmp.p)),levels=f.type) # 作成したベイジアンネットワークオブジェクトの仮のノード名p,gにFiとかG1,G2,...とかの名前を格納しなおす f.table[[i]]$vpa[1:2] <- paste(c("F","G"),i,sep="") } # 一つ見てみる f.table[[1]] ``` v,v(pa)というのは、v=本体ノード名、v(pa)=親ノード名のことであり、levelsは本体ノードがとりうる値、valuesは親ノードの値の条件付き確率をベクトルとして格納したものになっている。 次にジェノタイプに関する親子のつながりを作るための処理を行う。 すべての親子関係において、個々の座位ごとに片親のアレルコピー数がxでもう片親のアレルコピー数がyであったら、子のアレルコピー数が0,1,2になる確率はいくつか、というのは、0,1/4,1/2,1を用いて表される。 場合分けを丁寧にすれば間違うことはない。 座位数が複数になるとき、座位同士が独立であると仮定できれば、両親の複合ジェノタイプがどのような組み合わせのときに子のジェノタイプがどれにいくつの確率でなるかも丁寧に組み合わせを考えれば迷うことはない。 ただし、座位数が増えると場合の数はどんどん大きくなる。 以下では、このルールに基づき、「子のジェノタイプの確率が、親１のジェノタイプの場合ｘ親２のジェノタイプの場合のパターン別にいくつになるか」という確率を収めたベクトルを作るための処理である。 飛ばしてしまって構わない。 ```{r} # １座位につき両親の３ジェノタイプの組み合わせである９通りについて、３ジェノタイプがどのような割合で生じるかを３ｘ３ｘ３のアレイで作成する gs <- 0:2 one.locus <- array(0,rep(3,3)) for(i in 1:3){ for(j in 1:3){ if(i==1 & j==1){ one.locus[,i,j] <- c(1,0,0) }else if(i==1 & j==2){ one.locus[,i,j] <- c(1/2,1/2,0) }else if(i==1 & j==3){ one.locus[,i,j] <- c(0,1,0) }else if(i==2 & j==1){ one.locus[,i,j] <- c(1/2,1/2,0) }else if(i==2 & j==2){ one.locus[,i,j] <- c(1/4,1/2,1/4) }else if(i==2 & j==3){ one.locus[,i,j] <- c(0,1/2,1/2) }else if(i==3 & j==1){ one.locus[,i,j] <- c(0,1,0) }else if(i==3 & j==2){ one.locus[,i,j] <- c(0,1/2,1/2) }else{ one.locus[,i,j] <- c(0,0,1) } } } # ついで、ローカス数n.lについて両親の複合ジェノタイプ3^n.l x 3^n.l通りについて、それぞれが3^n.l通りの複合ジェノタイプをどのような確率でもたらすかのベクトルを作る n.l <- 2 g.complex <- expand.grid(rep(list(gs),n.l)) n.g <- length(g.complex[,1]) oyako <- array(0,rep(n.g,3)) g.comp.pair <- expand.grid(1:n.g,1:n.g) values <- c() for(i in 1:length(g.comp.pair[,1])){ oya1.compg <- g.comp.pair[i,1] oya2.compg <- g.comp.pair[i,2] oya1.g <- unlist(g.complex[oya1.compg,]) oya2.g <- unlist(g.complex[oya2.compg,]) per.locus <- list() for(k in 1:length(oya1.g)){ per.locus[[k]] <- one.locus[,oya1.g[k]+1,oya2.g[k]+1] } tmp <- as.matrix(expand.grid(per.locus),ncol=length(oya1.g)) values <- c(values,apply(tmp,1,prod)) } oya2.tandem <- values ``` 親子関係をつなぐ確率を収めたベクトルができたので、一人ずつ、ベイジアンネットワーク用のオブジェｸﾄを作る。 ソースを確認してほしい。 ```{r} # 集団のジェノタイプ頻度 gen.pop <- c(outer(p.g[1,],p.g[2,],"*")) #gen.pop <- c(outer(c(0.25,0.5,0.25),c(0.25,0.5,0.25),"*")) #gen.pop <- c(0.3^2,2*0.3*0.7,0.7^2) gen.pop <- gen.pop/sum(gen.pop) # 個人ごとのループ for(i in 1:n.s){ father <- mother <- 0 # 家系図情報の行列から父母IDを取り出し、それが何行目にあたるかを抜き出す father <- my.ped[i,2] if(father!=0){ father <- which(my.ped[,1]==father) } mother <- my.ped[i,3] if(mother!=0){ mother <- which(my.ped[,1]==mother) } # 親ID登録がある場合と、親登録IDがない場合とで場合分けする # 親登録がないときには、ダミーのノードG0.i.f, G0.i.mのように架空を親とするように作成し、それらにはジェノタイプ事前確率を与える(両親とも登録がなければ、本人ジェノタイプに事前確率を与える) if(father!=0){ if(mother!=0){ # 両親ともに登録がある場合 g.table[[i]] <- cptable(~g:gf+gm,values=oya2.tandem,levels=g.type) g.table[[i]]$vpa[1] <- paste("G",i,sep="") g.table[[i]]$vpa[2] <- paste("G",father,sep="") g.table[[i]]$vpa[3] <- paste("G",mother,sep="") }else{# 父親登録はあり母親登録はない場合 g.table[[i]] <- cptable(~g:gf+g0,values=oya2.tandem,levels=g.type) g.table[[i]]$vpa[1] <- paste("G",i,sep="") g.table[[i]]$vpa[2] <- paste("G",father,sep="") g.table[[i]]$vpa[3] <- paste("G",0,i,"m",sep="") g0.table[[g0.cnt]] <- cptable(~g,values=gen.pop,levels=g.type) g0.table[[g0.cnt]]$vpa[1] <- paste("G",0,i,"m",sep="") g0.cnt <- g0.cnt+1 } }else{ if(mother!=0){# 父親登録がなく母親登録がある場合 g.table[[i]] <- cptable(~g:g0+gm,values=oya2.tandem,levels=g.type) g.table[[i]]$vpa[1] <- paste("G",i,sep="") g.table[[i]]$vpa[2] <- paste("G",0,i,"f",sep="") g.table[[i]]$vpa[3] <- paste("G",mother,sep="") g0.table[[g0.cnt]] <- cptable(~g,values=gen.pop,levels=g.type) g0.table[[g0.cnt]]$vpa[1] <- paste("G",0,i,"f",sep="") g0.cnt <- g0.cnt+1 }else{# 父親登録も母親登録もない場合 #g.table[[i]] <- cptable(~g:g0+g02,values=oya2.tandem,levels=g.type) g.table[[i]] <- cptable(~g,values=gen.pop,levels=g.type) g.table[[i]]$vpa[1] <- paste("G",i,sep="") #g.table[[i]]$vpa[2] <- paste("G",0,i,"f",sep="") #g.table[[i]]$vpa[3] <- paste("G",0,i,"m",sep="") #g0.table[[g0.cnt]] <- cptable(~g,values=gen.pop,levels=g.type) #g0.table[[g0.cnt]]$vpa[1] <- paste("G",0,i,"m",sep="") #g0.cnt <- g0.cnt+1 #g0.table[[g0.cnt]] <- cptable(~g,values=gen.pop,levels=g.type) #g0.table[[g0.cnt]]$vpa[1] <- paste("G",0,i,"f",sep="") #g0.cnt <- g0.cnt+1 } } } ``` 以上ですべてのノードの作成は終了したので、今度は、これをベイジアンネットワークにまとめる処理をする。 ```{r} # ジェノタイプのオブジェクトとフェノタイプのオブジェクトを単一のリストに収める gf.list <- list() for(i in 1:length(g.table)){ gf.list[[i]] <- g.table[[i]] } cnt <- length(gf.list) if(length(g0.table)>0){ for(i in 1:length(g0.table)){ gf.list[[cnt+i]] <- g0.table[[i]] } cnt <- length(gf.list) } for(i in 1:length(f.table)){ gf.list[[cnt+i]] <- f.table[[i]] } # ベイジアンネットワーク化する # 個別に作ったノード関係のつじつまを確認し cptlist <- compileCPT(gf.list) # ベイジアンネットワークに登録した条件付き確率によって、各ノードの確率を計算する pn <- grain(cptlist) ``` プロットすると、ジェノタイプノードとフェノタイプノードとの関係がわかる。 ```{r} plot(pn) ``` この段階では、親のいないノードに事前確率を投入しただけで、それ以外の情報（個人のジェノタイプとフェノタイプの情報）は一切使っていない。 この状態で、それぞれのノードのジェノタイプ確率を見てみる。すべてのノードのジェノタイプ確率は同じである。 ```{r} gen.pop querygrain(pn) ``` フェノタイプノードは、年齢の影響で確率が異なる。ここで現れるフェノタイプ確率は、発病者は「既知の発病年齢のそれ」であり見発病者は「登録時年齢のそれ」である。なぜなら、条件付き確率をネットワークに登録したときの情報がそのようになっているからである。 さて、ここで、すべての個人のフェノタイプが確定し、見発病者は登録時年齢がわかっており、発病者は発病年齢がわかっているとする。 また、ジェノタイプ情報は一部がわかっているものとする。 gRainパッケージでは以下の要領で観察データ・エビデンスを投入することができて、投入するたびに、事後確率が再計算される。 ```{r} # フェノタイプ情報のノード名ev.1とフェノタイプ情報ev.2 ev.1 <- ev.2 <- c() for(i in 1:n.s){ ev.1 <- c(ev.1,paste("F",i,sep="")) ev.2 <- c(ev.2,as.character(aff.pedigree$aff[i])) #setEvidence(pn, paste("F",i,sep=""),as.character(aff.pedigree$aff[i])) } # 情報を投入して新しいオブジェクトを作る pn.F <- setEvidence(pn, ev.1,ev.2) # この段階でG1の事後確率がどのように変わったかを確認する # 事前 querygrain(pn.F,nodes=c("G1")) # 事後 querygrain(pn,nodes=c("G1")) # 次に２番以降の個人のうち５人のジェノタイプが既知として投入する g.available <- sample(2:n.s,5) # ev.g.1はそのラベル、ev.g.2はその値 ev.g.1 <- ev.g.2 <- c() for(i in 1:length(g.available)){ ev.g.1 <- c(ev.g.1,paste("G",g.available[i],sep="")) tmp.g <- aff.pedigree$genotype[g.available[i],] tmp.g.2 <- sum((tmp.g)*3^((length(tmp.g)-1):0))+1 ev.g.2 <- c(ev.g.2,paste("g",tmp.g.2,sep="")) } # フェノタイプ情報を持つオブジェクトに投入して pn.F.G <- setEvidence(pn.F,ev.g.1,ev.g.2) # G1の値を比較する # フェノタイプ・ジェノタイプ情報の投入前 querygrain(pn.F.G,nodes=c("G1")) # フェノタイプ情報のみ投入したとき querygrain(pn.F,nodes=c("G1")) # ジェノタイプ情報も投入したとき querygrain(pn,nodes=c("G1")) ``` これを用いると、ジェノタイプが確定していない個人のジェノタイプの保有パターンを事後確率として求めることができ、また、それがわかれば、ジェノタイプごとの発病年齢分布が仮定してあるので、それに応じて、年齢ごとの発病確率や、ある年齢までの累積発病確率も計算できる。
