GRAph ALigner Algorithm (GRAAL) を真面目にやる - 驚異のアニヲタ社会復帰の予備

MikuHatsune2013-12-17

GRAALをもう一回真面目にやってみる(Cancer Inform, J R Soc Interface, Bioinformatics)。
　
アルゴリズムとしては、ふたつのネットワーク $G(V,E)$ と $H(U,F)$ について( $G$ のほうがノード数 $V$ が少ない)、 $G$ 内のすべてのノード $v\in V$ を、 $H$ 内のノード $u\in U$ に完全に一対一対応させることを目的とする(global alignment)。
1: graphlet内にある、トポロジー的にユニークなノード(orbit)がどれくらいあるかを数える。
2: orbitのユニーク性により重み付けしながら、すべてのorbitについて相同性スコアを足し合わせる。graphletは5つまでならなんとか手で描くことができる程度の数(30通り)に収まり、orbitも73個まで数えられるから、5つまで考えることにしているっぽい。
3: ノードの密度を次数から算出して、先の相同性スコアと適当な係数比 $\alpha$ で足すことでコスト行列 $C$ を計算する。これはすべてのノードの組み合わせについて計算するので行数 $V$ 、列数 $U$ である。
4: コスト行列 $C$ のなかで最小のコストになるノードの組み合わせ $C_{min}(v\in V,\hspace{3}u\in U)$ をseedとして最初にアラインする。最小となる組み合わせが複数存在するときはランダムで選ばれるため、試行によって結果は変わりうる。
5: seedの周辺に距離 $r=1$ の球(shpere)を作り、その距離内にあるノードをアラインする。sphereの半径 $r$ は3まで大きくしてアラインしていく。
6: これでもまだ余っているノードが $G$ に存在する場合は、 $G^p(V,\hspace{3}E^p:p=\{1,\hspace{3}2,\hspace{3}3\})$ という新しくエッジを張り直したネットワークを作成する。つまり、 $G$ 内のノード $v\in V$ について、距離 $p$ にあるノードについてひたすらエッジを張り直す。感覚的には迂回路ができるようなもので、生物学的には塩基欠失や挿入に相当する。 $p=1$ のときは距離1のエッジなので、 $G^1=G$ である。この状態で4, 5を繰り返す。
　
すると $G$ のノードは $U$ 内のノードに完全に一対一に対応する。対応したノードに対して対応したエッジが推定されるので、アラインメントされたエッジの割合をEdge Correctness (EC)として統計量にする。また、この統計量を用いて、アラインメントしたネットワークがランダムに得られる確率を帰無仮説として、GO解析的に超幾何分布に従うp値を計算できる。
　
ということで、linkcommパッケージに入っているinteractomeを例にGRAALをやってみる。

library(linkcomm)
library(igraph)

p <- getLinkCommunities(pp_rnapol)
h <- getLinkCommunities(human_pp)

n1 <- p$igraph
n2 <- h$igraph

n1 <- graph.data.frame(get.data.frame(n1), directed=FALSE)
n2 <- graph.data.frame(get.data.frame(n2), directed=FALSE)

GRAALはエッジリストが必要なのでそのようにデータ変換する。

e1 <- get.data.frame(n1)
e2 <- get.data.frame(n2)

write.table(e1, "edge1.txt", row.names=FALSE, col.names=FALSE, quote=FALSE, sep=" ")
write.table(e2, "edge2.txt", row.names=FALSE, col.names=FALSE, quote=FALSE, sep=" ")

これらを gw format に変換する。

# ターミナルで
./list2leda edge1.txt >> graph1.gw
./list2leda edge2.txt >> graph2.gw

# GRAAL を行う。
# alpha は 0.8 推奨。ネットワーク1, ネットワーク2, 出力ファイル名
./GRAALRunner.py 0.8 graph1.gw graph2.gw result

# R に戻って
ealn <- read.table("result.ealn", header=FALSE)
aln  <- read.table("result.aln", header=FALSE)

newlist <- rbind(as.matrix(get.data.frame(n2)), as.matrix(ealn[,3:4])) # アラインメントしたエッジを追加した新しいエッジリスト
rownames(newlist) <- NULL

n3 <- graph.data.frame(newlist, directed=FALSE) # 新しいネットワーク
notnaidx <- !is.na(match(get.data.frame(n3, "vertices")[,1], aln$V2))

nodecol <- "yellow" # 対応するノードの色
edgecol <- "red"    # 対応するエッジの色
V(n3)$label <- NA
#V(n3)$frame.color <- replace(rep("black", nrow(get.data.frame(n3, "vertices"))), notnaidx, "red")
V(n3)$color <- replace(rep("black", nrow(get.data.frame(n3, "vertices"))), notnaidx, nodecol)
V(n3)$size <- replace(rep(0, nrow(get.data.frame(n3, "vertices"))), notnaidx, 2)
E(n3)$color <- unlist(mapply(rep, c(grey(0.6), edgecol), c(nrow(get.data.frame(n2)), nrow(ealn))))
E(n3)$width <- unlist(mapply(rep, c(1, 3), c(nrow(get.data.frame(n2)), nrow(ealn))))
par(mar=c(0,0,0,0))
# lay <- layout.auto(n3)
plot(n3, layout=lay)
legend("topleft", legend=paste("Aligned", c("node", "edge")), col=c("black", edgecol), lty=c(-1, 1), pch=c(21, -1), bty="n", merge=TRUE, lwd=3, pt.cex=2, cex=1.8, pt.bg=nodecol)