近邻传播聚类算法(affinity propagation clustering AP算法)^[1]是2007年由Frey等提出并发表在《Science》上的聚类算法。与其他的聚类算法相比，AP算法有较高的精度，在应用方面AP算法有较少的约束条件，但在适应不同形式的数据方面相对不足，肖宇等^[2]提出SAP(semi-supervised Affinity propagation)算法，利用数据中的先验信息，调整相似性矩阵，在先验性信息较多的情况下，算法适应性较好。AP算法特点在于，参考度的值与聚类结果相关，因此本文融合ECC^[3]与CD^[4]方法对AP算法进行改进，动态设置参考度，不再简单的设置为所有数据中位数的值。

我们进入了大数据时代，在大数据中挖掘出有用的信息，已经成为科技发展的主流趋势^[5]。单一的串行模式耗时长、效率低。为了解决这个问题，一些大数据计算平台随之产生。其中基于内存计算的Spark^[6]分布式平台倍受关注。它的生态圈十分丰富，包括特定场景下的计算库(Streaming、SQL、MLlib、Graphx)。此外其数据保存在内存中，减少了多次计算中磁盘的I/O开销。AP算法聚类的过程中，需要在矩阵之间不断地进行迭代替换，算法运算时间较长，面对海量数据，单机串行实现AP算法，不能满足实验计算需求。为了改进对处理海量数据的适应能力。本文提出基于Spark并行运算改进的AP算法SIAP(spark based on improved AP)，提高算法的运行速率。

蛋白质复合物是由2个以上功能相关且相互作用连接在一起的蛋白质组成^[7]。目前，蛋白质复合物的识别手段分为生物实验和计算方法两大类，生物实验方法持续周期长、实验成本高，而计算方法可以更好的弥补生物实验技术的不足。通过不断的发展，出现很多蛋白质复合物识别算法，MCODE算法检测网络中蛋白质链接相对密集的簇为蛋白质复合物^[8]，PCP采用FS加权的方法构造有权图，然后在加权网络中通过寻找稠密子图预测蛋白质复合物^[9]。本文应用AP算法识别蛋白质复合物。首先，PPI数据形式多样^[10]，本文对数据进行评分处理，并应用Spark平台对AP算法进行并行化处理，并给出加速比图示。最后分别在黑腹果蝇(dm)和人类(hs)2个真实PPI数据集上对本文改进的AP算法与原始的AP算法、clusterone^[11]、OH-PIN^[12]、IPCA^[13]等算法进行比较，最后引入F-Measure^[14]评价指标对结果进行评估，并给出实验效果图。

1 基于ECC与CD改进的AP算法 1.1 AP算法

AP算法的基本思想基于数据点，在数据点之间进行消息的传递，通过迭代运算自动识别聚类中心。AP的计算过程是通过不断更新吸引度矩阵和归属度矩阵实现的。

吸引度r(i, k)中, i为数据对象，k表示作为候选代表点的数据对象，r(i, k)则代表从点i发送到候选聚类中心k的数值消息，反映的是k点是否适合作为i点的聚类中心^[15]。吸引度的定义为：

$ \mathit{\boldsymbol{r}}(i,k) = \mathit{\boldsymbol{s}}(i,k) - \mathop {\max }\limits_{{k^\prime };{k^\prime } \ne k} \left\{ {\mathit{\boldsymbol{s}}\left( {i,{k^\prime }} \right) + \mathit{\boldsymbol{a}}\left( {i,{k^\prime }} \right)} \right\} $

(1)

式中：s(i, k)表示i点与k点之间的相似度；a(i, k)代表归属度，表示从候选中心k发送到i的数值消息，反映i点是否选择k作为聚类中心^[15]。归属度的定义为：

$ \mathit{\boldsymbol{a}}(i,k) = \left\{ {\begin{array}{*{20}{l}} {\sum\limits_{{i^\prime }:{i^\prime } \ne i} {\max } \left\{ {0,\mathit{\boldsymbol{r}}\left( {{i^\prime },k} \right)} \right\},k = i}\\ {\min \left\{ {0,\mathit{\boldsymbol{r}}(k,k) + \sum\limits_{{i^\prime }:{i^\prime } \notin \{ i,k\} } {\max } \left\{ {0,\mathit{\boldsymbol{r}}\left( {{i^\prime },k} \right)} \right\}} \right\},}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} k \ne i} \end{array}} \right. $

(2)

在算法的运算过程中，产生的聚类结果可能不稳定，加入震荡系数，用于调节算法的稳定性：

$ {{\mathit{\boldsymbol{r}}_{t + 1}}(i,k) = \lambda {\mathit{\boldsymbol{r}}_t}(i,k) + (1 - \lambda ){\mathit{\boldsymbol{r}}_{t + 1}}(i,k)} $

(3)

$ {{\mathit{\boldsymbol{a}}_{t + 1}}(i,k) = \lambda {\mathit{\boldsymbol{a}}_t}(i,k) + (1 - \lambda ){\mathit{\boldsymbol{a}}_{t + 1}}(i,k)} $

(4)

式中：t和t+1分别代表上一次和本次的迭代结果；λ为震荡系数。

基于式(3)、(4)，AP算法的运算过程如下：

1) 输入数据生成相似度矩阵，并设置参考度，在原始的算法中一般设为数据中位数的值。

2) 输入所需的最大迭代次数和震荡变量λ的值，构造初始为零矩阵的吸引度与归属度矩阵，并进行矩阵的迭代，更新数据。

3) 确定聚类中心，以r(k, k)+a(k, k)>0为准则，输出聚类中心以及每个点所属的类别，算法结束。

1.2 基于ECC与CD改进的AP算法

实验中应用的数据是无权的PPI数据，为了达到更好的实验结果，将PPI数据生成的网络中的边标记权重^[16]，ECC和CD是2种很有效的加权方法^[3-4]：

$ {{{{\mathop{\rm ECC}\nolimits} }_{(u,v)}} = \frac{{\left| {{N_u} \cap {N_v}} \right|}}{{\min \left\{ {\left| {{N_u}} \right|,\left| {{N_v}} \right|} \right\}}}} $

(5)

$ {{\rm{C}}{{\rm{D}}_{(u,v)}} = \frac{{2\left| {{N_u} \cap {N_v}} \right|}}{{\left| {{N_u}} \right| + \left| {{N_v}} \right|}}} $

(6)

式中N_u和N_v分别为蛋白质u和v的邻居数量。

本文将2种方法结合，使加权的精度更高。以此构造相似度矩阵。权重公式定义为：

$ {W_{(u,v)}} = \beta {\rm{EC}}{{\rm{C}}_{(u,v)}} + (1 - \beta ){\rm{C}}{{\rm{D}}_{(u,v)}} $

(7)

在构造相似度矩阵时，要设定参考度的值，传统的AP算法将对角线设为统一的定值，本文应用权重结果，动态分配对角线的值，使邻接边多的并且具有更高权重值的数据点，赋值越大，加大成为聚类中心的几率。参考度为：

$ \begin{array}{l} \mathit{\boldsymbol{s}}(i,i) = \frac{{\sum \mathit{\boldsymbol{s}} (i,k)}}{n} + {\mathop{\rm ave}\nolimits} \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} i,k,n \in m,s(i,k) > 0 \end{array} $

(8)

式中：s(i, k)为i行中大于0的值；n为当前行s(i, k)的个数; ave为所有边权重的均值。

改进的AP算法的运算步骤如下：

1) 根据式(7)、(8)，对输入的数据进行打分赋值，生成相似度矩阵，并得出参考度的值。

2) 输入最大迭代次数maxIter和震荡变量λ的值，更新迭代吸引度矩阵和归属度矩阵。

3) 以r(k, k)+a(k, k)>0为标准，输出聚类中心以及每个点所属的类别，算法结束。

2 基于Spark的改进AP算法

Spark是一个高速运算级别的，具有强大功能的开源式计算引擎。有着Hadoop的MapReduce^[17]中的优点，并扩展了MapReduce的计算模式，Spark的分布式计算更高效。Spark整体架构基于4大部分Spark Streaming、Spark SQL、Mlib、Graphx, 如图 1所示。

Spark的主要核心是弹性分布式数据集，即RDD^[18](resilient distributed datasets), RDD是一个只读的，缓存在内存中的数据集。在并行计算中应用RDD可以减少内外存的读写操作。并且支持map、reduce、collect等数据操作。此外单机模式，应用到大数据中会有很多限制，达不到运算要求。Spark以集群的模式进行部署，Spark子集群结构如图 2所示。集群中一共有2类节点，master负责资源管理，Slave则负责计算。

基于以上理论基础，本文AP算法进行并行化处理。此外在算法并行的过程中内存会受到限制，RDD应用的数据结构存储方式至关重要，本文将list类型与tuple类型相结合，以三元组的形式存储数据，生成最小内存的RDD，存储结构形式如下：

List[tuple(rowid, rowValues)], 其中rowid为行号，rowValues为第rowid行的numpy数组。

读取数据，2条数据之间有边连接，则将数据的编号和边的权值存储在List里，2条数据之间没有边连接则不存储。将相似性矩阵s, 吸引度矩阵r，归属度矩阵a数据构建RDD的操作如下所示：

s =sc.parallelize(l, NUM_PARTITIONS)；s = s.map(lambda x : f0(x, count))

其中f0为数据格式转换函数，读List进行操作，转换成数据矩阵中行的形式，并应用map函数生成一个新的RDD。

a =sc.parallelize(range(count), NUM_PARTITIO-NS)；

a = a.map(lambda t:(t, np.zeros(count)))

r =sc.parallelize(range(count), NUM_PARTITIO-NS)；

r = r.map(lambda t:(t, np.zeros(count)))

其中NUM_PARTITIONS是并行分配的参数, count为数据点的个数。

基于Spark平台SIAP算法具体执行过程如下所示：

输入：原始的PPI数据集，打分调节变量β，震荡变量λ，最大迭代次数maxIter.

1) 初始化：s [count][count]=0，r [count][count]=0，a [count][count]=0

2) 根据式(7)计算相似度，为边加权重, 并存储到s中

3) 根据式(8)更新参考度的值s[i, i]

4) 构建s矩阵，r矩阵，a矩阵RDD

5) 初始化i=1

6) 根据式(3)更新吸引度矩阵r；

根据式(4)更新归属度矩阵a；

7) l=i+1, 如果l < maxIter, 执行6)，否则执行8)

8) 如果a(i, i)+r(i, i)>0:

则第i个数据为聚类中心，存储在数组centers中

9) 数据分类，计算找出a(i, i)+r(i, i)每行所属聚类中心的最大值，将该点分配到此聚类中心。返回RDD所有元素

10) 输出结果，算法结束。

3 实验结果与分析

本文将SIAP算法应用于蛋白质复合物识别中，以此来验证SIAP聚类的加速效果，并根据F-measure将SIAP算法与AP、Cytoscape^[19]软件中的clusterone、OH-PIN、IPCA这4种算法进行对比，实验中分别采用Drosophila melanogaster(dm)^[20]数据集与Homo sapiens(hs)^[21]数据集作为测试对象。参照集本文使用Vinayagam等的蛋白质复合物数据集^[22]。

3.1 实验环境

本文选择Spark2.4.0作为数据计算框架，操作系统为Ubuntu16.04, python版本为3.5.2，Anaconda版本为4.2.0，分布式环境由5台物理机组成，master为主节点，slave1~slave4为worker节点，各节点的配置如表 1所示。

3.2 实验结果

在并行计算中采用加速比来验证效果，加速比是指，同一任务在单处理系统和并行处理系统中运行消耗时间的比率。加速比定义为：

$ S = \frac{{{T_{{\rm{serial }}}}}}{{{T_{{\rm{parallel }}}}}} $

(9)

实验应用黑腹果蝇(dm)和人类(hs)2个真实数据集分别在2、4、8、16个节点上运行SIAP算法，算法运行时间如图 3所示。

由图 3可知，针对2种数据集，随着计算节点数的增加，2种数据集的运行时间均呈现下降趋势，并且接近线性趋势。说明算法有良好的拓展性，适用于处理不同的数据集，未呈现线性的原因是在算法运行过程中，需要将数据分片处理，在不同的节点间，要进行信息交流，数据传递需要消耗一部分时间。为了使对比更加明显，针对2种数据集的算法加速比如图 4所示。

从图 4中可以看出，随着计算节点数的增加，加速比的增加趋势形式正相关，但与线性趋势有一定差距，随着节点数的增多，差距略微变大。未呈现线性的原因，同样是由于节点之间的信息传递需要一定的额外时间。通过Spark平台的应用，缩短了单线程运算的时间，充分体现的算法的并行化的速度优势。

为了验证算法的准确率，引入综合评价指标F值(F-Measure)对识别的蛋白质复合物进行评估, 准确率(precision)和召回率(recall)是F值的重要参数。准确率是指算法识别的蛋白质复合物中被匹配部分所占比例，召回率是指已知蛋白质复合物中被匹配部分所占比例，具体为：

$ {\rm{ precision }} = \frac{{{\rm{TP}}}}{{{\rm{TP}} + {\rm{FP}}}},\quad {\rm{ recall }} = \frac{{{\rm{TP}}}}{{{\rm{TP}} + {\rm{FN}}}} $

(10)

式中：TP是指算法识别出的蛋白质复合物与已知复合物匹配的数值；FP是指算法识别出的蛋白质复合物总数减去TP数的值；FN表示不能被算法识别出的蛋白质复合物匹配的已知复合物的数值。

F值的定义为：

$ F = \frac{{2 \times {\rm{ precision }} \times {\rm{ recall }}}}{{{\rm{ precision }} + {\rm{ recall }}}} $

(11)

F值是综合正确率与召回率的一种评估指标，反映整体的综合效果。F值越大，代表所得到的聚类结果越好。基于2个真实dm数据集与hs数据集，SIAP算法与原始的AP算法、clusterone、OH-PIN、IPCA这4个聚类算法进行对比，实验结果如图 5所示。

通过图 5分析可知SIAP与未改进得AP算法在两种数据集中，都有良好的聚类效果。clusterone在2种数据集中，结果相差不大，也表明了clusterone聚类算法对数据集的适应能力。针对不同的数据集SIAP算法的精度高于其他聚类算法，证明了SIAP算法的有效性。

4 结论

1) 引入ECC和CD这2种加权方法，对无权的蛋白质相互作用网络数据进行加权处理。

2) 设置动态的参考度，将参考度与边的紧密型关联起来。

3) 应用Spark平台对算法进行并行化计算，缩短了算法的运行时间，提高了算法的运行速率。

实验中将SIAP与原始的AP算法以及clusterone、OH-PIN、IPCA算法进行对比，实验结果表明，SIAP算法在不同数据集上均有最高的F值，达到了良好的聚类精度证明了算法的有效性。在接下来的研究中，将考虑生物意义，将聚类算法与生物学信息相结合，提出更有效的蛋白质复合物识别方法。