R语言构建蛋白质网络并实现GN算法

1.蛋白质网络的构建

我们使用与人类HIV相关的蛋白质互作数据hunam-HIV PPI.csv来构建这个蛋白质互作网络。

想要读取csv文件，我们需要：

• 设置工作目录
  
• 读取CSV文件

代码如下：

setwd("/Users/.../Documents/...")
data <- read.csv("HIV-human PPI.csv")  

这样，我们就得到了蛋白质互作数据并存储在了data中。

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算法的基本思想：

可以看到，这个算法的思想非常简单。但是，这个算法什么时候终止，才能使得社群划分的结构最优？在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该方法通过直接优化模块度来发现模块。

一个例子：

cfg <- cluster_fast_greedy(g)
plot(cfg, g)  

2）GN算法：edge.betweenness.community该方法通过中间中心度找到网络中相互关联最弱的点，删除它们之间的边，并以此对网络进行逐层划分，就可以得到越来越小的模块。在适当的时候终止这个过程，就可以得到合适的模块划分结果。

调用这个方法并将其图形展示和保存的代码如下：

##
#• 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

##
#• 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 画图函数

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.模块发现算法综述