Directed Random Walk on Graph - ryamadaの遺伝学・遺伝統計学メモ

"Topologically inferring risk-active pathways toward precise cancer classification by directed random walk"という論文で、遺伝子ネットワークにおける部分パスウェイにスコアを与えている This paper gave scores to subgraph pathways in gene network.
その中で、観測データから、個々の遺伝子ノードにスコアを与えた上で、ネットワーク上でDirected Random Walkをさせることによって、ノイズを減らし(減ると思って)、遺伝子ノードのスコアのつけ直しを行っている In the procedure, the authors used directed random walk so that they re-scored individual genes from the initial scores based on experiments to post-random walk scores, which the authors claimed reduced the noise of the initial scores.
そのDirected Random WalkをRでやらせてみよう
- n nodes are scored with t-test at time 0 (W(0))
- W(t): Scores of node 1,2,…,n at time t
- M: Row-normalized adjacency matrix of graph
- r: (0,1),Restart probability to respect scores at time t

library(igraph)
# 適当に小さい有向グラフを作る
el <- matrix(c(1,2,1,3,2,4,3,4,4,5,3,5,2,6,3,6,1,6),byrow=TRUE,ncol=2)
g1 <- graph.edgelist(el)
plot(g1)
# その隣接行列
adj <- as.matrix(get.adjacency(g1))
adj
# どん詰まりノードや、湧き出しノードがあるとRandom walkさせるのに都合がわるいので
# その有無を確認する
No.ins <- which(apply(adj,2,sum)==0)
No.outs <- which(apply(adj,1,sum)==0)
No.ins
No.outs
# バーチャルノードを付け加えて、湧き出しノード・どん詰まりノードをなくす
# そんな隣接行列を作る
adj2 <- adj
adj2 <- cbind(adj2,rep(1,length(adj2[,1])))
adj2 <- rbind(adj2,rep(1,length(adj2[1,])))
diag(adj2) <- 0
adj2
# 隣接行列からグラフオブジェクトへ
g2 <- graph.adjacency(adj2)
plot(g2)
# 湧き出し、どん詰まりが無いことを確認
No.ins2 <- which(apply(adj2,2,sum)==0)
No.outs2 <- which(apply(adj2,1,sum)==0)
No.ins2
No.outs2

# M is row-normalized adj matrix
# 隣接行列を使って、確率的にエッジ選択をさせたい
# エッジを等確率で選ばせるにあたって、行和が１となるように補正した行列を作る
M <- adj2/apply(adj2,1,sum)
M
# Normalized?
# 行和が１であることを確認
apply(M,1,sum)
# ランダムウォークにあたり、「ランダムに歩かせ」つつ、「初期値」に戻るのとを
# 一定の割合にする(ランダムウォークだけをすると、完全平衡に達してしまうから)
r <- 0.9
# 何回繰り返す？
n.iter <- 100
# ノード数を数えておく
n.node <- length(adj2[,1])
# 各ノードのスコアを時刻ごとに記録
X <- matrix(0,n.iter,n.node)
# 初期スコアを適当に決める
library(MCMCpack)
X[1,1:(n.node-1)] <- c(rdirichlet(1,rep(1,n.node-1)))

# 繰り返す
for(i in 2:n.iter){
	X[i,] <- (1-r)*t(M) %*%as.matrix(X[i-1,],ncol=1) + r*X[1,]
}
# 収束の様子を見る
matplot(X,type="l")