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)