基于最近邻区间的不完整基因表达数据多目标聚类算法

常 巧 珍，曹 隽 喆，顾 宏，李 丹

(大连理工大学 控制科学与工程学院，辽宁 大连 116024 )

0 引 言

随着高通量DNA微阵列检测技术的发展，数量庞大的基因相关数据相应而生．基因表达数据反映了直接或间接测量得到的基因转录产物mRNA 在细胞中的丰度[1]，阐明隐藏在这些数据中的模式，从中获取细胞的生理状态、基因表达调控信息以及基因功能，对功能基因组学的研究有着重要的意义．然而，数量庞大的基因和复杂的生物网络成为理解和解释这些数据的巨大挑战，基因聚类能够有效识别共表达基因，推断尚未确定功能基因的表达模式，进而有助于理解基因功能、基因调控及细胞过程[2-3]．

在基因表达数据的获取过程中，受设备、实验环境、采集方法等因素影响，很多数据不可避免地存在缺失值[4]，其填补准确度在一定程度上影响了最终的聚类效果．现有的针对不完备基因表达数据的聚类算法通常为“两阶段”算法[5]，即将缺失值填补作为数据预处理过程，在填补后的数据集上进行聚类，是基因表达数据集聚类分析的常用方法．基因表达数据缺失值预处理的常用方法有：采用均值法计算缺失值对应样本下所有完整表达值的均值作为填补值(Meanimpute)[6]；根据表达值不完整基因的k个完整近邻基因进行缺失值加权估计填补(k-nearest neighbor impute，KNNimpute)[7]；Oba等利用贝叶斯主成分分析法(Bayesian principal component analysis，BPCA)处理基因表达数据中的缺失值[8]；Buuren等则将多重填补法(multivariate imputation by chained equations，MICE)应用于基因表达数据集[9]；Kim等依据皮尔逊相关系数提出了采用多元线性回归模型的局部最小二乘法填补缺失值[10]；Yu等提出了自动估计不同近邻基因权重的自动加权局部最小二乘填补法[11]．除上述几种代表性方法以外，相关文献还利用高斯混合聚类估算法、缺失值多重并行估算法、相关向量机回归估算法等进行基因表达数据缺失值填补．

针对基因表达数据维数高、结构复杂等特点，近年来相关文献提出了基因表达数据的多目标聚类算法．如Bandyopadhyay等提出了一种以度量类内紧密度的Jm[12]和类间分离度的Jxb[13]为目标函数的多目标聚类算法[14]，通过设置不等长编码确定类别数并实现聚类划分；Faceli等将聚类集成思路引入多目标聚类问题，通过初始种群及交叉算子设计实现了多种聚类算法的集成[15]；针对目标函数的自适应选取问题，Mukhopadhyay等提出了多目标交互式聚类算法[16]；Maulik等则提出了将多目标聚类与SVM相结合的算法[17]；为识别形状对称的基因簇，Saha等提出了采用基于对称距离的对称指标及Jxb为目标函数的多目标聚类算法[18]；针对数据集维数高的问题，Liu等提出了利用参考向量划分子空间的多目标聚类算法[19]．上述多目标聚类问题大多以具有低复杂度、高效性及灵活性等特点的NSGA-Ⅱ[20]为多目标优化框架．NSGA-Ⅱ能够在整个解空间内搜索得到一组平衡各目标函数的解，并利用拥挤距离保持种群多样性，广泛应用于基因表达数据的多目标聚类问题．

“两阶段”算法能够实现对缺失基因的后续处理，但其未考虑聚类与缺失值填补的相互影响，造成聚类效果不佳．为此，本文在NSGA-Ⅱ框架下，提出一种基于最近邻区间的不完整基因表达数据多目标聚类算法(multi-objective clustering algorithm based on the nearest neighbor interval，MOC-NNI)．所提算法从近邻相似性角度出发，首先计算缺失值的最近邻区间，进而利用该最近邻区间将缺失值的搜索限定在合理范围内，并在NSGA-Ⅱ框架下实现聚类及缺失值填补的一体化求解，通过两者的协同进化提高缺失值填补准确度及聚类效果．

1 算法介绍

1.1 缺失值的最近邻区间

1.2 目标函数

基因表达数据的聚类问题中，度量聚类结果类内紧密度的Jm和类间分离度的Jxb是常用目标函数[14-17]，本文采用Jm和Jxb作为MOC-NNI的目标函数，Jm和Jxb越小则表示聚类效果越好．

(1)

(2)

式中：uik为基因gi隶属于第k类的程度，∀i，k：uik∈[0，1]，m∈[1，∞)为模糊指数，K为类别数，vk为第k类的聚类中心，D(vk，gi)为基因gi与聚类中心vk的欧几里得距离．uik的计算公式如下：

(3)

缺失值的最近邻区间充分利用基因表达数据中的近邻统计信息，并将缺失值的搜索限定在合理范围内约束聚类中心进化方向，影响基因隶属度，提升聚类效果．

1.3 编码方式及初始种群设置

Eh(t)=(vh11…vh1d…vhK1…vhKdeh1…ehc)

(4)

对于初始种群的设置，个体的eh1，…，ehc可在对应缺失值的最近邻区间内随机生成；聚类中心部分vh11，…，vh1d，…，vhK1，…，vhKd则采用密度峰值法[22]选择K个局部密度高且距其他高密度基因远的基因，将其表达值作为初始聚类中心．上述初始种群设置方法能够将基因表达数据的近邻及密度信息引入初始种群．

在混合编码的基础上实现聚类中心以及缺失值的协同进化，有利于提高NSGA-Ⅱ框架下遗传搜索的收敛速度及优化能力．缺失值的填补值影响目标函数Jm及Jxb的值，进而影响聚类中心进化方向．算法比较目标函数值得到个体非支配排序等级，选择非支配排序等级最高的个体作为子代，在NSGA-Ⅱ框架下进行聚类中心和缺失值的交叉变异，同时利用距离保持种群多样性，利用精英保留策略得到下一代父代种群以及Pareto最优前沿，利用投影相似性指标从前沿中选择最终聚类结果和缺失值的填补值．

1.4 MOC-NNI算法流程

Step 2初始化代数t=0，迭代数Tmax，设定聚类类别数K、种群规模F，选择操作的预定义常数α、交叉算子β、变异概率Pm，在缺失值的最近邻区间中随机生成初始填补值，采用密度峰值法生成个体的初始聚类中心，按照式(4)的混合编码方式产生初始种群．

Step 3由初始种群获得聚类中心以及填补值，根据式(3)得到隶属度矩阵，计算目标函数Jm及Jxb．

Step 4根据Jm及Jxb计算初始种群的拥挤距离和非支配排序等级．

Step 5对第t代种群，采用轮盘赌进行选择[23]，后代竞争择优策略[24]进行交叉，从交叉后代中选择非支配排序等级高且拥挤距离最大的2个个体作为子代．

Step 6对第t代种群，个体中的每个位点以概率Pm发生变异．预处理中对数据进行了max-min归一化，因此变异个体的聚类中心部分为[0，1]内的随机值，填补值部分为相应最近邻区间内的随机值．

Step 7根据子代个体的填补值恢复数据集，依据式(1)、(2)更新子代个体的Jm及Jxb．

Step 8父子代个体融合，依据精英保留策略得到下一代的父代种群以及Pareto最优前沿．

Step 9设置t=t+1，若t

1.5 最终解选取策略

多目标聚类算法终止后，需要从最优前沿Ps中确定最终优化解．采用聚类的内部有效性指标选取最终解不免与算法中的两个目标函数有一定的重合[23]，因此本文采用投影相似性指标[23，25]完成最终解的选取，通过下式度量各类内基因之间的相似性：

(5)

其中nk为划分到第k类的基因数．

(6)

pij=gij×b

(7)

其中b为投影区间分割数．可见，SPSVIndex从基因表达值出发，依据投影坐标衡量基因在各个样本下的表达值的相似性，进而得到基因间相似性，其值越小，表明同一类内的基因越相似．因此本文采用投影相似度指标能够实现从Ps中选取聚类效果最好的解．

1.6 时间复杂度分析

MOC-NNI最坏时间复杂度为O(TmaxFnKd+FN2Kd)，详细分析如下：

(1)目标函数Jm及Jxb的计算时间复杂度均为O(FnKd)．

(2)对于每一次进化操作，交叉和变异操作的时间复杂度分别为O(F(Kd+c))和O(PmF(Kd+c))．

(3)非支配排序时间复杂度为O(2F2)，2为目标函数个数．

(4)从Ps中选取最优解时间复杂度为O(FN2Kd)，N=max(nk)．

K通常远小于n，因此MOC-NNI的时间复杂度由目标函数的复杂度支配，可求得MOC-NNI总迭代数为Tmax的最坏时间复杂度为O(TmaxFnKd+FN2Kd)．

2 结果与讨论

2.1 评价指标

在缺失值填补方面，采用标准化均方根误差E度量填补值与真实表达值之间的偏差：

(8)

(9)

2.2 实验结果与对比分析

实验选取了4个公开的基因表达数据集：拟南芥数据集(Arabidopsis Thaliana)，酵母细胞数据集1(Yeast Cell Cycle_384)，酵母细胞数据集2(Yeast Cell Cycle_237)，人体纤维细胞血清数据集(Serum)．

2.2.1 填补准确度分析 图1所示为各算法在4个基因表达数据集上得到的E值．

可以看出MOC-NNI在各数据集的各种缺失率下均得到了更小的E值，表明所提算法中设计的缺失值与聚类结果的协同优化方法得到了更接近真实表达值的填补结果．相比于MOC-NNI，Meanimpute在填补过程中未考虑数据集中其他基因反映的缺失值分布信息，导致填补效果不理想；KNNimpute的填补结果则易受到k值及权重的影响；BPCA及MICE引入了概率分布模型，通过统计分析在一定程度上提升了缺失值填补效果，但其结果易受到分布模型类型和缺失值不确定性的影响，造成填补效果不佳．MOC-NNI无须引入概率分布模型，充分利用数据集隐含的模式相似性将缺失值的填补限定在一个合理范围内，进而在缺失值与聚类结果的协同优化过程中得到更为准确的填补结果．

(a)Arabidopsis Thaliana数据集

2.2.2 聚类性能分析 表1～4为各算法在基因表达数据集上得到的S值，加粗部分为相同缺失率下的最优值，下划线部分为次优值．可以看出，MOC-NNI除个别情况下取得次优值外均取得了最优S值，表明MOC-NNI中提出的在缺失值最近邻区间约束下进行聚类和缺失值填补协同求解的方法较“两阶段”算法得到了更好的聚类结果，并且MOC-NNI适用于对不同基因表达数据集进行聚类，表现出的鲁棒性较好．结合表1～4以及图1可以看出缺失值填补准确度在一定程度上对聚类结果具有正向影响，与文献[4-5]分析一致．

表1 拟南芥数据集在不同缺失率下的轮廓系数均值Tab.1 Mean values of silhouette index in Arabidopsis Thaliana under different missing rates

图2所示为MOC-NNI在Yeast Cell Cycle_384数据集5%缺失率下得到的聚类热力图及表达谱图．图2(a)中，红色及绿色分别表示高、低表达水平，黑色表示无差异表达值，可见Yeast Cell Cycle_384分成5类，且具有相似颜色排列的基因均被分到同一类中，表明MOC-NNI实现了将表达值相似的基因划分到同一类中．图2(b)中，绿色曲线为每类基因相对于各样本的归一化基因表达值，黑色线条为每类基因的平均表达值与标准差，可以看出同一类内的基因表达谱相似，而不同类的基因表达谱差异较大，表明MOC-NNI对不完整基因表达数据集具有良好的聚类性能．

表2 酵母细胞数据集1在不同缺失率下的轮廓系数均值Tab.2 Mean values of silhouette index in Yeast Cell Cycle_384 under different missing rates

表3 酵母细胞数据集2在不同缺失率下的轮廓系数均值Tab.3 Mean values of silhouette index in Yeast Cell Cycle_237 under different missing rates

表4 人体纤维细胞血清数据集在不同缺失率下的轮廓系数均值Tab.4 Mean values of silhouette index in Serum under different missing rates

(a)聚类热力图

2.3 统计学检验

为了检验MOC-NNI得到的聚类结果是统计显著的，本文还进行了Wilcoxon rank-sum检验．表5所示为4个基因表达数据集在5%缺失率下，MOC-NNI与其他算法所得S值在5%显著性水平下的p值．零假设为MOC-NNI与其他算法所得S值不存在显著差异，备择假设为存在显著差异．

表5 Wilcoxon rank-sum检验所得p值Tab.5 p Values of Wilcoxon rank-sum test

可见，检验所得p值均远小于0.05，表明所提MOC-NNI得到的更优S值在统计学上是显著的，即不是偶然发生的．在其他基因表达数据集的不同缺失率下得到的检验结果类似．

3 结 语

数据不完整问题广泛存在于基因表达数据中，本文从提升缺失值填补准确度出发，提出了一种基于最近邻区间的不完整基因表达数据多目标聚类算法．算法利用最近邻规则挖掘基因表达数据蕴含的统计信息，进而引入最近邻区间描述缺失值的合理搜索范围，在NSGA-Ⅱ框架下通过混合编码实现了缺失值填补与聚类结果的协同进化．在多个基因表达数据集上的实验均表明，本文算法在聚类性能和填补效果方面优于同类算法，能够实现对不完整基因更为可靠的分析及功能推断．