蛋白质相互作用系列：GN算法

 随着酵母双杂交、基于质谱的串联亲和纯化、蛋白质芯片等高通量生物实验方法的发展和文本挖掘技术在蛋白质相互作用预测中的广泛应用，所以蛋白质相互作用数据库日益丰富 ，并形成了越来越多的蛋白质相互作用网络(PPI)。
 所以面对大量可用的PPI网络数据，如何快速、有效地识别各种具有生物学功能的功能模块（也称蛋白质复合物）就成为蛋白质组学研究中一项极为关键的科学问题。并且在校期间的一个科研项目就正是关于蛋白质相互作用，所以决定结合目前的文献，新开一个蛋白质相互作用的系列，用来介绍相关的内容。
 传统挖掘网络结构的方法有层次聚类法(hierarchical clustering),其方法是基于各个节点之间连接的相似性，把网络自然的划分为各个子网络，根据从网络中移除边和添加边，将算法分为分裂和凝聚两类。
 在介绍GN算法前，先介绍凝聚法中一个传统的算法。其方法是先计算网络中每对顶点i,j的权重w_ij,这个权重代表顶点连接的紧密程度，然后单独拿出所有的顶点，按权重的顺序依次在两个顶点对之间添加边，从具有最大权重的边到最小权重的边。然后结果生成一个系统树状图(dendrogram)。

 
 该图是树状图，从图中可明显的看出，其层次聚类的结果。综上可知，我们要用此法时候必须是知道权重的，然而我们通常会遇到无权重的图，所以我们得想方法将该图赋予权重。其中一个方法是计算节点无关的路径数，该定义：连接同一对顶点的两条路径除了它们起始顶点和最终顶点之外不共享任何相同的顶点，则这称有两个独立路径数目。
 结合上面所述的传统方法，我们构建一个网络，来手动模拟一下该算法。
 
 上图是给定的网络图，顶点内是顶点的编号。
 
 以上表格是根据上文所提的方法来确定权重，即计算节点不相关的路径。为了方便起见，单元格内容空白是表示权重为1。所以根据构建的方法，先构建一个含有13个顶点的无边图，之后我们依次选取权重为3，权重为2，权重为1的边加入该图。为了彰显聚类的效果，我们暂将权重为3、2的边加入无边图。如下所示,明显有两个聚类中心：
 
 以上就是一种传统的层次聚类算法，而Newman等人^[1]提出了GN算法，可以很好避免传统方法的一些缺陷。在了解该算法前，我们得事先知道“介数”(betweennes)。参考Freeman（1977）^[2]的文献，他解释了顶点的介数，即穿过该顶点的连接其他顶点对的最短路径数。说白了，顶点的介数是衡量了顶点的中心性。何为中心性？如果i,j的所有最短路径都经过了k,那么认为k处于中心的位置，即k掌握着从i到j传递的信息，或者说k是一个中心枢纽。如果我们统计所有经过k的最短路径数目，那么此数目就衡量了k是有多么重要。而此数目就是k的介数。
 
 上图k的介数就是3（只考虑起点为i,终点为j）
 而GN算法是用到了边的介数。边的介数类似于顶点的介数，只不过是统计的是通过该边的最短路径数。该算法是属于分裂算法，它是不断从图中移除介数最大的边，来揭示图的模块结构。它的算法流程为：
  1.计算网络中所有边的介数;
 2.找到介数最大的边并将它从网络中删除;
 3.重新计算受删除边影响的边的介数;
  4.重复步2，直到所有边都被删除
手动模拟如下：

 由图可知，箭头所指的边是具有最大的介数：17，所以我们下一步就要删除这条边。

 截至此步，我们就可以发现该图的两大模块了。
 假设一个有m条边，n个顶点的图，则计算一次所有边的介数的时间复杂度为O(mn),而经过第三步的迭代，则最终的时间复杂度为O(m²n)。所以在面对大规模的ppi网络，该算法的时间复杂度就不尽如人意了。那么由此我们可能提出一个解决方案，即一次算出所有的边介数，然后将介数最大的边到介数最小的边依次删除，所以就只对边介数进行一次计算，其时间复杂度为O(mn)。
 但是该方案在某种情况下是行不通的。如下图，

 显然边ik的介数是最大为17，而边ig或边gk的介数为1,边kj的介数为5。如果按上述的减少时间复杂度的方法的话，先删除了边ik,而后删除边kj,最后才删除了边ig和gk。

 显然这是不对的。因为按GN算法来执行时，是删除边ik后，再次计算剩余边的介数，发现边ig和边gk为目前最大介数边，所以第二步就删除了边ig和边gk。如下图：

 故上述提出复杂度为O(mn)的方法在该种情况下的效果是不理想的。
 那么接下来我们来探讨一下如何加速该算法。首先，我们从计算边的介数的角度来考虑。如何让计算机执行该步骤？首先，我很快想到了可以使用广度优先搜索算法（BFS），它可以很好的求解单源最短路径问题。可以由该算法形成一个包含最短路径的树。如图：

 上图中的数字就是对应的边介数。（从底部开始计算，每向上一层就加一。）
 但是上图的树是最简单的一种树，当有多条最短路径经过同一顶点时，这时候的计算就有些棘手了。如图：

 上图中的顶点d就有两条不同的最短路径经过。综合Freeman（1977）^[2]、Newman(2004)^[3]Brandes(2001)^[4]等人的论文，将此情况下计算边介数的方法个人理解如下：
 假设信息要从a点流到d点，它有两个最短路径选择，即a->b->d或a->c->d（d点有两条最短路径通过）。所以流经bd边的概率为二分之一，所以该边被利用的概率为二分之一。我们继续观察f,可知有三条最短路径结束于f。有两条最短路径由d点到f,一条最短路径由e到f。所以边df的被用到的概率为三分之二，边fg的被用到的概率为三分之一。
 所以认为边df的介数为三分之二，边ef的介数为三分之一，而边bd的介数为$\frac{1}{2}\times \left(1+\frac{2}{3}\right)=\frac{5}{6}$，其中“1+”是边bd是b点到d点的一个最短路径。当计算边ab的介数时，我们发现b点仅位于一条最短路径上，所以边ab的被利用概率为1，故 $1\times \left(1+\frac{5}{6}\right)=\frac{11}{6}$
  总之，运用BFS后，我们需要先计算出各个顶点的介数，即经过该顶点的最短路径数。如图：

 上图中的数字就是对应点的介数，我们将其设置为点的权重，如w_f=3、w_d=2。计算边介数时，我们从树的顶部向上依次计算。首先计算“叶节点”所邻接的边。边df=w_d/w_f=2/3、边eg=w_e/w_g=1/1、边ef=w_e/w_f=1/3。其次计算上一层的边，如边bd=w_b/w_d*(1+边df的介数)=$\frac{1}{2}\times \left(1+\frac{2}{3}\right)=\frac{5}{6}$。最后依此公式计算上一层的边。最后结案结果如图：

 综上，GN算法可以运用BFS算法解决，当对所有顶点均使用一次BFS时，就得到本次的所有边的介数，然后删除最大介数边，然后再对所有的顶点使用BFS,不断的迭代，直至算法结束。所以根据此算法的设计思想来考虑加速问题。如果一个图中有50个顶点，在计算所有边的介数时候，我们得运行50次BFS，而BFS是可以独立、同步进行的，如果我们利用50个处理器并行计算，那么速度应该会大大加快。在Lonardi（2007）^[5]等人研究中发现：“并行实现使在包含7,000个顶点和20,000条边的图上运行群集​算法成为可能，如果在​16处理器上运行，则不到7小时；如果在32个处理器上运行，则不到5小时；如果在单个处理器上运行，则几乎需要3​天”。
 那根据上面的思路，这算法啥时候停止呢？如果我们输入一个我们未知的网络图，我们如何衡量算法输出层次结构的好坏？参考Newman(2004)^[3]的论文可知，他们介绍了两种方法。
 第一种方法：考虑一种划分将原网络划分才成k个社团（community)，构造一个k*k的对称矩阵e，e[ij]代表从i社团的节点到j社团的节点的所有边数占所有总边数的比例。这样e(ii)就代表i社团的所有内部边数占总边数的比例。若这个矩阵的迹（即e[11]+e[22]+…）最大，则代表了社团的内部联系最大，这就是符合社团特征的结果，可以作为社团的评价条件。但这种方法有缺陷，比如将整个网络作为一个社团，那自然e是最大的，但这显然不代表整个网络是一个社团。
 第二种方法：定义矩阵e的行和或列和为$a_i={\Sigma }_j{e}_{ij}$，这个和表示社区i中的边数占总边数的比例。所以有$e_{ij}=a_i{a}_j$。之后其定义了标准-模块性(Moldularity),$Q={\sum }_i\left(e_{ii}-a_i^2\right)=\mathrm{Tr}\mathrm{e}-\Vert {\mathrm{e}}^2\Vert$。Q的取值为0~1，并且Q越大，越能代表真实的分层情况。
 下面来探讨一下Q值是从何而来的。上文中的Q值定义是参考Newman(2004)^[3]，原文中并未对Q值的由来进行解释。为了弄清Q值，我又继续查询了Newman(2006)^[6]的文献，里面对Q值形成的过程进行了介绍。
 在此论文中^[6]中有:

Only if the number of between-group edges is sig­nificantly lower than would be expected purely by chance ​can we justifiably claim to have found significant community ​structure.

 即Newman认为将一个社区划分为几个社团时，只有当社团与社团间连接的边数低于纯粹出于偶然的预期边数时，才认为当前是一个好的划分。 （经实验后，这被证明是非常有效的）
 在总边数一定的情况下，上面的那句话等价于：当社团内的边数高于预期边数时，才认为这是一个好的划分。于是就引出了Q值的定义：

Q = (number of edges within communities) ​− (expected number of such edges).
社区内的边数减去此类边的预期数目

 接下来我们用式子表示社区内的边数和预期的边数。
$A_{ij}$=1或0，当i,j间有边，则为1，否则为0。
 注意：我们认为网络是无向的，若是有向网络，则将其对称处理变成无向的。
 总边数$m=\frac{1}{2}\Sigma A_{ij}$。
定义$c_i$为i所在的社团，则函数$\delta$($c_i$,$c_j$)定义为：若i,j处于同一社团，即$c_i$=$c_j$。所以社团内的边数为${\sum }_{ij}{A}_{ij}\delta \left(c_i,c_j\right)$
 关于期望边数，文献中提到了建立一个零模型(在之前关于打分矩阵的博客中，也涉及了零模型)，它本质是一种概率模型，该模型网络图与实际网络图有着相同的顶点和一样的总边数m以及每个顶点度数与原图顶点度数一样。令$p_{ij}$为顶点i，j有边的概率，即边的数目的期望值。可得：
${\sum }_{i,j}{P}_{ij}={\sum }_{i,j}{A}_{ij}=2m$

${\sum }_j{P}_{ij}=k_i$ （顶点i的度数）
令f($k_i$)为选到顶点i的概率。
$p_{ij}$=f($k_i$)f($k_j$)
$f\left(k_i\right){\sum }_{j}f\left(k_j\right)=k_i$
因为由于 ${\sum }_{j}f\left(k_j\right)$ 井不包含i, 可以看出 $f\left(k_i\right)$=$C{k}_i$,C是常量。总上有：
${\sum }_{i,j}f\left(k_i\right)f\left(k_j\right)=2m$

${\sum }_{i,j}C{k}_{i}C{k}_j=2m$
又因为：$2m=k_1+k_2+\dots +k_n$

$(2m{)}^2={\sum }_{ij}{k}_i{k}_j={\left(k_1+k_2+\dots +k_n\right)}^2$

所以$C^2{\sum }_{ij}{k}_i{k}_j=2m$
$C^2(2m{)}^2=2m$
$C=\frac{1}{\sqrt{2m}}$
$P_{ij}=f\left(k_i\right)f\left(k_j\right)=C^2{k}_i{k}_j=\frac{k_i{k}_j}{2m}$
 为了与04年定义的Q值式子一致，我们将社团内的边数转换为社团内边数占总边数的比例，即 $\frac{1}{2m}$${\sum }_{ij}{A}_{ij}\delta \left(c_i,c_j\right)$
则$Q=\frac{1}{2m}{\sum }_{ij}\left(A_{ij}-\frac{k_i{k}_j}{2m}\right)\delta (i,j)$
接着对Q进行变形
首先令v表示社团，则${\sum }_v$表示社团的累加

$\begin{array}{rl}Q& =\frac{1}{2m}\sum _{ij}\left[A_{ij}-\frac{k_i{k}_j}{2m}\right]\sum _v\delta \left(c_i,v\right)\delta \left(c_j,v\right)\\ & =\sum _v[\frac{1}{2m}\sum _{ij}{A}_{ij}\delta \left(c_i,v\right)\delta \left(c_j,v\right)\\ -& \frac{1}{2m}\sum _i{k}_i\delta \left(c_i,v\right)\frac{1}{2m}\sum _j{k}_j\delta \left(c_j,v\right)]\\ & =\sum _i\left(e_{ii}-a_i^2\right)\end{array}$
 截至此，GN算法算是介绍的差不多了。本来以为五千字能够完结，想不到却写了将近一万字。。。。中间过程不断遇到疑问，不断的查找文献，看了文献后，又多了几个疑问，真是学无止境呀。

参考文献：
[1] Girvan M , Newman M E J . Community structure in social and biological networks[J]. PNAS, 2002, 99.

[2] Freeman L C . A Set of Measures of Centrality Based on Betweenness[J]. Sociometry, 1977, 40(1):35-41.

[3] Newman M E . Fast algorithm for detecting community structure in networks[J]. Phys Rev E Stat Nonlin Soft Matter Phys, 2004, 69(6 Pt 2):066133.

[4] Brandes U . A faster algorithm for betweenness centrality*[J]. Journal of Mathematical Sociology, 2001, 25(2):163-177.

[5]Yang Q , Lonardi S . A parallel edge-betweenness clustering tool for Protein-Protein Interaction networks.[J]. International Journal of Data Mining & Bioinformatics, 2007, 1(3):241.

[6]Newman M E J . Finding community structure in networks using the eigenvectors of matrices[J]. PhRvE, 2006, 74.