ネットワークのコミュニティ

ネットワークのコミュニティ検出をやろうやろうと思って長らく放置していたのでやる。
重複がないコミュニティはigraph, 重複があるコミュニティはlinkcommでできるので、重複があるほうをする。
linkcomm クラスというのができるので扱いにくいかと思ったら、中にigraphができるのでやりやすそうだった。しかし、linkcommのままプロットするのはパラメータの設定がめんどうだったので投げた。
と思ったら plotLinkCommGraph でできるらしい。
おなじみのラブライブの百合解析からデータをとるが、linkcommはエッジリストからのデータを受け取るようなので、エッジリストに強制変換して渡す。
コミュニティの検出は最尤推定っぽい何かで最もそれっぽいコミュニティ数にしてくれる。何かパラメータをいじれば変えられる。

# res にリンクの百合 matrix があるとして
weightmat <- res + t(res)
w <- weightmat[lower.tri(weightmat)] # weight の作成

# 隣接行列からエッジリストにする。
mat <- matrix(unlist(strsplit(unlist(mapply(function(i) mapply(function(x) paste(rownames(res)[i], x, sep=" "), rownames(res)[res[i,]>0]), seq(nrow(res)))), " ")), nc=2, byrow=TRUE)

library(linkcomm)
lc <- getLinkCommunities(mat)

   Checking for loops and duplicate edges... 100.00%
   Found and removed 9 duplicate edge(s)
   Calculating edge similarities for 21 edges... 100.00%
   Hierarchical clustering of edges...
   Calculating link densities... 100.00%
   Maximum partition density =  0.6063492 
   Finishing up...4/4... 100.00%
   Plotting...
   Colouring dendrogram... 100% 

　
オブジェクトには igraph が入っているので、こちらに慣れている人はこれでゴリ押しできる。

lc$igraph
edgecolor <- numeric(nrow(lc$edgelist))
for(i in seq(lc$cluster)){
	edgecolor <- replace(edgecolor, lc$cluster[[i]], i)
}
E(lc$igraph)$width <- w[w > 0]
E(lc$igraph)$color <- edgecolor
plot(lc$igraph)

　
コミュニティを反映してグラフを描く。味気ないのでまたまた png 重ねる。
ネットワークの次数(degree)・近接中心性(closeness)・媒介中心性(betweenness)もついでにお勉強。

stat <- round(cbind(closeness(lc$igraph), betweenness(lc$igraph)), 2)
plot(lc, type="graph")
s0 <- 0.0038 #拡大縮小率
lay0 <- matrix(unlist(locator(length(V(lc$igraph)))), nc=2) # ポチポチ座標を取得する。
for(i in seq(pics)){
	xleft=lay0[i, 1]*ra - xy0[2, i]/2*s0
	ybottom=lay0[i, 2]*ra - xy0[1, i]/2*s0
	xright=lay0[i, 1]*ra + xy0[2, i]/2*s0
	ytop=lay0[i, 2]*ra + xy0[1, i]/2*s0
	rasterImage(image=pics[[i]], xleft=xleft, ybottom=ybottom, xright=xright, ytop=ytop, xpd=TRUE)
	text(lay0[i,1]-0.1,lay0[i,2],paste(paste(c("Cl", "Be"), stat[i,]), collapse="\n"), adj=c(0,2))
}

　
こんな感じになるわけだが、もう少し真面目な解析をするなら、human_pp というHuman Protein Interactome (相互作用のビッグデータみたいなもの)がサンプルとしてあるので、これをやってみる。
が、プロットしたあとラベルが邪魔なのにどう消そうか面倒だったのでやる気があれば…

h <- getLinkCommunities(human_pp)
plot(h, type="graph")

R seminar on igraph