双精英遗传策略的基因聚类算法

贾瑞玉，宋飞豹,汤深伟

(安徽大学 计算机科学与技术学院，合肥 230601)

1 引 言

随着DNA微阵列技术的出现，为我们提供了很多的基因表达数据.如何通过基因表达数据来分析出样本的显著性特点是当前的主要研究内容之一.聚类算法是数据挖掘的重要算法[1-3]，也是分析基因表达数据的一个有效方法.王文俊等人[4]提出了核类别保留投影法,将其用于数据特征提取，提高了聚类的效率；黄健恒等人[5]将聚类算法镶嵌于MapReduce框架中，加快了算法的运行速度;刘展杰等人[6]为了提高算法的鲁棒性,利用流形学习将基因表达数据映射到一个局部子空间.虽然聚类算法可以很好的对基因数据进行分析，但传统的聚类算法依然存在初始中心点敏感、局部收敛等问题，常用的解决方法就是使用遗传算法与聚类结合[7].何宏等人[8]通过最大属性划分法来确定初始种群，通过两阶段的交叉变异来提高算法收敛速度.Mitchell M等人[9]通过添加扰动项和补类的策略增加了算法的多样性.

虽然遗传算法能够在一定程度上改善聚类算法的不足，但是遗传算法本身也存在着早熟、局部收敛等问题，为了解决这些问题，孟伟等人[10]通过为每次进化种群引入固定比例的随机个体来降低早熟影响.

近几年来，利用精英策略来改进遗传算法的研究越来越多.慕彩红等人[11]将精英个体的数目提升为M个，通过M个精英个体的协同遗传操作来提高算法的全局求优能力.刘全等人[12]将进化种群以精英集为基础划分成两个子种群,分别进行不同的遗传操作，从而提升了算法的效率.王福才等人[13]提出了一种混合精英策略，兼顾了算法的局部搜索能力和全局搜索能力，黄永青等人[14]提出了基于精英保留策略的蚁群遗传算法,从而提升了算法的效率.

在此基础上，本文提出一种双精英遗传聚类算法(DoubleEliteGeneticClusteringAlgorithm，简称DEGCA)，算法采用保优进化和两个种群协同进化策略，同时又加入了自适应交叉变异率和差异度权重.首先，通过最大属性范围划分得到初始聚类中心的初值，同时采用双种群的遗传算法操作,在两个种群均各自独立运行一定进化代数后，将两个种群进行种群间的交叉操作，使得两个种群能够协同进化，优势互补，完成种群的进化.最后，通过实验验证了该算法的有效性.

2 数据预处理

根据基因表达数据本身的特点，我们在对基因表达数据进行聚类分析的过程中，一般需要进行一定的数据清洗，包括对缺失数据的处理，数据标准化，排除表达不明显的基因数据及数据降维.

2.1 缺失值处理

由于基因表达数据是通过基因芯片技术获取的，在获取数据的过程中可能存在芯片缺陷，像素点不明显等因素而出现部分缺失值.聚类算法是以距离为基础的，由于缺失值的存在使得该基因无法进行聚类分析.所以我们在求解过程中必须对这些缺失数据进行处理.本文利用以下规则对缺失数据进行了处理：

1.当数据缺失太严重时，一般认为超过总数的1/3，则删除该条基因数据；

2.当数据缺失数目在总数1/3之内时，使用平均值填充法填充之.

2.2 数据标准化

在基因表达数据中，由于测量过程中存在系统偏差和标签差异，使得不同基因之间存在数值范围差异.数值范围的差异使得不同基因在聚类中的影响也会有明显的差异.因此数据标准化也是必不可少的操作.本文将基因表达数据映射到0-1范围内，既能减少数值差异带来的影响，又能降低计算复杂度.

2.3 方差分析

基因表达数据通常具有数以万条的基因数据，但并不是每条基因都对基因数据分析做出了等同的贡献.实际上大多数基因在测量中的表达水平变化并不显著，而这些基因的存在大大增加了数据分析的难度.方差分析作为测量数据变异程度的常用手段，我们可以通过方差分析来过滤掉基因数据中数据波动较小的基因.基因在测量中的数据波动越大，表明在实验样本中作用越大.过滤掉数据波动小的基因，一方面可以降低数据维度，减少计算复杂度，另一方面，也可以将富含信息的基因保留下来，使得分析更为准确.

2.4 特征提取

主成分分析(Principal Component Analysis，简称PCA)算法[15]是一种无监督的特征提取方法.利用数据映射把高维数据转化成只含有少数几个综合属性(即主成分)的数据，每个主成分都能够反映出原始数据的大部分信息，并且所含信息互不重复.通过PCA降维的方法既减少了基因表达数据中预测变量的个数，也使得每个主成分的属性相互独立，从而实现降噪的功能.

3 双精英聚类算法

3.1 初始种群的产生

聚类初始中心点的选择对聚类结果的影响是很大的，目前为止最常用的生成初始中心点的方法分别有基于K-means算法的随机生成与基于K-means++算法的最远距离生成，但是两种生成方式并不适用于包含大量冗余数据的数据集，而基因表达数据中含有大量的冗余基因，我们在聚类过程中需要降低这类基因对聚类结果的影响，因此常用的两种初始中心点生成方式并不适用于基因表达数据集.

为了一定程度上降低聚类初始中心点敏感以及冗余基因(即属性)影响的问题，本文提出了一种基于波动最大属性的划分法.

1)根据实验数据集中每个数据的每个属性值建立一个矩阵Un*l，其中n为数据个数，l为属性个数，uij为第i个数据在第j个属性上的值.

2)按公式(1)计算数据集中每个属性的方差.

(1)

3)根据每个属性的方差降序排序,选取前m个属性(m与n相关)，使用自底而上的层次聚类算法对具有m个属性的数据集进行聚类；层次聚类过程如下：

①令每个数据为一个簇类.

②根据当前数据集中簇群的个数k重新计算并建立一个距离矩阵Dk*k.

(2)

公式(2)中其中I，J分别为一个簇类集合，dis(I,J)为簇类I和J的欧式距离.

(3)

③寻找矩阵中距离最小且不为0的元素，即可得到距离最短的两个不同簇类I和J.如果dIJ小于事先设定的阈值thre，则合并两个簇类I和J，并返回步骤②；否则转向4).

4)依次从层次聚类得到的k个簇类中随机挑选一个数据点组成一个初始中心点.

5)重复4)过程N次，由此得到初始种群，其中N为初始种群规模.

选取波动最大的m个属性是为了降低冗余属性对聚类结果的影响，同时也降低了算法过程中的计算量；我们使用层次聚类将数据集划分为不同簇类，那是因为层次聚类算法不受初始中心的影响且不需要事先设定簇类个数.

3.2 适应度函数

设数据集D=(x1,x2,…,xn)共有k个数据中心，其k个中心点为C=(c1,c2,…,ck).则聚类所得的聚类总代价如下：

(4)

其中CLi是根据中心点ci划分的簇，o是簇CLi中的元素，dist为欧式距离.E代表了数据集在中心点C=(c1,c2,…,ck)作用下所求出的聚类代价，其值越小，求得的簇内距离越小.为了体现聚类代价E的作用，本文所采用的适应度函数如下：

(5)

其中Exi代表数据集A在进化个体Xi下的聚类代价.f(Xi)的值越大，聚类代价越小，则聚类的效果越好.f(Xi)的最大值即为本文需要求得的最优聚类方案.

3.3 双精英进化机制

精英个体对于种群的进化具有重要的引导作用，在精英策略下种群可以较快地收敛到全局最优.DEGCA是双种群的进化模式，为了充分发挥精英策略的作用，本文采用了精英库策略.

精英库共有两个精英个体，开始是由两个初始种群种群A和种群B中的最优个体EliteA和EliteB组成的.EliteA和EliteB分别引导种群A和种群B各自进化一定代数后，进行种群间的交叉，同时在精英库中更新EliteA和EliteB的值.然后，种群A和种群B重新从精英库中选择精英个体引导进化，直至完成整个进化过程.

3.3.1 种群A的进化

种群A的精英个体EliteA初始值是从种群自身产生的，之后在种群内的进化过程中随着遗传进化不断地自我更新，在进行种群间的交叉操作之后，再将EliteA与精英库中适应度较差的个体进行对比，若EliteA较大，则将EliteA替换掉精英库中适应度较差的个体.然后在精英库中挑选适应度最高的个体作为ELiteA继续引导种群A的进化.

种群A的进化目的是种群的收敛速度与多样性的平衡.

3.3.1.1 自适应交叉变异率调整

为了加快种群收敛速度，种族A在进化中的交叉率和变异率将随着进化个体的适应度而不断调整，使得适应度较高的个体具有较低的交叉变异率，从而可以保护种群中的优秀个体.而种群中适应度较低的个体具有更高的交叉变异率，从而使得适应度低的个体产生更多的变化.其交叉率Pc及变异率Pm为：

(6)

(7)

其中，fmax为种群中的最优个体的适应度，favg是种群的平均适应度，f′是交叉个体中适应度较高的个体，f为将要变异个体的适应度，k1,k2,k3,k4是种群进化的可控参数.由公式(6)、公式(7)可得，当进化个体适应度越差时，交叉率变异率就越大，越容易产生优秀的个体，而当进化个体的适应度较高时，交叉率变异率较低，从而可以减少对优秀个体的破坏.

传统的遗传算法常常会出现优秀个体经过交叉变异之后，得到的子代反而处于一个较低的水平.为了使进化个体的交叉变异尽量朝着一个正确的方向进行，同时使个体中优秀基因能够保留下来.本文提出一种保优进化的策略来保障优秀个体的权益.

定义 1.(保优进化)：当优秀个体(适应度高于种群的平均适应度)经过交叉变异后的子代个体，其适应度低于种群平均适应度时，保留父代个体作为新种群的成员，同时该父代个体将不再参与后续的进化；若子代个体的适应度高于种群平均值，则保留子代，父代个体依然可以参与后续的交叉变异中.

通过保优进化的策略可以最大程度的发挥优秀个体中优良基因的作用，而抑制了较差基因的作用，从而加快种群的收敛速度.

3.3.1.2 自适应差异度调整

为了保证种群在加快收敛过程中的多样性，本文提出了基于差异度的交叉个体选择策略，即根据与目标个体的差异度选择交叉个体，算法过程中根据种群的不同进化程度自适应调整差异度权重kd.即当种群中适应度大于平均适应度的个体超过2/3，说明该种群整体较优，因此若个体xi的适应度大于平均适应度，我们应该为其挑选差异度较大的个体；当种群中适应度大于平均值的个体少于1/3，说明该种群整体较差，因此对于适应度大于平均值的个体xi，我们应该为其挑选差异度较小的个体.根据这种思想，我们提出了差异度权重调整方式：

(8)

其中kd(j)为第j个个体的差异度权重，fj为已选择的交叉个体(第j个)适应度，favg为种群的平均适应度，sumgood为适应度大于平均值的个体数目，sum为种群规模.

得到差异度权重kd之后，我们就可以计算个体ai与aj之间的差异度D(i,j)：

(9)

L为遗传算法中二进制编码长度，ail为个体ai在第l位的编码，取值为0或1.

在交叉操作中，在选择第一个交叉个体ai之后，第二个个体xj被选择的概率为：

F(xj)=D(i,j)*f(xj)

(10)

其中f(xj)为个体xj的适应度，F(xj)为个体xj被选择的概率；适应度越大，D(i,j)越大，个体被选择的概率就越大，这有效地避免了无效交叉(即个体被交叉的基因完全相同).

3.3.2 种群B的进化

种群B是以精英个体EliteB为主导的，EliteB的产生方式与ELiteA相似，只是每次从基因库中要选择与EliteA相异的个体.种群B的进化目的是为了扩大种群的多样性，避免早熟现象，为了实现这一目的，DEGCA引入一定数目随机个体来替换掉每次遗传操作后种群B的进化团队中适应度较差的个体.为了使随机个体的数目与种群整体优秀程度相关，对引入随机个体的数目rnum由公式(11)和公式(12)来确认：

(11)

rnum=(τ0×PA+τ1)·N

(12)

其中，fbest是种群A经历过的最优个体的适应度，fmin是种群A当前的最差个体的适应度，N为种群个体数目.PA表征种群A的进化程度，PA越大，种群A的整体适应度则越小，也需要更多的随机个体为种群产生更多的变化.公式(11)、(12)保证了随机个体的数目将随着种群的整体适应度而不断调整，且在一定范围之内.

3.4 种群间的交叉

当种群A和种群B各自独立进化r代后，分别从种群A和种群B中随机选择m个个体进行交换，这一步骤相当于交换种群A和种群B中的m个个体，使得种群A和种群B得以优势互补，协同进化.交叉之后的两个种群各自独立的进行进化，直至满足终止条件为止.

3.5 DEGCA算法描述

输入：数据集A，中心点个数k，种群内独立进化代数r

输出：聚类中心点集合

算法步骤如下：

1)数据预处理

2)初始化，使用波动最大属性划分法产生数目相同的两个初始种群.并根据初始种群来初始化精英库

3)计算种群中每个个体的适应度

①使用k-means方法更新质心

②使用公式(5)计算个体的适应度

③更新精英库，若两个种群中的最优个体适应度高于精英库的最差个体，则替换之.更新精英库后，将精英库中的数据作为EliteA和EliteB

4)种群A的进化

①使用轮盘赌选择法挑选进化个体

②使用公式(6)-公式(10)计算交叉变异率并选取交叉个体进行交叉操作

③对交叉变异后的个体，采用保优进化的策略进行筛选

5)种群B的进化

①使用公式(12)获取rnum个随机个体并替换掉种群B中适应度最差的rnum个个体，得到进化团队TempB

②将TeamB与EliteB一同完成进化

6)种群间的交互，当种群A和种群B各自独立进化r代后，分别从种群A和种群B中随机选择m个个体进行交换，同时更新精英库

7)判断是否满足终止条件，满足则输出聚类结果，否则返回步骤3)继续进化

4 实验仿真

4.1 评价标准

本文使用ARI(Adjust Rand index)指标，NMI(Normalized Mutual Information)指标和聚类纯度来作为评判聚类效果的评价标准.聚类纯度(purity)能够准确反映聚类结果与真实划分的符合程度[17].假设真实划分为S={s1,s2,…,sk}，聚类得到的簇为Sc={Sc1,Sc2,…,Sck}，则聚簇Sci的纯度可由公式(13)表示：

(13)

其中，|S′|为聚簇Sci中包含真实划分类Si的个数，其最大值就是该聚簇占比最多的类.若聚类生成的某一簇与真实划分的某一类纯度较高，而与其它类的较低，说明其与该类的符合程度较高.所有聚簇划分的总体纯度为：

(14)

其中，k为划分簇的个数，Pt的取值范围在0和1之间.Pt值越大则表示聚类的结果越符合真实划分.ARI和NMI为聚类中常用的评价标准，这里不作详细介绍.

4.2 实验结果

为了验证DEGCA算法的有效性，本文将DEGCA与K-means，Kmeans++，以及标准遗传聚类算法(Standard Genetic Clustering Algorithm，简称SGCA)进行对比，从基因表达综合数据库(Gene Expression Omnibus，简称GEO)上选取两个经典数据集GDS5401和GDS5403进行测试.其中GDS5401出版于2014年2月，是Woetzel等人[16]在柏林测量正常人(ormal)，类风湿关节炎患者(Rheumatoid)和骨关节炎患者(Osteoarthritis)的滑膜组织基因数据，用以区分正常人，类风湿关节炎患者和骨关节患者.2014年3月，Woetzel等人又发布了耶拿人群的正常，类风湿关节炎患者和骨关节炎患者的滑膜组织基因数据，即GDS5403.具体数据如表1所示.

表1 数据集的描述

Table 1 Description of data set

数据集属性值正常类风湿骨关节炎GDS540122283101010GDS540322283101310

4种算法各自运行10次，得到的NMI和ARI如表2和表3所示.从表中可得，K-means++在两个数据集的ARI和NMI均高于K-means，这说明K-means++在一定程度上解决了K-means初始中心点敏感的问题.SGCA聚类得到的ARI和NMI均优于k-means，这说明了遗传算法在一定程度上提高了聚类的精度，在GDS5401中SGCA的性能优于K-means++，而在GDS5403中SGCA的ARI指标高于K-means++，而NMI要略低K-means++，这说明SGCA相对于K-means++没有明显的优势.DEGCA在AMI和ARI指标上均优于其他三个聚类算法，这说明DEGCA聚类结果更为准确.

表2 各个算法求得的ARI

Table 2 ARI obtained by each algorithm

数据集函数均值最大值最小值GDS5401k-means0.5230.6010.423k-means++0.5840.6670.553SGCA0.6370.8020.517DEGCA0.7830.8980.667GDS5403k-means0.5510.6010.462k-means++0.5810.6030.575SGCA0.5780.6010.545DEGCA0.6520.6670.601

表3 各个算法求得的NMI

Table 3 NMI obtained by each algorithm

数据集函数均值最大值最小值GDS5401k-means0.5930.6850.452k-means++0.6570.7460.632SGCA0.6830.7990.585DEGCA0.8220.8980.746GDS5403k-means0.5960.6350.534k-means++0.6130.6840.575SGCA0.6090.6350.604DEGCA0.7010.7190.634

4个算法在这10次实验中的纯度如表4所示，相对于k-means，在两个数据集中k-means++的平均纯度分别提升了0.064和0.067，SGCA分别提升了0.102和0.046，而DEGCA提升幅度最大，为0.159和0.121，这说明DEGCA的聚类结果最接近真实划分的结果，性能提升也最为明显.

表4 各个算法求得的纯度

Table 4 Purity obtained by each algorithm

数据集函数均值最大值最小值GDS5401k-means0.7620.8330.700k-means++0.8260.8670.733SGCA0.8640.9330.767DEGCA0.9210.9670.867GDS5403k-means0.7360.8180.727k-means++0.8030.8480.788SGCA0.7820.8180.758DEGCA0.8570.8790.818

4.3 多样性度量

文献[12]提出了在以二进制编码为基础的遗传算法中，如何度量种群多样性的计算公式.

(15)

其中L为个体的编码长度，∑alj0表示所有个体在第l位为0的个数，∑alj1表示所有个体在第l位为1的个数，max(∑alj0,∑alj1)-min(∑alj0,∑alj1)代表了种群0和1的分布情况，其值越大则表示种群中0和1的分布越不均衡，从而种群的多样性越低.

图1 种群多样性在GDS5401的变化

图2 种群多样性在GDS5403的变化

SGCA和DEGCA求得种群多样性变化图如图1和图2所示.从图1和图2可知SGCA运行到一定代数后，种群多样性将会下降到一个较低水平，从而使算法陷入早熟收敛，而DEGCA由于种群B引入了部分随机个体，种群多样性一直维持在一定的水平，从而减少了早熟现象.

4.4 DEGCA与SGCA收敛速度对比

SGCA设置最大运行代数为100，种群数目为200，DEGCA种群内进化代数为5，种群间交叉最多20次.种群A和种群B的种群规模均为100.实验结果见表5.

表5 两种算法收敛的平均代数和运行时间

Table 5 Average algebra and running time of convergence of the two algorithms

数据集SGCADEGCA种群A种群B总种群 收敛代数46283430GDS5401每代耗时(s)252.18131.17127.25258.42 收敛耗时(s)115963672.84326.57752.6 收敛代数53323735GDS5403每代耗时(s)284.52148.96146.74295.7 收敛耗时(s)150804766.75429.410350

表5表明，在算法收敛代数上，DEGCA比SGCA缩短了将近1/3.在DEGCA的双种群中，种群A的收敛速度较快，而在种群A的带动下，种群B也能够快速达到收敛.在算法效率上看，DEGCA在每代耗时上不如SGCA，但是在收敛耗时上，GDS5401提高了33.1%；GDS5403提高了31.4%.

5 结 论

本文提出一种基于双精英遗传策略的基因聚类算法，在传统遗传算法基础上采用两个种群协同保优进化策略，同时又提出了自适应调整的交叉变异率和差异度选择交叉，有效地提高了算法的收敛速度.相对于k-means,k-means++和SGCA，DEGCA聚类结果更接近真实划分，聚类更加准确.由于算法中加入太多进化策略，导致运行时间也相对较长，因此我们下一步的工作旨在如何在保证聚类结果的同时降低算法的时间复杂度.