R語言構建蛋白質網路並實現GN演算法 1.蛋白質網路的構建 我們使用與人類HIV相關的蛋白質互作數據hunam HIV PPI.csv來構建這個蛋白質互作網路。 在R中，我們可以從存儲在R環境外部的文件讀取數據。還可以將數據寫入由操作系統存儲和訪問的文件。 R可以讀取和寫入各種文件格式，如：csv，e ...
R語言構建蛋白質網路並實現GN演算法
1.蛋白質網路的構建
我們使用與人類HIV相關的蛋白質互作數據hunam-HIV PPI.csv來構建這個蛋白質互作網路。
在R中，我們可以從存儲在R環境外部的文件讀取數據。還可以將數據寫入由操作系統存儲和訪問的文件。 R可以讀取和寫入各種文件格式，如：csv，excel，xml等。
想要讀取csv文件，我們需要：
- 設置工作目錄
- 讀取CSV文件
代碼如下：
setwd("/Users/.../Documents/...")
data <- read.csv("HIV-human PPI.csv") 這樣，我們就得到了蛋白質互作數據並存儲在了data中。
接下來，我們使用
igraph包來構建該網路。（因為數據中只有兩列表示兩個有連接的頂點，因此我沒有構建數據幀用於存放頂點的特征）
edges <- data.frame(from=data[,1],to=data[,2])
g <- graph.data.frame(edges, directed = FALSE) graph.data.frame(也可寫作graph_from_data_frame)函數有許多參數，具體內容如下:
graph_from_data_frame(edges,direced,vertices) 現在，我們已經建立了圖形g，如果你想看看它的樣子，可以簡單地通過plot(g)來做到。
2.生物網路的模塊發現方法
在許多複雜網路中，對於模塊（或稱為社區）的劃分是非常有意義的。模塊發現，或稱為社群發現主要有五種模型。
社群結構特點：社群內邊密度要高於社群間邊密度，社群內部連接相對緊密，各個社群之間連接相對稀疏。
| 社群模型 | 概念 | 效果 |
|---|---|---|
| 點連接 | 某點與某社群有關係就是某社群的 | 最差，常常是某一大類超級多 |
| 隨機游走 | 利用距離相似度，用合併層次聚類方法建立社群 | 運行時間短，但是效果不是特別好，也會出現某類巨多 |
| 自旋玻璃 | 關係網路看成是隨機網路場，利用能量函數來進行層次聚類 | 耗時長，適用較為複雜的情況 |
| 中間中心度 | 找到中間中心度最弱的刪除，並以此分裂至到劃分不同的大群落 | 耗時長，參數設置很重要 |
| 標簽傳播 | 通過相鄰點給自己打標簽，相同的標簽一個社群 | 跟特征向量可以組合應用，適用於話題類 |
其中，中間中心度的模型，為Gievan-Newman(GN)演算法的思想相同。其餘模型的詳細情況不作更多介紹，此處，參考了R語言︱SNA-社會關係網路—igraph包（社群劃分、畫圖）。
下麵，我們介紹GN演算法的基本思想：
1.計算網路中所有邊的中介中心性；
2.去除中介中心性最高的邊；
3.重新計算去除邊後的網路中所有邊的中介中心性；
4.跳至步驟2，重新計算，直至網路中沒有邊存在。
可以看到，這個演算法的思想非常簡單。但是，這個演算法什麼時候終止，才能使得社群劃分的結構最優？在Newman and Girvan 2004中，他們提出了Modularity Q（全局模塊度）的概念，進一步完善了這個演算法。一般認為，Q的取值在0.3～0.7之間最優，但是，也需具體情況具體考慮。
3.模塊發現方法實現和圖形展示
現在模塊劃分有非常多的演算法，很多都已集成在igrah中。在library("igraph")之後，我們可以調用許多包中已實現的函數對網路g劃分模塊。
| 演算法 | 作者 | 年份 | 複雜度 |
|---|---|---|---|
| GN | Newman & Girvan | 2004 | |
| CFinder | 2005 | ||
| 隨機游走方法 | Pons & Latapy | 2005 | |
| 自旋玻璃社群發現 | Reichardt & Bornholdt | 2006 | |
| LPA(標簽傳播演算法) | Raghavan et al | 2007 | O(m) |
| Fast Unfolding | Vincent D. Blondel | 2008 | |
| LFM | 2009 | O(n^2) | |
| EAGLE | 2009 | O(s*n^2) | |
| GIS | 2009 | O(n^2) | |
| HANP(Hop Attenuation & Node Preferences) | Lan X.Y. & Leung | 2009 | O(m) |
| GCE | 2010 | O(mh) | |
| COPRA | 2010 | ||
| NMF | 2010 | ||
| Link | 2010 | ||
| SLPA/GANXis(Speaker-listener Label Propagation) | Jierui Xie | 2011 | |
| BMLPA(Balanced Multi-label Propagation) | 武志昊（北交大） | 2012 | O(n*logn) |
1)基於點連接的模塊發現：cluster_fast_greedy該方法通過直接優化模塊度來發現模塊。
cluster_fast_greedy(graph, merges = TRUE, modularity = TRUE,
membership = TRUE, weights = E(graph)$weight)
graph待劃分模塊的圖。
merges是否返回合併後的模型。
modularity是否將每次合併時的模塊度以向量返回。
membership是否在每次合併時考慮所有可能的模塊結構，對應最大的模塊度計算成員向量。
weights如果非空，則是一個邊權重的向量。
return一個communities對象。
一個例子：
cfg <- cluster_fast_greedy(g)
plot(cfg, g) 2）GN演算法：edge.betweenness.community該方法通過中間中心度找到網路中相互關聯最弱的點，刪除它們之間的邊，並以此對網路進行逐層劃分，就可以得到越來越小的模塊。在適當的時候終止這個過程，就可以得到合適的模塊劃分結果。
member <-edge.betweenness.community(g.undir,weight=E(g)
$weight,directed=F)
有預設的邊權重weight，並且預設邊是無向的，directed=T時代表有向。
調用這個方法並將其圖形展示和保存的代碼如下：
##
#• Community structure in social and biological networks
# M. Girvan and M. E. J. Newman
#• New to version 0.6: FALSE
#• Directed edges: TRUE
#• Weighted edges: TRUE
#• Handles multiple components: TRUE
#• Runtime: |V||E|^2 ~稀疏:O(N^3)
##
ec <- edge.betweenness.community(g)
V(g)$size = 1 #我將大部分頂點的大小設置為1
V(g)[degree(g)>=300]$size = 5 #但度很大的頂點更大
png('/Users/.../Documents/.../protein.png',width=1800,height=1800)# 指明接下來要做的圖形的格式和長寬
plot(ec,g)
dev.off() # 關閉圖形設備
print(modularity(ec)) 這樣，圖片保存為了protein.png，還輸出了模塊度。
3)walktrap.community 利用隨機游走模型
##
#• Computing communities in large networks using random walks
# Pascal Pons, Matthieu Latapy
#• New to version 0.6: FALSE
#• Directed edges: FALSE
#• Weighted edges: TRUE
#• Handles multiple components: FALSE
#• Runtime: |E||V|^2
##
system.time(wc <- walktrap.community(g))
print(modularity(wc))
#membership(wc)
plot(wc , g) 4）Newman快速演算法：leading.eigenvector.community
Newman快速演算法將每個節點看作是一個社團，每次迭代選擇產生最大Q值的兩個社團合併，直至整個網路融合成一個社團。整個過程可表示成一個樹狀圖，從中選擇Ｑ值最大的層次劃分得到最終的社團結構。該演算法的總體時間複雜度為Ｏ（ｍ（ｍ＋ｎ））
##
#• Finding community structure in networks using the eigenvectors of matrices
# MEJ Newman
# Phys Rev E 74:036104 (2006)
#• New to version 0.6: FALSE
#• Directed edges: FALSE
#• Weighted edges: FALSE
#• Handles multiple components: TRUE
#• Runtime: c|V|^2 + |E| ~N(N^2)
##
system.time(lec <-leading.eigenvector.community(g))
print(modularity(lec))
plot(lec,g) 5)fastgreedy.community
##
#• Finding community structure in very large networks
# Aaron Clauset, M. E. J. Newman, Cristopher Moore
#• Finding Community Structure in Mega-scale Social Networks
# Ken Wakita, Toshiyuki Tsurumi
#• New to version 0.6: FALSE
#• Directed edges: FALSE
#• Weighted edges: TRUE
#• Handles multiple components: TRUE
#• Runtime: |V||E| log |V|
##
system.time(fc <- fastgreedy.community(g))
print(modularity(fc))
plot(fc, g) 6)Fast unfolding演算法：multilevel.community
##
#• Fast unfolding of communities in large networks
# Vincent D. Blondel, Jean-Loup Guillaume, Renaud Lambiotte, Etienne Lefebvre
#• New to version 0.6: TRUE
#• Directed edges: FALSE
#• Weighted edges: TRUE
#• Handles multiple components: TRUE
# Runtime: “linear” when |V| \approx |E| ~ sparse; (a quick glance at the algorithm \
# suggests this would be quadratic for fully-connected graphs)
system.time(mc <- multilevel.community(g, weights=NA))
print(modularity(mc))
plot(mc, g) 7）標簽傳播演算法：label.propagation.community
##
#• Near linear time algorithm to detect community structures in large-scale networks.
# Raghavan, U.N. and Albert, R. and Kumara, S.
# Phys Rev E 76, 036106. (2007)
#• New to version 0.6: TRUE
#• Directed edges: FALSE
#• Weighted edges: TRUE
#• Handles multiple components: FALSE
# Runtime: |V| + |E|
system.time(lc <- label.propagation.community(g))
print(modularity(lc))
plot(lc , g) 8）自旋玻璃社群發現：spinglass.community
member<-spinglass.community(g.undir,weights=E(g.undir)$weight,spins=2)
#需要設置參數weights，因為無預設值9）為了更好地理解GN演算法，我們當然要嘗試自己實現一個GN演算法。
4.附錄：igraph中常用函數
1）plot 畫圖函數
plot(g, layout = layout.fruchterman.reingold, vertex.size = V(g)$size+2,vertex.color=V(g)$color,vertex.label=V(g)$label,vertex.label.cex=1,edge.color = grey(0.5), edge.arrow.mode = "-",edge.arrow.size=5)
layout設置圖的佈局方式
layout、layout.auto、layout.bipartite、layout.circle、layout.drl、layout.fruchterman.reingold、layout.fruchterman.reingold.grid、layout.graphopt、layout.grid、layout.grid.3d、layout.kamada.kawai、layout.lgl、layout.mds、layout.merge、layout.norm、layout.random、layout.reingold.tilford、layout.sphere、layout.spring、layout.star、layout.sugiyama、layout.svd
vertex.size設置節點的大小
de<-read.csv("c:/degree-info.csv",header=F)
V(g)$deg<-de[,2]
V(g)$size=2
V(g)[deg>=1]$size=4
V(g)[deg>=2]$size=6
V(g)[deg>=3]$size=8
V(g)[deg>=4]$size=10
V(g)[deg>=5]$size=12
V(g)[deg>=6]$size=14
vertex.color設置節點的顏色
color<-read.csv("c:/color.csv",header=F)
col<-c("red","skyblue")
V(g)$color=col[color[,1]]
vertex.label設置節點的標記
V(g)$label=V(g)$name
vertex.label=V(g)$label
vertex.label.cex設置節點標記的大小
edge.color設置邊的顏色
E(g)$color="grey"
for(i in 1:length(pa3[,1])){
E(g,path=pa3[i,])$color="red"
}
edge.color=E(g)$color
edge.arrow.mode設置邊的連接方式
edge.arrow.size設置箭頭的大小
E(g)$width=1設置邊的寬度
2)聚類分析
邊的中介度聚類
system.time(ec <- edge.betweenness.community(g))
print(modularity(ec))
plot(ec, g,vertex.size=5,vertex.label=NA) 隨機游走
system.time(wc <- walktrap.community(g))
print(modularity(wc))
#membership(wc)
plot(wc , g,vertex.size=5,vertex.label=NA) 特征值（個人理解覺得類似譜聚類)
system.time(lec <-leading.eigenvector.community(g))
print(modularity(lec))
plot(lec,g,vertex.size=5,vertex.label=NA) 貪心策略
system.time(fc <- fastgreedy.community(g))
print(modularity(fc))
plot(fc, g,vertex.size=5,vertex.label=NA) 多層次聚類
system.time(mc <- multilevel.community(g, weights=NA))
print(modularity(mc))
plot(mc, g,vertex.size=5,vertex.label=NA) 標簽傳播
system.time(lc <- label.propagation.community(g))
print(modularity(lc))
plot(lc , g,vertex.size=5,vertex.label=NA) 文件輸出
zz<-file("d:/test.txt","w")
cat(x,file=zz,sep="\n")
close(zz) 查看變數數據類型和長度
mode(x)
length(x)參考鏈接
1.易百R語言教程
2.R語言igraph包構建網路圖——詳細展示構建圖的基本過程
4.官方R語言手冊
6.模塊度(Modularity)與Fast Newman演算法講解
