R語言構建蛋白質網路并實作GN演算法

目錄

• R語言構建蛋白質網路并實作GN演算法
  • 1.蛋白質網路的構建
  • 2.生物網路的模塊發現方法
  • 3.模塊發現方法實作和圖形展示
    • 1) 基于點連接的模塊發現：cluster_fast_greedy該方法通過直接優化模塊度來發現模塊，
    • 2) GN演算法：edge.betweenness.community
    • 3) 隨機游走：walktrap.community
    • 4）Newman快速演算法：leading.eigenvector.community
    • 5) Newman Fast Greedy:fastgreedy.community
    • 6) Fast unfolding演算法：multilevel.community
    • 7）標簽傳播演算法：label.propagation.community
    • 8）自旋玻璃社群發現：spinglass.community
    • 9）自實作GN演算法
  • 4.附錄：igraph中常用函式
    • 1）plot 畫圖函式
      • layout 設定圖的布局方式
      • vertex.size 設定節點的大小
      • vertex.color 設定節點的顏色
      • vertex.label 設定節點的標記
      • vertex.label.cex 設定節點標記的大小
      • edge.color 設定邊的顏色
      • edge.arrow.mode 設定邊的連接方式
      • edge.arrow.size 設定箭頭的大小
      • E(g)$width=1 設定邊的寬度
    • 2) 聚類分析
      • 邊的中介度聚類
      • 隨機游走
      • 特征值（個人理解覺得類似譜聚類)
      • 貪心策略
      • 多層次聚類
      • 標簽傳播
      • 檔案輸出
      • 查看變數資料型別和長度
  • 參考鏈接

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劃分模塊，

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畫社交關系圖示例

標籤：R

上一篇：pip安裝總是失敗怎么辦？ 3個方法幫你解決

下一篇：R語言dai xie