随着多种生命体全基因组序列的破译和功能基因组研究的开展，生物学家将关注点逐渐转移到对蛋白质组学的研究^[1]. 蛋白质是一切生命的物质基础，是机体细胞的重要组成部分，是人体组织更新和修补的主要原料. 蛋白质并不是以个体为单位产生作用的，而是通过多种方式构成特定的蛋白质功能模块来完成特定的生命功能. 生物体内所有蛋白质的相互作用被称为蛋白质相互作用网络（protein-protein interaction network，PPIN）^[2]. 从PPIN中挖掘蛋白质功能模块不仅能够为揭示蛋白质网络的组成结构、预测蛋白质功能、解释特定的生物进程提供实验基础和理论依据^[3]，而且能够在疾病预测、药物开发等方面发挥重要的作用^[4].

为了检测生命体中的蛋白质功能模块，生物学家设计了许多生物实验方法. 近年来随着高通量生物技术的发展，海量的蛋白质相互作用数据不断涌现^[5]，生物实验方法由于检测周期长、开销大等局限性已无法满足时代发展的需求. 21世纪科技的进步带动了基于机器学习、数据挖掘等计算方法的迅速兴起，通过聚类算法高效识别PPIN中功能模块的方法逐渐发展起来^[6]. Wu等^[7]提出基于核心-附属关系的聚类方法（core-attachment based method，COACH），该方法首先检测出模块的核心蛋白质，然后将附属蛋白质逐个分派到核心蛋白质所在的模块中. Aldecoa等^[8]提出的Jerarca算法基于层次聚类思想，计算蛋白质节点间的距离权重，利用距离权重构建层次结构树，然后根据模块内和模块间节点的连接分布进行层次划分得到功能模块. Nepusz等^[9]设计了利用重叠邻域进行扩展的聚类方法（clustering with overlapping neighborhood expansion，ClusterONE），该算法能够检测潜在的重叠蛋白质功能模块，还可以识别一个功能模块完全被另一个功模块包含的情况. Liu等^[10]基于显露模式(emerging patterns，EPs)方法提出有监督聚类算法ClusterEPs，该算法通过对比正负数据类发现真实模块的特征信息，定义基于EPs的聚类评分来识别当前未知的蛋白质功能模块. 近年来，群集智能算法已成功应用于PPIN功能模块检测领域. Ji等^[11]基于蚁群寻优思想提出PPIN功能模块检测方法NACO-FMD，该方法融合了PPIN的拓扑信息和功能信息，并设计了有效的启发函数和一系列新的后处理策略来提高功能模块的检测质量. Ji等^[12]设计了适合检测大规模PPIN的多粒度模型，此模型通过缩小PPIN的规模来缩短求解时间，并在消粒度过程中利用边界节点细化机制处理模块重叠问题. Yang等^[13]提出基于细菌觅食机理的PPIN功能模块检测算法BFO-FMD，利用细菌觅食的4种生物行为设计了4种相应的有效进化机制. 由上可见，群智能算法已成为解决PPIN功能模块检测问题的有效手段.

为了提高功能模块的检测质量，探索更加高效的群智能算法，本研究模拟自然界中蝙蝠利用超声波捕食的行为，提出基于蝙蝠算法的PPIN功能模块检测方法（bat algorithm for functional module detection in PPIN，BA-FMD），通过调节声波脉冲频度协调每个蝙蝠个体的移动方向，将全局搜索和局部搜索相结合以增强算法的寻优能力和鲁棒性，期望得到更优的蛋白质网络功能模块划分.

蝙蝠是唯一能振翅飞翔的哺乳动物，大多数蝙蝠是食虫类动物. 蝙蝠是著名的“夜行侠”，虽然它的视力较差，但它却拥有超常的回声定位能力，可在黑暗中导航觅食. 蝙蝠在初始搜索猎物的过程中，发射的脉冲响度大而频度较低，随着自身与猎物距离的缩小，蝙蝠会增大脉冲发射频度并减小脉冲响度以便于精确定位猎物的空间位置^[14]. Yang^[15]将蝙蝠利用超声波捕食的行为与目标函数优化问题相联系，提出启发式随机搜索的蝙蝠算法（bat algorithm，BA）. 该算法由于具有模型简单、易于实现、潜在的并行性和分布式等特点，被广泛用于工程优化^[16-17]、聚类^[18-19]、网络模型优化^[20]、病理检测^[21]、粒子滤波器优化^[22]等领域.

蝙蝠算法有3个基本要素：脉冲频率 ${f_i}$、脉冲频度 ${r_i}$、脉冲响度 ${A_i}$. 算法将蝙蝠个体位置映射为搜索空间中的点，将优化过程模拟成蝙蝠个体为搜寻猎物而进行的移动过程，用求解问题的目标函数来衡量蝙蝠个体所处位置的优劣.

算法的具体步骤如下.

1）在目标函数优化问题中，设定义域为d维搜索空间，每个蝙蝠个体根据不同的脉冲频率来更新飞行速度和位置：

(1) $ {f_i} = {f_{\min }} + ({f_{\max }} - {f_{\min }}){\beta }, $

(2) $ {{v}}_i^t = {{v}}_i^{t - 1} + ({{x}}_i^{t - 1} - {{{x}}_*}){f_i}, $

(3) $ {{x}}_i^t = {{x}}_i^{t - 1} + {{v}}_i^t. $

式中： ${f_i}$为蝙蝠i搜索猎物时的脉冲频率； ${f_{\min}}$、 ${f_{\max}}$分别为频率的最小值和最大值；β为服从均匀分布的随机数，β∈[0,1.0]； ${{x}}_i^t$、 ${{v}}_i^t$分别为第t代第i只蝙蝠的空间位置和飞行速度， ${{x}}_i^t \!=\! {[x_{i1}^t,x_{i2}^t, \cdots,x_{id}^t]^{\rm T}}$， ${{v}}_i^t = [v_{i1}^t,$ $v_{i2}^t, \cdots ,v_{id}^t{]^{\rm T}}$， $x_{ik}^t$、 $v_{ik}^t$为第k维数值； ${{{x}}_*}$为当前种群中最优蝙蝠位置.

2）生成均匀分布随机数rand1，若 ${\rm{rand1}}>r_i^t$（第t代第i只蝙蝠的脉冲频度），则蝙蝠i在最好解周围产生新解，进行局部搜索；否则在搜索空间随机扰动，进行全局搜索. 局部搜索、全局搜索的位置更新公式分别为

(4) $ {{{x}}_{{\rm{new}}}} = {{{x}}_*} + \varphi \overline {{A^t}} , $

(5) $ {{{x}}_{{\rm{new}}}} = {{x}}_i^t + \varphi \overline {{A^t}} . $

式中： ${{{x}}_{{\rm{new}}}}$为更新后的新位置；φ∈[−1.0，1.0]，为随机数； $\overline {{{A}^t}} $为所有蝙蝠在t时刻的平均响度.

3）生成均匀分布随机数rand2，若 ${\rm{rand2}}<A_i^t$（第t代第i只蝙蝠的脉冲响度），且更新之后的位置有所改善，则接受新解并利用式（6）、（7）来减小响度、增大脉冲频度：

(6) $ A_i^{t + 1} = \alpha A_i^t, $

(7) $ r_i^{t + 1} = r_i^0[1 - \exp\; ( - \gamma t)]. $

式中：α、γ为常数，分别为响度衰减系数和脉冲频度增强系数， $0 < \alpha < 1,\gamma > 0$； $ r^0_i$为最大脉冲频度. 由式（6）、（7）可以看出，当迭代次数t趋于无穷时， $A_i^t$越来越小逐渐趋于0， $r_i^t$逐渐增大最终趋于 $r_i^0$.

4）重复上述过程，直至满足某种收敛准则从而得到全局最优解和最优值.

蝙蝠算法的基本流程如算法1所示.

算法1：蝙蝠算法

输入：目标函数 ${\rm{goal}}({{x}})$， ${{x}} ={({x_1},{x_2}, \cdots ,{x_d})^{\rm T}}$；

输出： ${\rm{goal}}({{x}})$在定义域上的最优解和最优值（此处假设为最小值）；

初始化：

设置参数值：种群数量M、最大响度 $A_i^0$、 $r_i^0$、 ${f_{\min}}$、 ${f_{\max}}$、α、γ、最大迭代次数T；

随机初始化种群，确定每个个体的 ${{{x}}_i}$、 ${{{v}}_i}$；由式（6）、（7）初始化 ${A_i}$、 ${r_i}$；

根据 ${\rm{goal}}({{x}})$得到种群中最优蝙蝠位置 ${{{x}}_*}$；

while （t < T）

1. for i=1 to M

2. 根据式（1）~（3）调整频率创建新解；

3. if ${\rm{rand1}}>r_i^t$

4. then 根据式（4）执行局部搜索；

5. else 根据式（5）执行全局搜索；

6. end if

7. 计算蝙蝠所处位置的 ${\rm{goal}}{({{{x}}_i}^t)}$；

8. if ${\rm{rand2}}<A_i^t$且 ${\rm{goal}}{({{{x}}_i}^t)} < {\rm{goal}}{({{{x}}_i}^{t - 1})}$

9. then接受新解并按照式（6）~（7）更新 $r_i^t$、 $A_i^t$；

10. end if

11. end for

12. 所有蝙蝠个体按照 ${\rm{goal}}({{x}})$排序，并找到种群中最优蝙蝠位置 ${{{x}}_*}$；

13. t++；

end while

输出最优蝙蝠个体所在的位置 ${{{x}}_*}$以及目标函数的最小值 ${\rm{goal}}({{{x}}_*})$.

BA-FMD算法利用蝙蝠算法的寻优机制完成PPIN功能模块检测. 1）在搜索空间中采用随机游走的方法初始化蝙蝠种群，每只蝙蝠所处位置代表一个候选功能模块划分，同时每个蝙蝠个体均被赋予较大的脉冲响度和较低的脉冲频度；2）在种群优化阶段，每个蝙蝠个体根据脉冲频度动态选择定向局部扰动操作或随机扰动操作，若移动后功能模块划分的效果有所改善，则接受更新之后的解并减小脉冲响度、增大脉冲频度；3）基于当前解与最优解之间的距离，以不同的脉冲频率执行自适应变异操作；4）根据自然界中优胜劣汰的进化机制执行自然选择操作，通过定向局部扰动、随机扰动、基于距离和频率的自适应变异、自然选择这4种机制的结合不断优化种群，从而得到最优或次优的功能模块划分；5）利用合并和过滤2个后处理操作对最优或次优的功能模块划分做进一步的修正.

BA-FMD中每个蝙蝠个体的位置代表一种候选的功能模块划分，将其编码为有N个有向边的图，x_i={(1→x_i1), (2→x_i2), ···, (k→x_ik), ···, (N→x_iN)}，其中，k为节点编号，x_ik为节点k的连接节点，N为PPIN中的节点总数. 采用随机游走的方式对蝙蝠个体位置进行初始化，具体步骤如下：

1）对PPIN中所有的节点从1到N进行编号，如图1所示为包括9个节点的PPIN示例图，对其从1到9进行编号；将每个节点的所有直接邻居节点按编号顺序排列建立邻居有序表，过程如图2所示.

2）随机选择节点k，根据邻居有序表和式（8）计算连接每个节点的概率，利用轮盘赌机制选择节点k的连接节点 ${x_{ik}}$，再将 ${x_{ik}}$作为起始节点，寻找 ${x_{ik}}$的连接节点. 若当前没有符合要求的节点，此节点指向自身结束游走，然后再重新选择一个未被周游过的节点重复上述操作，直到PPIN中所有节点被遍历，过程如图3（a）所示. 将 ${{{x}}_i}$中每一维元素 $(k \to {x_{ik}})$按照k节点从1到N的顺序进行排列最终形成一个蝙蝠个体位置的编码，结果如图3（b）所示.

3）对该蝙蝠个体位置进行解码得到一种功能模块划分，其中相连接的节点在同一个功能模块，不相连接的节点不在同一个模块，结果如图4所示.

在蝙蝠个体i中，假设随机选中的节点为k，则k要连接l节点的概率为

(8) $ p_{kl}^i = \left\{\begin{array}{*{20}{l}} \dfrac{{C({k,l})}}{{\displaystyle\sum\nolimits_{n \in U_k^i} {C(k,n)} }},&\;l \in U_k^i \;; \\ 0, &{\text{其他}}\;. \\ \end{array} \right. $

式中： $U_k^i$为节点集合，其中任一节点n满足以下条件：1)节点n为节点k的邻居节点；2)未被周游过；3)节点k与节点n的相似度 $C(k,n) \geqslant \varepsilon $（ε为设定的连接相似度阈值）. $C({k,l})$为节点k与节点l的相似度，是2个节点之间的结构相似性 $S({k,l})$和功能相似性 ${G}({k,l})$的平均值. 表达式如下：

(9) $ C({k,l}) = {{\left( {S({k,l}) + {G}({k,l})} \right)} / 2}, $

(10) $ S({k,l}) = {{\left| {\varGamma (k) \cap \varGamma ({l})} \right|}}/{{\sqrt {\left| {\varGamma (k)} \right|\left| {\varGamma ({l})} \right|} }}, $

(11) $ {G}({k,l}) = {{\left| {{g^k} \cap {g^l}} \right|}}/{{\left| {{g^k} \cup {g^l}} \right|}}. $

式中： $\varGamma (k)$为节点k的邻居节点集合， $\left| {\varGamma (k)} \right|$为 $\varGamma (k)$集合中元素的个数， ${g^k}$为节点k的基因本体论（gene ontology，GO）注释信息集合.

在BA-FMD算法中，为了突出蛋白质功能模块簇内连接紧密而簇间连接稀疏的特点，将模块密度函数作为蝙蝠个体i的适应度：

(12) $ {F_i} = \sum\limits_{h = 1}^H {\left[ {\frac{{{e_h}}}{{\left| E \right|}} - {{\left(\frac{{{d_h}}}{{2\left| E \right|}}\right)}^2}} \right]} . $

式中： ${F_i}$为蝙蝠个体i的适应度，适应度越大，说明当前蝙蝠所处的位置越好；H为蝙蝠个体所处位置 ${{{x}}_i}$解码后的功能模块数目； ${e_h}$为第h个功能模块中节点间的连接数； ${d_h}$为第h个功能模块中所有节点的度之和； $\left| E \right|$为整个网络中的连接数.

BA-FMD的优化过程主要包括定向局部扰动操作、随机扰动操作、基于距离和频率的自适应变异操作、自然选择操作4种机制. 其中，定向局部扰动操作和随机扰动操作的选择是由每个个体的脉冲频度高低所决定的，拥有较低脉冲频度的个体执行定向局部扰动操作，否则执行随机扰动操作.

若蝙蝠个体i的脉冲频度较低，说明它目前的进化过程并没有向着适应度增大的方向移动，且距离目标位置较远. 此时选择若干个较优蝙蝠，使蝙蝠个体i直接在较优蝙蝠附近产生新的位置，以快速提升蝙蝠个体i的适应能力. 将这种在较优解周围进行局部搜索的机制定义为定向局部扰动操作. 此策略的关键在于较优蝙蝠的选择数目. 在本研究中，较优蝙蝠数目的选择由该个体位置中每个节点与其连接节点之间的连接情况来决定. 在优化初期，蝙蝠种群分布较稀疏，脉冲频度也普遍较低，大多数个体中节点与其连接节点之间的相似度较差，此时在执行定向局部扰动操作时若选择较多的较优蝙蝠，有助于个体快速向较优解移动同时又保留了种群的多样性，并在一定程度上限制其陷入局部最优；在优化后期，随着所有个体基本都趋于较优的位置，且脉冲频度较高，大多数个体中节点与其连接节点之间的相似度较大，执行定向局部扰动操作的可能性大大降低，此时已没必要选择较多的较优蝙蝠，此时选择较少的较优蝙蝠，可以有效缩减计算量.

对于蝙蝠个体i的位置 ${{{x}}_i}$，若 ${{{x}}_i}$中节点k与其连接节点 ${x_{ik}}$之间的相似度小于相似度阈值ε，则将该节点视为连接较差节点：

(13) $ b_k^i = \left\{ \begin{array}{*{20}{l}} 1,&C({k,}{x_{{ik}}}) < {\varepsilon }\, ; \\ 0,&{\text{其他}}\, . \end{array} \right. $

式中： $ b_k^i =1$表示蝙蝠个体i的第k个节点是连接较差节点. 连接较差节点越多，该个体需要的较优蝙蝠个数越多.

定向局部扰动操作的具体步骤如下.

1）首先统计 ${{{x}}_i}$中连接较差节点的总数，按照种群数量与节点总数的比例计算该个体进行定向局部扰动操作所需的较优个体数，表达式为

(14) $ m({{{x}}_i}) = \left\lfloor {\frac{M}{N}\sum\limits_{k = 1}^N {b_k^i} } \right\rfloor \, . $

2）将m个较优个体提取出来，对于每一个节点k，统计在m个较优个体中k节点的连接节点是l的个体数，将其占m的比值作为k节点要连接l节点的概率 ${\rm{ploc}} _{kl}$，并利用轮盘赌机制决定节点k的连接节点j， ${\rm{ploc}} _{kl}$的表达式为

(15) $ {{\rm{ploc}} _{kl}} = {{\sum\nolimits_{i = 1}^m {u_{kl}^i} }/ m}, $

(16) $ u_{kl}^i = \left\{ \begin{array}{*{20}{l}} 1,&{x_{ik}} = l\, ; \\ 0,&{\text{其他}} \, . \end{array} \right. $

式中： $\displaystyle\sum\nolimits_{i = 1}^m {u_{kl}^i} $为在m个较优解中，第k个节点的连接节点是l的解的个数. 其中，对于蝙蝠个体i的位置，第k个节点的连接节点是l，则 $u_{kl}^i$为1，否则为0.

3）在j和 ${x_{ik}}$之间选择与节点k连接相似度较大的节点作为最终k的连接节点. 策略流程如算法2所示.

算法2：BA-FMD中的定向局部扰动操作

输入： ${{{x}}_i}$；

参数：N，M，ε；

输出：定向局部扰动操作后蝙蝠个体i的新位置 ${\rm{new}}\;{{{x}}_i}$；

由式（13）、（14）计算 ${{{x}}_i}$定向局部扰动所需要的较优蝙蝠个数m；

for k=1 to N

1. 由m个蝙蝠计算k节点连接每个节点的概率式(15)、(16)；

2. 根据概率利用轮盘赌机制得到k节点的最终连接节点j；

3. if（ $C(k,j) > C(k,{x_{ik}})$）then ${\rm{new}}\;{x_{ik}} = j$；

4. else ${\rm{new}}\;{x_{ik}} = {x_{ik}}$；

5. end if

end for

输出 ${\rm{new}}\;{{{x}}_i}$.

该操作的执行过程如图5所示. 图中， ${{{x}}_i}$中连接较差节点个数为4，假设M=8，N=9，则m约为3；将适应度排名前3的个体排列出来，对于节点“1”，其连接节点为“2”的概率为1/3，连接节点为“7”的概率为2/3，由轮盘赌机制决定节点“1”连接的节点为 “7”；假设C（1，2）>C（1，7），新位置 ${\rm{new}}\;{{{x}}_i}$中节点“1”的连接节点最终为节点“2”.

若蝙蝠个体i的脉冲频度较高，说明它之前的移动方向大体正确，目前也已经到达了相对较优越的位置，则只须在搜索空间随机扰动，观察是否能找到比当前状态更优的位置. 本研究将这种机制定义为随机扰动操作，其具体步骤如下：1）在 $\left( {1,{N^p}} \right)$之间生成随机数e，其中，p为（0，1.0）之间的常数，定义为随机扰动指数；2）将N维向量 ${{{v}}_1}$设为蝙蝠i的移动方向， ${{{v}}_1}$中包含e个1和N−p个0，对于 ${{{v}}_1}$中1位置对应的节点k在其邻居中随机搜索得到节点l；3）在l和 ${x_{ik}}$之间选择与节点k连接相似度较大的节点作为最终k的连接节点. 策略流程如算法3所示，执行过程如图6所示. 图中，向量 ${{{v}}_1}$的第2个位置是“1”，节点“2”在其邻居中随机搜索得到节点“7”，假设C（2，7）>C（2，1），则新位置 ${\rm{new}}\;{{{x}}_i}$中节点“2”的连接节点最终为节点“7”.

算法3：BA-FMD中的随机扰动操作

输入： ${{{x}}_i}$；

参数：N，p；

输出：随机扰动操作后蝙蝠个体i的新位置 ${\rm{new}}\;{{{x}}_i}$；

1. 生成随机数 $e = {\rm{random\;(1,}}\;{N^p}{\rm{)}}$；

2. 生成N维向量 ${{{v}}_1}$，将e个1随机分布在向量 ${{{v}}_1}$中，其余位置为0；

3. 对于 ${{{v}}_1}$中1对应的节点k，在其邻居中随机搜索得到节点l；

4. if（ $C(k,l) > C(k,{x_{ik}})$） then ${\rm{new}}\;{x_{ik}} = l$；

5. else ${\rm{new}}\;{x_{ik}} = {x_{ik}}$；

6. end if

输出 ${\rm{new}}\;{{{x}}_i}$

在种群中每个蝙蝠个体以不同的脉冲频率进行搜索，由式（2）可以看出，蝙蝠个体的飞行速度由脉冲频率和它与最优解之间的距离共同决定. 模拟此操作，将策略定义为基于距离和频率的自适应变异操作. 在每一代进化中，蝙蝠个体i与最优解的距离越大，说明该个体的进化空间越大，此时设置较大的自适应步伐有助于该个体向最优解快速迈进. 操作的具体步骤如下.

1）计算个体位置 ${{{x}}_i}$与当前最优解 ${{{x}}_*}$的距离：

(17) $ D({{{x}}_i},{{{x}}_*}) = \sum\limits_{k = 1}^N {({x_{ik}} \oplus {x_{*k}})} , $

(18) $ {x_{ik}} \oplus {x_{*k}} = \left\{ \begin{array}{*{20}{l}} 1,&{x_{ik}} = {x_{*k}}\, ; \\ 0,&{x_{ik}} \ne {x_{*k}}\, . \end{array} \right. $

式中： ${x_{ik}}$、 ${x_{*k}}$分别为个体i位置和最优个体位置的第k个节点所连接的节点.

2）在 ${f_{\min}}$与 ${f_{\max}}$之间随机产生脉冲频率 ${f_i}$，若 $D({{{x}}_i},{{{x}}_*}){f_i} > {\mu }N$，说明该个体距离最好解较远，则重新初始化该个体得到 ${\rm{new}}\;{{{x}}_i}$，其中 ${\mu } \in (0,1.0)$，为重新初始化系数；否则生成N维向量 ${{{v}}_2}$作为蝙蝠的飞行方向，其中包含 $\left\lfloor {D({{{x}}_i},{{{x}}_*}){f_i}} \right\rfloor $个1，其余位置为0， ${{{v}}_2}$中1位置对应节点k，在 ${x_{*k}}$与 ${x_{ik}}$之间选择与节点k连接相似度较大的节点作为最终k的连接节点.

具体流程如算法4所示.

算法4：基于距离和频率的自适应变异操作

输入： ${{{x}}_i}$；

参数：N， ${f_{\min}}$， ${f_{\max}}$，μ；

输出：基于距离和频率的自适应变异操作后蝙蝠个体i的新位置 ${\rm{new}}\;{{{x}}_i}$；

1. 由式（17）、（18)计算个体i与当前最优解的距离；

2. 设置个体i的频率 ${f_i}$，为（ ${f_{\min}}$， ${f_{\max}}$）之间的随机数；

3. if（ $D({{{x}}_i},{{{x}}_*}){f_i} > {\mu }N$）

4. then 重新初始化得到 ${\rm{new}}\;{{{x}}_i}$；

5. else 生成N维向量 ${{{v}}_2}$作为蝙蝠的飞行方向，将 $\left\lfloor {D({{{x}}_i},{{{x}}_*}){f_i}} \right\rfloor $个1随机散布在 ${{{v}}_2}$中，其余位置为0；

对于 ${{{v}}_2}$中1对应的节点k：

6. if（ $C(k,{x_{*k}}) > C(k,{x_{ik}})$）

7. then ${\rm{new}}\;{x_{ik}} = {x_{*k}}$；

8. else ${\rm{new}}\;{x_{ik}} = {x_{ik}}$；

9. end if

10. end if

输出 ${\rm{new}}\;{{{x}}_i}$.

将种群中的个体按照适应度进行排序，根据生物界的优胜劣汰规律，淘汰种群中η个较差个体，并将η个较优个体进行复制放入种群，这样不仅保证了种群规模，同时还提升了种群的适应能力. 将此策略定义为自然选择操作，其中η的表达式为

(19) $ {\eta } = 0.05M + {{0.2M} / t}. $

式中：t为目前迭代次数，η随着进化的深入而不断变小. 通常，种群在进化的前期具有较好的多样性，但是随机性较强，较大的淘汰数量可以迅速提升种群整体的适应性；在进化的后期，算法逐渐收敛，减小淘汰数量可以保留种群的多样性.

1）BA-FMD算法首先产生一个初始种群，种群中每个蝙蝠个体以较大的脉冲响度和较低的脉冲频度进行搜索，每只蝙蝠所处位置编码为一个候选功能模块划分. 2）在种群优化阶段，发射较低脉冲频度的蝙蝠执行定向局部扰动操作，即在选择的较好解周围进行局部搜索；发射较高脉冲频度的蝙蝠执行随机扰动操作，在搜索空间寻找更优的解；若更新之后所有功能模块的模块密度增大，则接受新解并减小脉冲响度、增大脉冲频度；3）采用基于距离和频率的自适应变异操作对候选解进行优化，根据当前解与最优解之间的距离设置向最优解迈进的自适应步伐，并以一定的可能性将该个体重新初始化为一个功能模块划分，从而有效避免种群陷入局部最优；4）对整个种群执行自然选择操作提升种群整体的适应能力，加快算法的收敛速度. 通过这4种机制的结合不断优化种群，得到最优或次优的功能模块划分. 最后利用合并和过滤2个后处理操作对最优或次优的功能模块划分做进一步的修正. BA-FMD算法的具体流程如算法5所示. 其中，G (V, E)中，G为PPIN，V为G中所有节点的集合，E为G中所有边的集合.

算法5：BA-FMD

输入：PPIN：G（V，E），N=|V|；GO信息；

输出：PPIN功能模块集合；

初始化：

设置参数值：M、最优解保持不变的最大迭代次数Q、p、 $r_i^0$、 $A_i^0$、α、γ、 ${f_{\min}}$、 ${f_{\max}}$、μ、ε、合并阈值λ、过滤阈值δ；

计算每2个节点的相似度；

初始化种群：

for i=1 to M

选择任意蛋白质节点kϵ[1，N]，根据式（8）选择l作为k的连接节点，再将l节点作为起始节点确定其连接节点，直到所有的节点被遍历，形成个体i的模块划分 ${{{x}}_i}$；

end for

蝙蝠算法优化种群：

计算每个蝙蝠个体的适应度 $F_i^0$；

最优个体适应度为 $F_{\max }^0$；

迭代次数t=1；由式（6）、（7）计算 $A_i^1$、 $r_i^1$；

最优个体适应度不变的迭代次数q=0；

while q<Q

1. for i=1 to M

2. if ${\rm{rand1}}>r_i^t$

3. then执行定向局部扰动操作；

4. else执行随机扰动操作；

5. end if

6. if ${\rm{rand2}}<A_i^t$且 $F_i^t > F_i^{t - 1}$

7. then接受新解，按照式（6）、（7）更新 $A_i^t$、 $r_i^t$；

8. end if

9. end for

10. 计算每个蝙蝠个体的适应度并排序，更新最优解；

11. for i=1 to M

12. 执行基于距离和频率的自适应变异操作；

13. end for

14. 执行自然选择操作；

15. if （ $F_{\max }^t = F_{\max }^{t - 1}$）then q++；

16. else q=0；

17. end if

18. t++；

end while

对优化得到的最好解进行合并和过滤2个后处理操作^[11]；

输出后处理之后的PPIN功能模块划分.

基于算法5的描述，对BA-FMD的时间复杂度进行分析. 在初始化过程中，计算每2个节点相似度的时间复杂度为 $O({d_{\max}} N)$，初始化种群需要 $O({d_{\max}} N M)$，其中， ${d_{\max}}$为节点中最大的度. 在种群进化过程中，时间复杂度为 $O(T (M m N + M({N^p} + N) + $ $ {M^2} + M + N))$，其中， $p \in (0,1.0)$，所以 ${N^p} < N$，可将 ${N^p} + N$看作2N，而 $m \gg 2$，所以进化过程的时间复杂度可近似为 $O(T M m N)$. 后处理过程的时间复杂度为 $O({H^2} + H) \approx $ $O({H^2})$，H为最终的检测模块数目. PPIN是具有无标度性和小世界性的复杂网络，有 ${d_{\max}} \ll N$， $H < N$，所以最终的时间复杂度主要取决于种群进化的时间复杂度，即 $O(T M m N)(m < M,$对于大规模的PPIN网络 $M \ll N$）.

本实验在处理器为Core（TM）i5-6500 CPU、RAM为4.00 GB、操作系统为Windows10的环境下，利用eclipse运行工具编写Java代码进行实现.

选用5个经典的酵母菌PPIN公共数据集：Gavin数据集、DIPcore数据集、Krogan数据集、Collins数据集和BioGRID数据集. 5个数据集节点和边的详细情况如表1所示. 算法中用到的GO信息可以通过网站bioDBnet（https://biodbnet-abcc.ncifcrf.gov/db/db2db.php）查询. 此外，为了检测本算法的有效性，利用文献[23]提供的由MIPS、Aloy、SGD这3个数据库整合而来的428个标准模块作为基准模块.

采用3组常见的评价指标来证明算法的有效性.

在检测功能模块的过程中，有可能将一些蛋白质节点丢弃或者过滤. 覆盖率（coverage）指标为检测的所有功能模块组成的蛋白质节点集合中元素个数占蛋白质节点总数的比例，表达式为

(20) $ {\rm{Cv}} = {\left.{\left| \scriptstyle{\bigcup\nolimits_{i = 1}^{\left| A \right|} {{V_i}} } \right|} \right/ N} \,. $

式中： $\left| A \right|$为检测的功能模块数目， ${V_i}$为第i个模块中的蛋白质节点集合. 覆盖率越小，说明算法丢弃的蛋白质节点越多；相反，覆盖率越大，说明保留了较多的蛋白质节点.

精度（precision）、召回率（recall）和F度量（F-measure）是检测PPIN功能模块问题中常用的评价指标. 在计算这些评价指标之前首先要计算邻域亲和评分（neighborhood affinity，NA），此评分衡量检测模块 $a = ({V_a},{E_a})$与标准模块 $b = ({V_b},{E_b})$的匹配程度（ $V$、E分别为模块中蛋白质节点集合和边连接集合），表达式为

(21) $ {\rm{NA}}({a,b}) = \frac{{{{\left| {{V_a} \cap {V_b}} \right|}^2}}}{{\left| {{V_a}} \right|\left| {{V_b}} \right|}}. $

若 ${\rm{NA}}({a,b}) \geqslant {\omega }$（一般取 ${\omega }=0.2$），则认为检测模块a和标准模块b相匹配. 令A为通过计算方法检测到的功能模块集合，B为标准模块集合，则A中至少与一个标准模块匹配的模块数量为 ${N_{a}} = \left| {\left\{ {a|a \in A,\exists b \in B,{\rm{NA}}\left( {a,b} \right) \geqslant \omega } \right\}} \right|$，B中至少与一个检测模块匹配的模块数量为 ${N_{b}} = {| {{\{ }}b\left| {b \in } \right.B, }$ ${\exists a \in A} , $ ${\rm{NA}}({a,b}) \geqslant {\omega }{{\} }} |$. 则精度和召回率的表达式分别为

(22) $ {\rm{Pr}} = {{{N_{a}}}}/{{\left| A \right|}}, $

(23) $ {\rm{Rc}} = {{{N_{b}}}}/{{\left| B \right|}}\, . $

式中：Pr为精度，Rc为召回率， $\left| B \right|$为标准数据集中的模块数目.

F度量是精度和召回率的调和平均值，表达式为

(24) $ {\rm{Fm}} = {{2 {\rm{Pr}}\times{\rm{Rc}}}}/({{{\rm{Pr + Rc}}}})\, . $

灵敏度Sn（sensitivity）、正的预测率PPV（positive predictive value）和准确度Acc（accuracy）是另一组评价功能模块检测方法的指标. 表达式分别为

(25) $ {\rm{Sn}} = {{\displaystyle\sum\nolimits_{i = 1}^{\left| B \right|} {\mathop {\max }\limits_j \{ {{T}_{ij}}\} } }}\Bigg/{{\displaystyle\sum\nolimits_{i = 1}^{\left| B \right|} {{N_i}} }}\, , $

(26) $ {\rm{PPV}} = {{\displaystyle\sum\nolimits_{j = 1}^{\left| A \right|} {\mathop {\max }\limits_i \{ {{T}_{ij}}\} } }}\Bigg/{{\displaystyle\sum\nolimits_{j = 1}^{\left| A \right|} {{{T.}_{j}}} }}\, , $

(27) $ {\rm{Acc}} = ( {{\rm{Sn}} \times {\rm{PPV}}})^{1/2} . $

式中： ${{T}_{ij}}$为第i个标准模块 ${b_i}$与第j个检测模块 ${a_j}$所共有的蛋白质数量， ${N_i}$为第i个标准模块 ${b_i}$包含的蛋白质数量， ${{T.}_{j}} = \displaystyle\sum\nolimits_{i = 1}^{\left| B \right|} {{T_{ij}}} $. 准确度利用几何平均数来度量灵敏度与正预测率的综合表现.

对于不同的PPIN数据集，以Cv、Pr、Rc、Fm、Sn、PPV、Acc作为评价指标，通过控制变量法对每一个数据集的参数进行多次调试，从而确定每个数据集的参数. 5个数据集共有的参数确定为：M=100、p=0.62、α=γ=0.95、 $A_i^0$=1、 $r_i^0$=100、 ${f_{\min}}$=0、 ${f_{\max}}$=2、 ${\mu } = 0.95$、Q=20、λ=0.2. 另外，对于ε、δ在不同数据集上的具体设置情况如表2所示.

对所有参数取值方法进行详细描述. 对于M、p、 ${\mu }$，利用最优个体的适应度 ${F_{\max}}$作为度量算法收敛时聚类结果的评价标准. 如图7所示为不同M下 ${F_{\max}}$的取值情况（其余参数保持不变）. 可以看出，当M=100时效果较好，当M>160时，表现也较好，但每一代的搜索时间会随着种群规模的增大而变长，因此，本研究最终取种群规模M=100；同样，p、 ${\mu }$的取值也采用此方法，由于p、 ${\mu }$对 ${F_{\max}}$的影响不明显，不再具体列出. 对于α、γ、 $A_i^0$、 $r_i^0$、 ${f_{\min}}$这5项与标准蝙蝠算法相关的参数，是在文献[15]的基础上进行取值， ${f_{\max}}$则根据自适应变异操作的特点取2. λ决定输出结果中的预测模块数目，如图8所示为在不同λ下6个评价指标的变化趋势（其余参数保持不变）可以看出，Pr、Rc、Fm 3项评价指标的取值对λ的变化不敏感，而随着λ的增加，Sn、Acc不断减小，PPV则有不太明显的增加趋势，因此最终取λ=0.2. ε与数据集中所有节点之间的平均相似度有关，对于平均相似度较小的数据集，ε较小，反之，ε较大. δ的取值与数据集的蛋白质节点总数相关，对于节点数目较多的BioGRID数据集，δ相应较大. 综合以上分析，在所有参数中，M、ε、λ、δ对算法性能的影响较大.

将每个数据集在以上相应参数下独立运行15次取平均值作为最终的实验结果. 为了验证BA-FMD算法的有效性，将其与NACO-FMD、BFO-FMD、Jerarca、COACH、ClusterONE、ClusterEPs这6个经典算法进行对比. 7种算法在5个不同规模数据集上的基本检测结果如表3所示，其中包括检测模块数、检测模块平均大小、 ${N_{a}}$、 ${N_{b}}$.

以DIPcore数据集为例，简单介绍表格的含义，在DIPcore数据集上，BA-FMD算法检测到的功能模块数为416个，其中有163个检测模块与252个标准功能模块相匹配，检测到的功能模块平均包含6个蛋白质. 可以看出，在5个数据集中，BA-FMD算法在 ${N_{a}}$指标上排名前三，在 ${N_{b}}$指标上，在Collins数据集上排名第三，在其余4个数据集上全部取得了第一的成绩. 说明BA-FMD算法检测到的功能模块中有较多的功能模块与标准模块相匹配.

为了进一步证明BA-FMD算法的优越性，如图9所示为7种算法在5个数据集上的7个不同评价指标的对比结果. 综合5个子图的结果可以看出，BA-FMD算法在Gavin、Krogan、Collins、BioGRID这4个数据集上均有4项指标排名第一，在DIPcore数据集上有3项指标排名第一，其余指标也大多排名前三.

1）对于Acc综合指标，BA-FMD算法在5个数据集上均取得了最好的结果. 2）BA-FMD算法在Cv指标上也占有绝对优势，在5个数据集中比排名第二的算法分别高出13.0%、23.0%、12.0%、5.0%、0.5%，说明BA-FMD算法在聚类过程中丢弃的蛋白质较少，能够将大多数蛋白质划分到功能模块中，从而有可能探测到几乎每个蛋白质在生命过程中的功能. 3）BA-FMD算法在Rc指标上表现出色，在Gavin、DIPcore、Krogan、BioGRID这4个数据集中都取得了第一的成绩；在Collins数据集中虽排名第三，但与第一、第二的BFO-FMD、Jerarca较接近. 4）BA-FMD算法在Sn指标上也表现较好，在Gavin、Krogan、Collins 这3个数据集中排名第一；在DIPcore数据集中排名第二，与排名第一的BFO-FMD表现基本齐平；在BioGRID大型数据集中与同类算法BFO-FMD、NACO-FMD相比有明显优势，但不如Jerarca和COACH的表现，这是因为本研究的初始化方法无法识别重叠的功能模块，致使功能模块分布较均匀，而Jerarca和COACH在BioGRID中识别出多个较大的功能模块，与标准模块共有蛋白质的数量较多，Sn指标较高. 5）BA-FMD算法在PPV上表现也不差，在Collins和BioGRID数据集中均取得了最好的结果；在Gavin、DIPcore、Krogan数据集中均略逊于NACO-FMD排名第二，说明PPIN所蕴含的网络特性不同，即使采用同类算法，由于优化机理不同，算法的检测性能也会有所不同.

综合以上5个指标，BA-FMD算法均取得了不错的结果. 原因如下：1）本研究利用动态变化的脉冲频度调节蝙蝠个体的移动方向，使每个蝙蝠个体根据自己当前移动方向情况智能变换下一步的移动策略；2）其次，基于距离和频率的自适应变异操作通过与最优个体位置距离的比较来决定自身位置变化的程度，同时，以一定的可能性重新初始化某些个体，也为种群注入了新鲜血液；3）依据生物界中优胜劣汰的生存法则设计的自然选择操作，大大加快了种群进化的速度，提升了种群的适应能力，使算法获得更优秀的解.

不过，BA-FMD算法在Pr指标上表现不佳，造成该现象的主要原因是BA-FMD算法舍弃的蛋白质节点较少，而在该指标上表现较好的算法大多舍弃了大量的蛋白质，在Pr指标上的表现与Cv指标完全相反. 值得一提的是，结合Pr指标和Rc指标，BA-FMD在4个数据集上均得到了具有竞争力的综合指标Fm，具体来说，该指标在Gavin数据集上以稍低于BFO-FMD算法的结果排名第二；在DIPcore、Collins、BioGRID数据集上排名第三；在Krogan数据集上因Pr指标表现不佳，Fm指标效果也不佳.

通过以上分析可知，BA-FMD算法与其他6种算法相比，在5个不同规模的PPIN数据集上均取得了较好的检测结果，表明BA-FMD算法是有效的功能模块检测方法.

基于蝙蝠觅食原理设计了PPIN功能模块检测的新方法BA-FMD. 该算法最主要的特点是利用标准蝙蝠算法的基本思想，提出了定向局部扰动操作、随机扰动操作、基于距离和频率的自适应变异操作，利用这些优化机制对搜索空间中每一个功能模块划分进行优化；此外，采用较好个体替代较差个体的自然选择操作来提升整个种群的功能模块划分效果. 在5个数据集上进行的实验充分表明BA-FMD算法是有效的检测PPIN功能模块的方法. 该方法的研究，不仅拓展了蝙蝠算法的应用领域，而且对其他复杂网络的分析与研究也有一定的借鉴意义. 本研究所提出的算法未识别出重叠模块，在未来的工作中，将充分考虑重叠模块的检测问题，并融合更多生物信息，获得更高的检测质量，将高质量的功能模块应用到具体的疾病预测中.