- その１「多因子の周期的定常状態」で登場した行列という話題を「グラフ理論」という話題に引き継いで、ネットワーク解析、生命科学におけるグラフ、少ないルールで複雑な構造を作るL-system、グラフの応用としてのベイジアンネットワークの臨床利用についての話題を提供します
- Rmdファイルです。html化、epub化できます(やり方はこちら)
- htmlファイル
- このシリーズの目次はこちら
```{r,echo=FALSE,warning=FALSE}
# やはり文字列で行こう
library(hash)
MakeSequence2<-function(initX,Nstep,replacer){
X<-as.list(rep(c(0),Nstep+1))
x<-initX
X[[1]]<-x
for(i in 2:(Nstep+1)){
tmp<-c()
for(j in 1:length(x)){
tmp<-unlist(strsplit(paste(paste(tmp,sep="",collapse=""), replacer[[x[j]]],sep=""),""))
}
x<-tmp
X[[i]]<-x
}
X
}
ResultTurtle2<-function(x,d=1,D=pi/2,initX=NULL,initA=NULL){
Niter<-length(x)
X<-matrix(0,Niter+1,3)
Y<-X
L<-rep(0,Niter)
if(is.null(initX)){
initX<-c(0,0)
}
if(is.null(initA)){
initA<-0
}
X[1,]<-c(initX,initA)
Stocker<-c()
for(i in 1:length(x)){
X[i+1,]<-X[i,]
if(x[i]=="F"){
X[i+1,1]<-X[i+1,1]+d*cos(X[i,3])
X[i+1,2]<-X[i+1,2]+d*sin(X[i,3])
L[i]<-1
}else if(x[i]=="f"){
X[i+1,1]<-X[i+1,1]+d*cos(X[i,3])
X[i+1,2]<-X[i+1,2]+d*sin(X[i,3])
L[i]<-0
}else if(x[i]=="+"){
X[i+1,3]<-X[i+1,3]+D
L[i]<-0
}else if(x[i]=="-"){
X[i+1,3]<-X[i+1,3]-D
L[i]<-0
}else if(x[i]=="["){
Stocker<-c(Stocker,i+1)
}else if(x[i]=="]"){
X[i+1,]<-X[Stocker[length(Stocker)],]
L[i]<-0
Stocker<-Stocker[1:(length(Stocker)-1)]
}
Y[i,]<-X[i+1,]
}
list(X=X,Y=Y,L=L,x=x)
}
DrawTurtle2<-function(ttl,col=FALSE){
xlim<-ylim<-range(ttl$X[,1:2])
plot(ttl$X[,1],ttl$X[,2],cex=0.05,asp=1,axes=FALSE,xlab="",ylab="")
for(i in 1:length(ttl$L)){
if(ttl$L[i]==1){
c<-1
if(col)c<-ttl$x[i]
segments(ttl$X[i,1],ttl$X[i,2],ttl$Y[i,1],ttl$Y[i,2],col=c)
}
}
}
elem<-c("F","[","]","+","-","X")
initX<-c("F")
replacer<-hash(elem,elem)
replacer[["F"]]<-"FF-[-F+F+F]+[+F-F-F]"
Nstep<-3
X<-MakeSequence2(initX,Nstep,replacer)
D<-22.5/180*pi
ttl<-ResultTurtle2(X[[Nstep+1]],D=D)
DrawTurtle2(ttl)
```
4要素が次のような連立微分方程式を満たすような関係にあるときに
$$
\begin{equation}
\frac{dx_1}{dt} = -x_2\\
\frac{dx_2}{dt} = -x_3\\
\frac{dx_3}{dt} = -x_4\\
\frac{dx_4}{dt} = x_1
\end{equation}
$$
こんな増減変化が現れ
```{r}
A <- matrix(c(0,-1,0,0,0,0,-1,0,0,0,0,-1,1,0,0,0),byrow=TRUE,4,4)
eigen.out <- eigen(A)
V <- eigen.out[[2]]
U <- solve(V)
s <- eigen.out[[1]]
t <- seq(from=0,to=1,length=1000)*10
X <- matrix(0,length(t),4)
x1 <- c(1,3,10,100)
for(i in 1:length(t)){
X[i,] <- V %*% diag(s^t[i]) %*% U %*% x1
}
matplot(Re(X),type="l")
```
そこには
$$
\begin{equation}
\begin{pmatrix}0 & -1 & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & -1 \\ 1 & 0 & 0 & 0 \end{pmatrix}
\end{equation}
$$
こんな行列が関わっているのでした。
この行列をよく見ると、$x_i$が$x_j$に影響を受けて減少するときには、i行j列に-1が立ち、逆に増加するときには、1が立つことを表しており、$x_i,x_j$の間に直接の関係がないときには0が立っていることを表しています。
生命現象を理解するときには、たくさんの要素の間の関係の有無を捉えて理解しようとすることがあります。
```{r}
n <- 10
M <- diag(rep(1,n))
M <- M[c(2:n,1),]
library(igraph)
g <- graph.adjacency(M)
plot(g)
```
これは代謝系に出てくる分子サイクルに見えます。
絵で描けばこのようになりますが、これを行列で表示すれば
```{r}
M
```
となります。
2つのつながり具合の異なるネットワークを描いてみます。
```{r}
g <- erdos.renyi.game(500, 1/200)
plot(g,vertex.size=5,vertex.label=NA)
g <- as.undirected(barabasi.game(500))
plot(g,vertex.size=5,vertex.label=NA)
```
どちらも、あっさりとした設定で作成したネットワークですが、見た印象はずいぶん違います。
生命現象を支えている分子ネットワークや機能ネットワークがどちらのパターンに近いのか、それともこれらとは異なるパターンなのかを明らかにすることは、構成要素の１つ１つを細かく調べることと同じくらい重要なことです。
１つ１つが分子である場合には、その構成要素を１つ１つ細かく調べるのは「分子生物学」です。
全体像を把握しに行くのは「システムバイオロジー」です。
全体像の把握にはそれに合ったアプローチがあり、ネットワーク理論やそれを支えるグラフ理論は重要です。
感染微生物はヒトの健康の大敵ですが、その微生物はたえず変化しています。その変遷の様子を捕まえるために、微生物の遺伝情報である核酸塩基配列を調べ、その継承・変遷関係を捕まえる必要があります。以下の例は、免疫低下を引き起こす疾患であるAIDSの原因ウイルス(HIVウイルス)の塩基配列変異に基づく系統分岐を描いたものです。
```{r}
library(ape)
# example tree in NH format (a string)
data("hivtree.newick")
# generate file "hivtree.phy" in working directory
cat(hivtree.newick, file = "hivtree.phy", sep = "\n")
tree.hiv <- read.tree("hivtree.phy") # load tree
unlink("hivtree.phy") # delete the file "hivtree.phy"
plot(tree.hiv,show.tip.label=FALSE)
```
ヒトの病気が家系内集積することがあり、それを引き起こしている遺伝因子や環境因子を突き止めることも重要で、家系図を解析することも有用です。
家系図もグラフの一種です。家系図では、夫婦に水平線を引き、その中点から子への線を引くので線が折れ線になりますが、父親と子、母親と子、の間に線を引くことにすれば、いわゆるグラフになります。
どちらのグラフも「絵」として眺めることも重要ですが、そこにどんな原理があるのかを解釈しようとするとき、「絵」の状態では難しいです。
しかしながら、グラフは「行列」表現できることをもう知っていますから、数学的に扱うことが可能です。
```{r,echo=FALSE}
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){
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>= prev)] <- 1
return(list(pedigree=trio,genotype=genotype,gen.risk.mat=R.mat,gen.risk=R.sum,risk=R.E,affected=affected))
}
n.init <- 10
n.iter <- 50
mean.kid <-2
n.loci <- 10
p.a <- runif(n.loci)*0.1
p.g <- cbind(p.a^2,2*p.a*(1-p.a),(1-p.a)^2)
r.a <- c(100,rep(1,n.loci-1))
r.g <- cbind(r.a*2+1,r.a+1,rep(1,n.loci))
herit <- 0.9
prev <- r.a[1]*0.9
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)
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.001)
}
plot.affected.pedigree(aff.pedigree)
```
グラフの中でもぐるりと閉じた経路がないものに「木」というタイプのグラフがあります。
多細胞生物であるヒトの体の構造には分岐木状の構造があります。
身体に張り巡らされた血管系、口から肺へと通じる気道・気管系がその例です。
単一の細胞である受精卵からだんだんに特徴を変化させながら特別な形や機能を持った細胞に分化していくことも、「多細胞が作る機能的『分岐木』」です。
生物の進化を系統樹で表したものも『分岐木』です。
植物の木が成長する様子も、体内で木構造が成長する様子も基本は同じです。
```{r,echo=FALSE,warning=FALSE}
Nstep<-3
X<-MakeSequence2(initX,Nstep,replacer)
D<-22.5/180*pi
par(mfrow=c(2,2))
for(i in 1:(Nstep+1)){
DrawTurtle2(ResultTurtle2(X[[i]],D=D))
}
par(mfrow=c(1,1))
```
こんな感じでだんだんに複雑になっていきます。
しかしながら、この複雑な構造を作るルールは単純です。
5つの記号"F","+","-","[","]"を用意し、
"F"をスタートとして、その後は
$$
\begin{equation}
"F" \to "FF-[-F+F+F]+[+F-F-F]"
\end{equation}
$$
と入れ替える、という処理を繰り返すだけで、ある特徴的な"F,+,-"の文字列を作成することができます。
できた文字列を２次元に表示する決まりを与えると、ここに示したような木構造の成長が作成できます。
入れ替え文字列のパターンを変えると、木の様子も変わります。
以下のような入れ替えパターンで描いてみます。
$$
\begin{equation}
"F" \to "F[+F]F[-F]F"\\
"F" \to "F[+F]F[-F][F]"\\
"F" \to "FF-[-F+F+F]+[+F-F-F]"\\
"F" \to "F[+X]F[-X]+X"\\
("F","X") \to ("FF","F[+X]F[-X]+X")\\
("F","X") \to ("FF","F[+X][-X]FX")\\
\end{equation}
$$
```{r,echo=FALSE,warning=FALSE}
elem<-c("F","[","]","+","-","X")
par(mfcol=c(2,3))
initX<-c("F")
replacer<-hash(elem,elem)
replacer[["F"]]<-"F[+F]F[-F]F"
Nstep<-4
X<-MakeSequence2(initX,Nstep,replacer)
D<-25.7/180*pi
ttl<-ResultTurtle2(X[[Nstep+1]],D=D)
DrawTurtle2(ttl)
initX<-c("F")
replacer<-hash(elem,elem)
replacer[["F"]]<-"F[+F]F[-F][F]"
Nstep<-4
X<-MakeSequence2(initX,Nstep,replacer)
D<-20/180*pi
ttl<-ResultTurtle2(X[[Nstep+1]],D=D)
DrawTurtle2(ttl)
initX<-c("F")
replacer<-hash(elem,elem)
replacer[["F"]]<-"FF-[-F+F+F]+[+F-F-F]"
Nstep<-3
X<-MakeSequence2(initX,Nstep,replacer)
D<-22.5/180*pi
ttl<-ResultTurtle2(X[[Nstep+1]],D=D)
DrawTurtle2(ttl)
initX<-c("X")
replacer<-hash(elem,elem)
replacer[["F"]]<-"FF"
replacer[["X"]]<-"F[+X]F[-X]+X"
Nstep<-6
X<-MakeSequence2(initX,Nstep,replacer)
D<-20/180*pi
ttl<-ResultTurtle2(X[[Nstep+1]],D=D)
DrawTurtle2(ttl)
initX<-c("X")
replacer<-hash(elem,elem)
replacer[["F"]]<-"FF"
replacer[["X"]]<-"F[+X][-X]FX"
Nstep<-6
X<-MakeSequence2(initX,Nstep,replacer)
D<-25.7/180*pi
ttl<-ResultTurtle2(X[[Nstep+1]],D=D)
DrawTurtle2(ttl)
initX<-c("X")
replacer<-hash(elem,elem)
replacer[["F"]]<-"FF"
replacer[["X"]]<-"F-[[X]+X]+F[+FX]-X"
Nstep<-4
X<-MakeSequence2(initX,Nstep,replacer)
D<-22.5/180*pi
ttl<-ResultTurtle2(X[[Nstep+1]],D=D)
DrawTurtle2(ttl)
par(mfcol=c(1,1))
```
L-systemと呼ばれる形式文法の１つで植物の成長プロセスを初めとした様々な自然物の構造を記述・表現できるアルゴリズムです。
限られた構成要素と少数の規則で複雑な構造を作り上げる仕組みの１つです。フラクタルなどもこの方法で構成することができます。
生命体・生命現象の多様性が限られた遺伝子の情報に詰め込むには、こんな仕組みも背景にありそうです。
```{r}
library(fields) # for tim.colors
munit<-10 # 解像度指定 (30は重い、20くらいが適当)
m = munit^2 # grid size
C = complex( real=rep(seq(-1.8,0.6, length.out=m), each=m ), imag=rep(seq(-1.2,1.2, length.out=m), m ) )
C = matrix(C,m,m)
L<-1
# 複素数zを0+0iから、1+0iに増やしたり
# |a+bi|一定で、偏角を0からpiに増やしたりする
#modul<-rep(0.2,L) # 偏角漸増(|a+bi|固定)
modul<-seq(from=0,to=1,length.out=L) # 実軸に長くする
#args<-seq(from=0,to=pi,length.out=L) # 偏角漸増)
args<-rep(4*pi/10,length.out=L) # 偏角固定
Y<-array(0,c(m,m,L))
for(i in 1:L){
Z =complex(m=modul[i],a=args[i]) # z(n+1)=z(n)^2+Cのz
K<-20
X = array(0, c(m,m,K))
for (k in 1:K) {
#X = array(0, c(m,m,20))
#for (k in 1:20) {
Z = Z^2+C
X[,,k] = exp(-abs(Z))
}
image(X[,,K], col=tim.colors(256)) # show final image in R
Y[,,i]<-X[,,K]
}
```
患者さんが病気の好発地域の出身者だったら、その病気かもしれない確率は高い。そんな患者さんが、その病気かもしれない症状で受診した。
検査をしてそれが異常か正常かで診断にたどり着きたい。
検査が異常になるか正常になるかは、病気の有無にも左右されるが、患者さんの生活習慣にも左右される。
また、同じような症状の別の病気でも同じような検査結果が出る。
さて、問診と検査で診断しよう。
これは普通の臨床診断プロセスです。
これにグラフを持ち込んでみる方法の１つがベイジアンネットワークです。
医療上の判断に計算機の支援を導入するのはまだ現実的とは考えられていませんが、膨大になる医学知識のすべてを活用することも１人の人間の能力を超えているのも確かです。
どのようにして計算機を活用していくかの入口を眺めておくことは決して悪いことではありません。
```{r}
library(gRain)
yn <- c("yes","no")
a <- cptable(~asia, values=c(1,99),levels=yn)
t.a <- cptable(~tub+asia, values=c(5,95,1,99),levels=yn)
s <- cptable(~smoke, values=c(5,5), levels=yn)
l.s <- cptable(~lung+smoke, values=c(1,9,1,99), levels=yn)
b.s <- cptable(~bronc+smoke, values=c(6,4,3,7), levels=yn)
e.lt <- cptable(~either+lung+tub,values=c(1,0,1,0,1,0,0,1),levels=yn)
x.e <- cptable(~xray+either, values=c(98,2,5,95), levels=yn)
d.be <- cptable(~dysp+bronc+either, values=c(9,1,7,3,8,2,1,9), levels=yn)
plist <- compileCPT(list(a, t.a, s, l.s, b.s, e.lt, x.e, d.be))
pn <- grain(plist)
pn
summary(pn)
plot(pn)
pnc <- compile(pn, propagate=TRUE)
```
「アジア出身かどうかで結核(tub)のリスクが影響される」
「喫煙者かどうか(smoke)で肺がん(lung)かどうかのリスクが影響される」
「喫煙者かどうかで慢性気管支炎(bronc)かどうかのリスクが影響される」
「結核か肺がんのどちらか(either)が気になるが、そのどちらかは結核であるかどうか、肺がんであるかどうかが決める」
「結核か肺がんかのどちらか(either)は呼吸困難感(dysp)に影響を与え、慢性気管支炎も影響を与える」
「結核か肺がんかのどちらか(either)は胸部レントゲン写真(xray)に影響を与える」
という個々の情報をグラフにしたのがこの図です。
このネットワークでは、「アジア出身かどうか」「喫煙者かどうか」「レントゲン写真の所見」などの情報を投入すると、どの診断がそれらしいかを数字にして返してくれます。
このベイジアンネットワークは、このような診断の補助的な情報提供にも使えますし、実験データから意味を読み取るための手法としても使えます。
