常 巧 珍，曹 雋 喆，顧 宏，李 丹

(大連理工大學(xué) 控制科學(xué)與工程學(xué)院，遼寧 大連 116024 )

0 引 言

隨著高通量DNA微陣列檢測(cè)技術(shù)的發(fā)展，數(shù)量龐大的基因相關(guān)數(shù)據(jù)相應(yīng)而生．基因表達(dá)數(shù)據(jù)反映了直接或間接測(cè)量得到的基因轉(zhuǎn)錄產(chǎn)物mRNA 在細(xì)胞中的豐度[1]，闡明隱藏在這些數(shù)據(jù)中的模式，從中獲取細(xì)胞的生理狀態(tài)、基因表達(dá)調(diào)控信息以及基因功能，對(duì)功能基因組學(xué)的研究有著重要的意義．然而，數(shù)量龐大的基因和復(fù)雜的生物網(wǎng)絡(luò)成為理解和解釋這些數(shù)據(jù)的巨大挑戰(zhàn)，基因聚類(lèi)能夠有效識(shí)別共表達(dá)基因，推斷尚未確定功能基因的表達(dá)模式，進(jìn)而有助于理解基因功能、基因調(diào)控及細(xì)胞過(guò)程[2-3]．

在基因表達(dá)數(shù)據(jù)的獲取過(guò)程中，受設(shè)備、實(shí)驗(yàn)環(huán)境、采集方法等因素影響，很多數(shù)據(jù)不可避免地存在缺失值[4]，其填補(bǔ)準(zhǔn)確度在一定程度上影響了最終的聚類(lèi)效果．現(xiàn)有的針對(duì)不完備基因表達(dá)數(shù)據(jù)的聚類(lèi)算法通常為“兩階段”算法[5]，即將缺失值填補(bǔ)作為數(shù)據(jù)預(yù)處理過(guò)程，在填補(bǔ)后的數(shù)據(jù)集上進(jìn)行聚類(lèi)，是基因表達(dá)數(shù)據(jù)集聚類(lèi)分析的常用方法．基因表達(dá)數(shù)據(jù)缺失值預(yù)處理的常用方法有：采用均值法計(jì)算缺失值對(duì)應(yīng)樣本下所有完整表達(dá)值的均值作為填補(bǔ)值(Meanimpute)[6]；根據(jù)表達(dá)值不完整基因的k個(gè)完整近鄰基因進(jìn)行缺失值加權(quán)估計(jì)填補(bǔ)(k-nearest neighbor impute，KNNimpute)[7]；Oba等利用貝葉斯主成分分析法(Bayesian principal component analysis，BPCA)處理基因表達(dá)數(shù)據(jù)中的缺失值[8]；Buuren等則將多重填補(bǔ)法(multivariate imputation by chained equations，MICE)應(yīng)用于基因表達(dá)數(shù)據(jù)集[9]；Kim等依據(jù)皮爾遜相關(guān)系數(shù)提出了采用多元線性回歸模型的局部最小二乘法填補(bǔ)缺失值[10]；Yu等提出了自動(dòng)估計(jì)不同近鄰基因權(quán)重的自動(dòng)加權(quán)局部最小二乘填補(bǔ)法[11]．除上述幾種代表性方法以外，相關(guān)文獻(xiàn)還利用高斯混合聚類(lèi)估算法、缺失值多重并行估算法、相關(guān)向量機(jī)回歸估算法等進(jìn)行基因表達(dá)數(shù)據(jù)缺失值填補(bǔ)．

針對(duì)基因表達(dá)數(shù)據(jù)維數(shù)高、結(jié)構(gòu)復(fù)雜等特點(diǎn)，近年來(lái)相關(guān)文獻(xiàn)提出了基因表達(dá)數(shù)據(jù)的多目標(biāo)聚類(lèi)算法．如Bandyopadhyay等提出了一種以度量類(lèi)內(nèi)緊密度的Jm[12]和類(lèi)間分離度的Jxb[13]為目標(biāo)函數(shù)的多目標(biāo)聚類(lèi)算法[14]，通過(guò)設(shè)置不等長(zhǎng)編碼確定類(lèi)別數(shù)并實(shí)現(xiàn)聚類(lèi)劃分；Faceli等將聚類(lèi)集成思路引入多目標(biāo)聚類(lèi)問(wèn)題，通過(guò)初始種群及交叉算子設(shè)計(jì)實(shí)現(xiàn)了多種聚類(lèi)算法的集成[15]；針對(duì)目標(biāo)函數(shù)的自適應(yīng)選取問(wèn)題，Mukhopadhyay等提出了多目標(biāo)交互式聚類(lèi)算法[16]；Maulik等則提出了將多目標(biāo)聚類(lèi)與SVM相結(jié)合的算法[17]；為識(shí)別形狀對(duì)稱(chēng)的基因簇，Saha等提出了采用基于對(duì)稱(chēng)距離的對(duì)稱(chēng)指標(biāo)及Jxb為目標(biāo)函數(shù)的多目標(biāo)聚類(lèi)算法[18]；針對(duì)數(shù)據(jù)集維數(shù)高的問(wèn)題，Liu等提出了利用參考向量劃分子空間的多目標(biāo)聚類(lèi)算法[19]．上述多目標(biāo)聚類(lèi)問(wèn)題大多以具有低復(fù)雜度、高效性及靈活性等特點(diǎn)的NSGA-Ⅱ[20]為多目標(biāo)優(yōu)化框架．NSGA-Ⅱ能夠在整個(gè)解空間內(nèi)搜索得到一組平衡各目標(biāo)函數(shù)的解，并利用擁擠距離保持種群多樣性，廣泛應(yīng)用于基因表達(dá)數(shù)據(jù)的多目標(biāo)聚類(lèi)問(wèn)題．

“兩階段”算法能夠?qū)崿F(xiàn)對(duì)缺失基因的后續(xù)處理，但其未考慮聚類(lèi)與缺失值填補(bǔ)的相互影響，造成聚類(lèi)效果不佳．為此，本文在NSGA-Ⅱ框架下，提出一種基于最近鄰區(qū)間的不完整基因表達(dá)數(shù)據(jù)多目標(biāo)聚類(lèi)算法(multi-objective clustering algorithm based on the nearest neighbor interval，MOC-NNI)．所提算法從近鄰相似性角度出發(fā)，首先計(jì)算缺失值的最近鄰區(qū)間，進(jìn)而利用該最近鄰區(qū)間將缺失值的搜索限定在合理范圍內(nèi)，并在NSGA-Ⅱ框架下實(shí)現(xiàn)聚類(lèi)及缺失值填補(bǔ)的一體化求解，通過(guò)兩者的協(xié)同進(jìn)化提高缺失值填補(bǔ)準(zhǔn)確度及聚類(lèi)效果．

1 算法介紹

1.1 缺失值的最近鄰區(qū)間

1.2 目標(biāo)函數(shù)

基因表達(dá)數(shù)據(jù)的聚類(lèi)問(wèn)題中，度量聚類(lèi)結(jié)果類(lèi)內(nèi)緊密度的Jm和類(lèi)間分離度的Jxb是常用目標(biāo)函數(shù)[14-17]，本文采用Jm和Jxb作為MOC-NNI的目標(biāo)函數(shù)，Jm和Jxb越小則表示聚類(lèi)效果越好．

(1)

(2)

式中：uik為基因gi隸屬于第k類(lèi)的程度，?i，k：uik∈[0，1]，m∈[1，∞)為模糊指數(shù)，K為類(lèi)別數(shù)，vk為第k類(lèi)的聚類(lèi)中心，D(vk，gi)為基因gi與聚類(lèi)中心vk的歐幾里得距離．uik的計(jì)算公式如下：

(3)

缺失值的最近鄰區(qū)間充分利用基因表達(dá)數(shù)據(jù)中的近鄰統(tǒng)計(jì)信息，并將缺失值的搜索限定在合理范圍內(nèi)約束聚類(lèi)中心進(jìn)化方向，影響基因隸屬度，提升聚類(lèi)效果．

1.3 編碼方式及初始種群設(shè)置

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

(4)

對(duì)于初始種群的設(shè)置，個(gè)體的eh1，…，ehc可在對(duì)應(yīng)缺失值的最近鄰區(qū)間內(nèi)隨機(jī)生成；聚類(lèi)中心部分vh11，…，vh1d，…，vhK1，…，vhKd則采用密度峰值法[22]選擇K個(gè)局部密度高且距其他高密度基因遠(yuǎn)的基因，將其表達(dá)值作為初始聚類(lèi)中心．上述初始種群設(shè)置方法能夠?qū)⒒虮磉_(dá)數(shù)據(jù)的近鄰及密度信息引入初始種群．

在混合編碼的基礎(chǔ)上實(shí)現(xiàn)聚類(lèi)中心以及缺失值的協(xié)同進(jìn)化，有利于提高NSGA-Ⅱ框架下遺傳搜索的收斂速度及優(yōu)化能力．缺失值的填補(bǔ)值影響目標(biāo)函數(shù)Jm及Jxb的值，進(jìn)而影響聚類(lèi)中心進(jìn)化方向．算法比較目標(biāo)函數(shù)值得到個(gè)體非支配排序等級(jí)，選擇非支配排序等級(jí)最高的個(gè)體作為子代，在NSGA-Ⅱ框架下進(jìn)行聚類(lèi)中心和缺失值的交叉變異，同時(shí)利用距離保持種群多樣性，利用精英保留策略得到下一代父代種群以及Pareto最優(yōu)前沿，利用投影相似性指標(biāo)從前沿中選擇最終聚類(lèi)結(jié)果和缺失值的填補(bǔ)值．

1.4 MOC-NNI算法流程

Step 2初始化代數(shù)t=0，迭代數(shù)Tmax，設(shè)定聚類(lèi)類(lèi)別數(shù)K、種群規(guī)模F，選擇操作的預(yù)定義常數(shù)α、交叉算子β、變異概率Pm，在缺失值的最近鄰區(qū)間中隨機(jī)生成初始填補(bǔ)值，采用密度峰值法生成個(gè)體的初始聚類(lèi)中心，按照式(4)的混合編碼方式產(chǎn)生初始種群．

Step 3由初始種群獲得聚類(lèi)中心以及填補(bǔ)值，根據(jù)式(3)得到隸屬度矩陣，計(jì)算目標(biāo)函數(shù)Jm及Jxb．

Step 4根據(jù)Jm及Jxb計(jì)算初始種群的擁擠距離和非支配排序等級(jí)．

Step 5對(duì)第t代種群，采用輪盤(pán)賭進(jìn)行選擇[23]，后代競(jìng)爭(zhēng)擇優(yōu)策略[24]進(jìn)行交叉，從交叉后代中選擇非支配排序等級(jí)高且擁擠距離最大的2個(gè)個(gè)體作為子代．

Step 6對(duì)第t代種群，個(gè)體中的每個(gè)位點(diǎn)以概率Pm發(fā)生變異．預(yù)處理中對(duì)數(shù)據(jù)進(jìn)行了max-min歸一化，因此變異個(gè)體的聚類(lèi)中心部分為[0，1]內(nèi)的隨機(jī)值，填補(bǔ)值部分為相應(yīng)最近鄰區(qū)間內(nèi)的隨機(jī)值．

Step 7根據(jù)子代個(gè)體的填補(bǔ)值恢復(fù)數(shù)據(jù)集，依據(jù)式(1)、(2)更新子代個(gè)體的Jm及Jxb．

Step 8父子代個(gè)體融合，依據(jù)精英保留策略得到下一代的父代種群以及Pareto最優(yōu)前沿．

Step 9設(shè)置t=t+1，若t

1.5 最終解選取策略

多目標(biāo)聚類(lèi)算法終止后，需要從最優(yōu)前沿Ps中確定最終優(yōu)化解．采用聚類(lèi)的內(nèi)部有效性指標(biāo)選取最終解不免與算法中的兩個(gè)目標(biāo)函數(shù)有一定的重合[23]，因此本文采用投影相似性指標(biāo)[23，25]完成最終解的選取，通過(guò)下式度量各類(lèi)內(nèi)基因之間的相似性：

(5)

其中nk為劃分到第k類(lèi)的基因數(shù)．

(6)

pij=gij×b

(7)

其中b為投影區(qū)間分割數(shù)．可見(jiàn)，SPSVIndex從基因表達(dá)值出發(fā)，依據(jù)投影坐標(biāo)衡量基因在各個(gè)樣本下的表達(dá)值的相似性，進(jìn)而得到基因間相似性，其值越小，表明同一類(lèi)內(nèi)的基因越相似．因此本文采用投影相似度指標(biāo)能夠?qū)崿F(xiàn)從Ps中選取聚類(lèi)效果最好的解．

1.6 時(shí)間復(fù)雜度分析

MOC-NNI最壞時(shí)間復(fù)雜度為O(TmaxFnKd+FN2Kd)，詳細(xì)分析如下：

(1)目標(biāo)函數(shù)Jm及Jxb的計(jì)算時(shí)間復(fù)雜度均為O(FnKd)．

(2)對(duì)于每一次進(jìn)化操作，交叉和變異操作的時(shí)間復(fù)雜度分別為O(F(Kd+c))和O(PmF(Kd+c))．

(3)非支配排序時(shí)間復(fù)雜度為O(2F2)，2為目標(biāo)函數(shù)個(gè)數(shù)．

(4)從Ps中選取最優(yōu)解時(shí)間復(fù)雜度為O(FN2Kd)，N=max(nk)．

K通常遠(yuǎn)小于n，因此MOC-NNI的時(shí)間復(fù)雜度由目標(biāo)函數(shù)的復(fù)雜度支配，可求得MOC-NNI總迭代數(shù)為T(mén)max的最壞時(shí)間復(fù)雜度為O(TmaxFnKd+FN2Kd)．

2 結(jié)果與討論

2.1 評(píng)價(jià)指標(biāo)

在缺失值填補(bǔ)方面，采用標(biāo)準(zhǔn)化均方根誤差E度量填補(bǔ)值與真實(shí)表達(dá)值之間的偏差：

(8)

(9)

2.2 實(shí)驗(yàn)結(jié)果與對(duì)比分析

實(shí)驗(yàn)選取了4個(gè)公開(kāi)的基因表達(dá)數(shù)據(jù)集：擬南芥數(shù)據(jù)集(Arabidopsis Thaliana)，酵母細(xì)胞數(shù)據(jù)集1(Yeast Cell Cycle_384)，酵母細(xì)胞數(shù)據(jù)集2(Yeast Cell Cycle_237)，人體纖維細(xì)胞血清數(shù)據(jù)集(Serum)．

2.2.1 填補(bǔ)準(zhǔn)確度分析 圖1所示為各算法在4個(gè)基因表達(dá)數(shù)據(jù)集上得到的E值．

可以看出MOC-NNI在各數(shù)據(jù)集的各種缺失率下均得到了更小的E值，表明所提算法中設(shè)計(jì)的缺失值與聚類(lèi)結(jié)果的協(xié)同優(yōu)化方法得到了更接近真實(shí)表達(dá)值的填補(bǔ)結(jié)果．相比于MOC-NNI，Meanimpute在填補(bǔ)過(guò)程中未考慮數(shù)據(jù)集中其他基因反映的缺失值分布信息，導(dǎo)致填補(bǔ)效果不理想；KNNimpute的填補(bǔ)結(jié)果則易受到k值及權(quán)重的影響；BPCA及MICE引入了概率分布模型，通過(guò)統(tǒng)計(jì)分析在一定程度上提升了缺失值填補(bǔ)效果，但其結(jié)果易受到分布模型類(lèi)型和缺失值不確定性的影響，造成填補(bǔ)效果不佳．MOC-NNI無(wú)須引入概率分布模型，充分利用數(shù)據(jù)集隱含的模式相似性將缺失值的填補(bǔ)限定在一個(gè)合理范圍內(nèi)，進(jìn)而在缺失值與聚類(lèi)結(jié)果的協(xié)同優(yōu)化過(guò)程中得到更為準(zhǔn)確的填補(bǔ)結(jié)果．

(a)Arabidopsis Thaliana數(shù)據(jù)集

2.2.2 聚類(lèi)性能分析 表1～4為各算法在基因表達(dá)數(shù)據(jù)集上得到的S值，加粗部分為相同缺失率下的最優(yōu)值，下劃線部分為次優(yōu)值．可以看出，MOC-NNI除個(gè)別情況下取得次優(yōu)值外均取得了最優(yōu)S值，表明MOC-NNI中提出的在缺失值最近鄰區(qū)間約束下進(jìn)行聚類(lèi)和缺失值填補(bǔ)協(xié)同求解的方法較“兩階段”算法得到了更好的聚類(lèi)結(jié)果，并且MOC-NNI適用于對(duì)不同基因表達(dá)數(shù)據(jù)集進(jìn)行聚類(lèi)，表現(xiàn)出的魯棒性較好．結(jié)合表1～4以及圖1可以看出缺失值填補(bǔ)準(zhǔn)確度在一定程度上對(duì)聚類(lèi)結(jié)果具有正向影響，與文獻(xiàn)[4-5]分析一致．

表1 擬南芥數(shù)據(jù)集在不同缺失率下的輪廓系數(shù)均值Tab.1 Mean values of silhouette index in Arabidopsis Thaliana under different missing rates

圖2所示為MOC-NNI在Yeast Cell Cycle_384數(shù)據(jù)集5%缺失率下得到的聚類(lèi)熱力圖及表達(dá)譜圖．圖2(a)中，紅色及綠色分別表示高、低表達(dá)水平，黑色表示無(wú)差異表達(dá)值，可見(jiàn)Yeast Cell Cycle_384分成5類(lèi)，且具有相似顏色排列的基因均被分到同一類(lèi)中，表明MOC-NNI實(shí)現(xiàn)了將表達(dá)值相似的基因劃分到同一類(lèi)中．圖2(b)中，綠色曲線為每類(lèi)基因相對(duì)于各樣本的歸一化基因表達(dá)值，黑色線條為每類(lèi)基因的平均表達(dá)值與標(biāo)準(zhǔn)差，可以看出同一類(lèi)內(nèi)的基因表達(dá)譜相似，而不同類(lèi)的基因表達(dá)譜差異較大，表明MOC-NNI對(duì)不完整基因表達(dá)數(shù)據(jù)集具有良好的聚類(lèi)性能．

表2 酵母細(xì)胞數(shù)據(jù)集1在不同缺失率下的輪廓系數(shù)均值Tab.2 Mean values of silhouette index in Yeast Cell Cycle_384 under different missing rates

表3 酵母細(xì)胞數(shù)據(jù)集2在不同缺失率下的輪廓系數(shù)均值Tab.3 Mean values of silhouette index in Yeast Cell Cycle_237 under different missing rates

表4 人體纖維細(xì)胞血清數(shù)據(jù)集在不同缺失率下的輪廓系數(shù)均值Tab.4 Mean values of silhouette index in Serum under different missing rates

(a)聚類(lèi)熱力圖

2.3 統(tǒng)計(jì)學(xué)檢驗(yàn)

為了檢驗(yàn)MOC-NNI得到的聚類(lèi)結(jié)果是統(tǒng)計(jì)顯著的，本文還進(jìn)行了Wilcoxon rank-sum檢驗(yàn)．表5所示為4個(gè)基因表達(dá)數(shù)據(jù)集在5%缺失率下，MOC-NNI與其他算法所得S值在5%顯著性水平下的p值．零假設(shè)為MOC-NNI與其他算法所得S值不存在顯著差異，備擇假設(shè)為存在顯著差異．

表5 Wilcoxon rank-sum檢驗(yàn)所得p值Tab.5 p Values of Wilcoxon rank-sum test

可見(jiàn)，檢驗(yàn)所得p值均遠(yuǎn)小于0.05，表明所提MOC-NNI得到的更優(yōu)S值在統(tǒng)計(jì)學(xué)上是顯著的，即不是偶然發(fā)生的．在其他基因表達(dá)數(shù)據(jù)集的不同缺失率下得到的檢驗(yàn)結(jié)果類(lèi)似．

3 結(jié) 語(yǔ)

數(shù)據(jù)不完整問(wèn)題廣泛存在于基因表達(dá)數(shù)據(jù)中，本文從提升缺失值填補(bǔ)準(zhǔn)確度出發(fā)，提出了一種基于最近鄰區(qū)間的不完整基因表達(dá)數(shù)據(jù)多目標(biāo)聚類(lèi)算法．算法利用最近鄰規(guī)則挖掘基因表達(dá)數(shù)據(jù)蘊(yùn)含的統(tǒng)計(jì)信息，進(jìn)而引入最近鄰區(qū)間描述缺失值的合理搜索范圍，在NSGA-Ⅱ框架下通過(guò)混合編碼實(shí)現(xiàn)了缺失值填補(bǔ)與聚類(lèi)結(jié)果的協(xié)同進(jìn)化．在多個(gè)基因表達(dá)數(shù)據(jù)集上的實(shí)驗(yàn)均表明，本文算法在聚類(lèi)性能和填補(bǔ)效果方面優(yōu)于同類(lèi)算法，能夠?qū)崿F(xiàn)對(duì)不完整基因更為可靠的分析及功能推斷．