R語言構建蛋白質網絡並實現GN算法

目錄

• R語言構建蛋白質網絡並實現GN算法
  • 1.蛋白質網絡的構建
  • 2.生物網絡的模塊發現方法
  • 3.模塊發現方法實現和圖形展現
  • 4.附錄：igraph中經常使用函數
  • 參考連接

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) Newman Fast Greedy: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算法，咱們固然要嘗試本身實現一個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包構建網絡圖——詳細展現構建圖的基本過程

3.官方R語言igraph說明文檔

4.官方R語言手冊

5.R包igraph探究

6.模塊度(Modularity)與Fast Newman算法講解

7.模塊發現算法綜述

8.R語言的igraph畫社交關係圖示例