陳佳瑜，李 梁，羅 云

CHEN Jiayu,LI Liang,LUO Yun

重慶理工大學 計算機科學與工程學院，重慶 400054

College of Computer and Engineering,Chongqing University of Technology,Chongqing 400054,China

1 引言

隨著生物技術(shù)的迅速發(fā)展，出現(xiàn)了基因芯片和DNA微陣列等高通量檢測技術(shù)，人們已經(jīng)可以對某一特定狀態(tài)下的組織或細胞進行基因表達譜的測序，而具有相似表達譜的基因之間很有可能存在某種聯(lián)系。目前，基因表達數(shù)據(jù)分析最常用的方法是聚類，將表達水平相近的基因聚在同一個類中，從而可以發(fā)現(xiàn)未知基因的功能和具有生物科學價值的基因信息。但是傳統(tǒng)的聚類在分析基因數(shù)據(jù)時會受到很大的限制，于是Cheng和Church[1]首次提出了雙向聚類的概念，并通過貪婪策略逐步增刪使得雙向聚類的均方殘差達到最優(yōu)。但算法中使用隨機數(shù)替換已找到的雙向聚類，這樣會改變原來的數(shù)據(jù)造成誤差，并且容易陷入局部最優(yōu)。2003年時被提出的FLOC[2]算法先產(chǎn)生一定數(shù)量的種子，再通過增加或刪除行列的方式提高雙向聚類的質(zhì)量。此后又產(chǎn)生了耦合雙向聚類[3]、結(jié)合模糊C均值的雙向聚類算法[4]、貝葉斯雙向聚類模型[5]、并行雙向聚類算法[6]、基于格子模型的雙向聚類[7]、SAMBA算法[8]、使用模擬退火模型尋找雙向聚類[9]等雙向聚類算法。

雙向聚類算法實際上是一種多目標優(yōu)化的局部搜索算法，對于繁雜的基因數(shù)據(jù)矩陣，如果直接對其進行雙向聚類，顯然是不合理的。本文采用的方法是先對整個基因矩陣全局尋優(yōu)，將找到包含最優(yōu)解的子矩陣作為種子，在對其進行雙向聚類。

2 量子粒子群優(yōu)化算法

粒子群（PSO）算法[10]是Kennedy和Eberhart于1995年提出的一種智能優(yōu)化算法，自被提出以來就受到了廣泛的關(guān)注，并在解決實際問題中顯示了其優(yōu)越性。但是在目前已有的PSO算法中，粒子只能沿著特定的軌跡搜索，從而不能保證以概率1收斂到全局最優(yōu)，甚至不能收斂到局部最優(yōu)[11]。文獻[12]根據(jù)人類群體自組織性和協(xié)同性等特點，認為每個個體都可以用量子空間中的一個粒子來描述，建立量子的δ勢阱模型，開發(fā)了一種新的群體智能算法——具有量子行為的粒子群優(yōu)化（Quantum-behaved Particle Swarm Optimization，QPSO）算法。

在QPSO基因聚類問題中，給定一個M×N的基因矩陣A，將M條基因劃分到K個聚類中去。通過粒子在搜索空間中的飛行搜索，選出代表問題解決方案的最優(yōu)粒子。每一個粒子 Xi=(Mi1,…,Mij,…MiK)代表K個聚類的質(zhì)心向量，Mij表示第i個粒子中第 j個聚類的中心向量。Pi(t)為個體最好位置，G(t)為全局最好位置，且G(t)=Pg(t)，g為處于全局最好位置的粒子的下標。

在QPSO算法中，粒子的更新方程為：

QPSO算法不需要粒子的速度信息，β為收斂-發(fā)散因子，是除了群體規(guī)模和迭代次數(shù)外唯一的一個參數(shù)，相對于PSO算法，QPSO具有進化方程簡單，控制參數(shù)少，收速度快，運算簡單等優(yōu)點。

3 雙向聚類與FLOC算法

如圖1所示，雙向聚類拋棄了傳統(tǒng)聚類只從單方向上尋找相似特征，而是同時從行和列兩個方向進行聚類，這樣可以找到隱藏在局部數(shù)據(jù)中的生物信息。

對于一個m行n列的矩陣A，AIJ為A的一個子矩陣，其中I為行（基因）的子集，J為列（條件）的子集，若AIJ的均方殘差[1]值滿足一定的條件，則稱AIJ為一個雙向聚類；其中，均方殘差的計算公式為：

圖1 傳統(tǒng)聚類和雙向聚類

其中，aij為AIJ的有效元素，aiJ、aIj和aIJ分別為行平均值、列平均值和矩陣平均值。

FLOC（Flexible Overlapped Biclustering）算法是一種經(jīng)典的雙向聚類算法。該算法定義了一個收益函數(shù)衡量插入或刪除某行（列）后矩陣質(zhì)量的變化，通過插入或刪除具有最高收益的行（列）的方式來提高雙向聚類的質(zhì)量，進而改進整個雙向聚類種群的質(zhì)量，最終得到最優(yōu)雙向聚類集。它可以在避免隨機干涉的情況下，同時產(chǎn)生多個可能相互重疊的雙向聚類。

FLOC算法的流程為：

步驟1計算基因矩陣各行（列）的收益值。

步驟2對各行（列）的最高收益進行排序。

步驟3按順序插入或者刪除這些行（列）。

步驟4判斷執(zhí)行動作后，雙向聚類的平均容量是否增大，若增大則保存最優(yōu)雙向聚類作為下一次迭代的初始值，回到步驟1繼續(xù)迭代；否則，算法結(jié)束，輸出當前最優(yōu)雙向聚類。

4 基于多樣性選擇的QPSO雙向聚類算法

4.1 基于多樣性選擇的QPSO算法

QPSO算法在搜索開始時，多樣性相對比較高。但是在隨后的搜索過程中，由于粒子的逐漸收斂，群體的多樣性不斷下降，導致大部分粒子都在做局部搜索，而全局搜索能力越來越弱。此時，如果全局最好位置位于局部最優(yōu)區(qū)域或次優(yōu)區(qū)域，算法就會陷入早熟。2002年，Ursem[13]提出了一個多樣性引導的進化算法（DGEA）。同年，Riget等人[14]在PSO中引入Ursem多樣性控制的思想，提出了吸引-排斥粒子群算法（ARPSO）。受其啟發(fā)，提出了一種多樣性選擇的量子粒子群雙向聚類算法（Diversify-Optional QPSO，DOQPSO）。

4.1.1 相似性度量

DOQPSO算法中已知基因表達矩陣A=(X,Y)，包含M條基因和N個樣本條件，在使用歐幾里德距離作為相似性度量的基礎(chǔ)上，將該數(shù)據(jù)集合劃分為指定的K類，使所得到的聚類劃分能使總體類間差異[15]（TWCV）最小。

用Ak表示第k個基因聚類，xin表示在當前聚類中第i個基因在第n個樣本條件下的表達水平；ZK為Ak中質(zhì)心向量的個數(shù)，則 Ak的質(zhì)心向量為Ck=(ck1,ck2,…,ckn)，且。定義聚類Ak中類內(nèi)適應(yīng)度函數(shù)WCV為：

則TWCV可以表示為：

該適應(yīng)度函數(shù)是每種聚類的類內(nèi)適應(yīng)度函數(shù)WCV的和，WCV的值越小說明類內(nèi)各個對象彼此更緊密，聚類質(zhì)量越好。

4.1.2 多樣性度量

ARPSO算法根據(jù)粒子的多樣性，可以在全局搜索和局部搜索之間進行切換，增強了在多峰優(yōu)化問題中尋找全局最優(yōu)解的概率。粒子多樣性度量的公式如下[13-14]：

其中，S表示粒子群，||L表示求解空間最長對角線的長度，N表示問題的維數(shù)。

在DOQPSO算法中，單個粒子的收斂性是受收斂-發(fā)散因子β控制的。粒子進入收斂狀態(tài)，種群的多樣性減小，粒子進入發(fā)散狀態(tài)，種群的多樣性增多。粒子群在初始化后會進入收斂狀態(tài)，這時β會從1.0線性減小到 0.5，如公式（9）所示：

文獻[16]中已證明，以公式（1）為進化方程的QPSO算法中，單個粒子收斂的充分必要條件是β＜1.782。所以存在 β0=1.782，當 β＜β0時，粒子收斂；當 β≥β0時，粒子進入發(fā)散狀態(tài)。

所以算法中另外定義了一個多樣性下限值d_low，如果粒子群的diversity(S)低于d_low，說明粒子的多樣性過少，需進入發(fā)散狀態(tài)，此時重新設(shè)置β的值，且β＞1.782，這樣可以增加粒子的多樣性，直到diversity(S)大于d_low，粒子重新進入收斂狀態(tài)。

算法流程為：

輸入：聚類的個數(shù)K，基因表達數(shù)據(jù)集A，最大迭代次數(shù)GMAX，d_low和微粒數(shù)N。

輸出：最后的聚類質(zhì)心向量，K個聚類。

步驟1初始化每個粒子，使每個微粒中包含K個聚類質(zhì)心。

步驟2初始化粒子群的局部最優(yōu)位置pbest和全局最優(yōu)位置gbest。

步驟3當?shù)螖?shù)小于GMAX時，計算粒子群的平均最好位置C和多樣性函數(shù)diversity(S)。

步驟4粒子群進入收斂狀態(tài)，根據(jù)公式（9）計算β值。

步驟5如果diversity(S)小于d_low，則設(shè)置β＞1.782。

步驟6根據(jù)公式（6）計算適應(yīng)度函數(shù)TWCV；由公式（1）更新粒子的 pbest和種群的gbest。

步驟7回到步驟3，直到達到最大迭代次數(shù)。

步驟8將全局最優(yōu)位置gbest映射為K個聚類質(zhì)心，對數(shù)據(jù)進行劃分，得到K個聚類。

4.2 多目標優(yōu)化的FLOC算法

FLOC算法有兩個目標，一是雙向聚類的剩余值盡可能小，二是雙向聚類的容量盡可能大，這是一個多目標優(yōu)化問題[17-18]。但是FLOC算法判斷進步的依據(jù)僅為容量是否增大，只要有一組的雙向聚類的平均容量比初始化的要大，則認為這一次的執(zhí)行序列執(zhí)行是有進步的，迭代繼續(xù)，這明顯是不合理的。

多目標優(yōu)化算法本質(zhì)上是通過某種方式將復雜的多目標問題轉(zhuǎn)化成容易處理的單目標問題[19]，然后利用單目標優(yōu)化技術(shù)進行求解。常用的多目標優(yōu)化算法有多種求解方法，這里采用評價函數(shù)法[20]。

用評價函數(shù)F(k,A)來計算FLOC算法在一次迭代中產(chǎn)生的雙向聚類的質(zhì)量，F(xiàn)(k,A)的計算公式如下：

F(k,A)表示原基因矩陣A的第k個雙向聚類矩陣的質(zhì)量，rk是它的平均平方殘差，vk是它的容量，λ1和λ2分別是均方殘差和容量的權(quán)值。F(k,A)越小，就表示這個雙向聚類矩陣的質(zhì)量越好，相反，F(xiàn)(k,A)的值越大就表示這個雙向聚類矩陣的質(zhì)量越差。

4.3 基于DOQPSO的FLOC算法思想

算法首先使用DOQPSO算法對原矩陣進行劃分；然后對劃分后的子矩陣進行FLOC聚類，使用公式（10）作為評價函數(shù)。其算法流程如圖2所示。

基于DOQPSO的FLOC算法流程：

步驟1數(shù)據(jù)初始化，使用多樣性選擇的QPSO算法產(chǎn)生k個初始矩陣。

步驟2分別計算這k個雙向聚類矩陣每一行和每一列的動作收益。

圖2 DOQPSO雙向聚類算法執(zhí)行流程圖

步驟3比較收益值的大小，得到各行各列的最優(yōu)動作。

步驟4對這些最優(yōu)動作排序，并按順序依次執(zhí)行。

步驟5執(zhí)行動作后，根據(jù)公式（10）判斷雙向聚類的質(zhì)量是否有所提高，若有，則保存最優(yōu)雙向聚類作為下一次迭代的初始值，回到步驟1繼續(xù)迭代；否則，算法結(jié)束，輸出當前最優(yōu)雙向聚類。

5 實驗及討論

5.1 實驗數(shù)據(jù)與環(huán)境

為了評價雙向聚類的質(zhì)量，采用均方殘差（MSR）、矩陣規(guī)模（V）和行變動值（Var）作為基因矩陣表達一致性的度量標準[1]。目標是找到的雙向聚類均方殘差要盡可能小，規(guī)模和行變動值要盡可能大。行變動值（Var）的計算方法如公式（11）所示：

實驗使用急性白血病基因表達[21]（Leukaemia）數(shù)據(jù)集和酵母菌基因表達數(shù)據(jù)集[22]（yeast）。

急性白血病數(shù)據(jù)集共有38個急性白血病樣本和7 129個基因，其中，27個被診斷為急性淋巴性白血?。ˋLL），11個被診斷為急性骨髓性白血病（AML），所以聚類數(shù)為兩類。

酵母菌數(shù)據(jù)集包含2 884個基因，17個樣本。標準類的劃分是根據(jù)基因表達水平達到峰值的時相進行的，這些基因都是在單個時相達到峰值，細胞周期共分為5個時相，所以將這些基因劃分為5個標準類。

實驗在裝有Win10系統(tǒng)的PC機上運行，CPU為Intel酷睿i5，4 GB內(nèi)存，實驗環(huán)境為MATLAB R2014a。

5.2 實驗結(jié)果及分析

算法將與QPSO、多目標優(yōu)化的FLOC算法，以及未改進的FLOC算法進行比較。粒子群規(guī)模為40，迭代次數(shù)300次；將急性白血病數(shù)據(jù)集聚為兩類，酵母菌數(shù)據(jù)集聚為5類；通過實驗確定公式（10）中參數(shù) λ1取1，λ2取1.8；d_low 設(shè)置為0.000 5[16]。

表1、表2分別是FLOC算法、多目標優(yōu)化的FLOC算法（FLOC2）、QPSO算法和DOQPSO雙向聚類算法在兩個數(shù)據(jù)集上10次實驗的平均值。

表1 四種算法在Leukaemia上的實驗結(jié)果

表2 四種算法在yeast上的實驗結(jié)果

從表中前兩行可以看出，多目標優(yōu)化的FLOC算法的均方殘差值比FLOC算法分別減小了19%和6%，改進后的FLOC算法雖然在規(guī)模上略有減小，但相差不大，所以用了多目標優(yōu)化的FLOC算法表現(xiàn)結(jié)果更佳；從表中后兩行可以看出DOQPSO和QPSO算法的均方殘差值相較于FLOC算法有明顯減小，這是因為算法先對數(shù)據(jù)集進行了劃分，從繁雜的基因數(shù)據(jù)劃分為若干個表達波動相似的塊，剔除了不相關(guān)甚至會產(chǎn)生干擾的數(shù)據(jù)，當然這也導致了在第二階段得到的聚類矩陣的容量有所減少，但是和FLOC算法相比，得到的聚類矩陣相似性更集中，總體上質(zhì)量有很大幅度提高；DOQPSO和QPSO算法相比，雖然規(guī)模相差不大，但是DOQPSO均方殘差明顯減小。說明DOQPSO相較于QPSO算法在性能上有所提高，聚類效果更佳。

6 總結(jié)

基因表達數(shù)據(jù)繁雜且沒有規(guī)律，對其進行聚類分析，得到表達波動行一致的子基因矩陣又是一個多目標優(yōu)化問題，針對這些問題，本文提出了一種DOQPSO雙向聚類算法。算法先利用有多樣性控制的量子粒子群算法處理基因數(shù)據(jù)，根據(jù)多樣性度量函數(shù)的值判斷粒子的狀態(tài)，以提高算法的全局搜索能力；再結(jié)合FLOC算法，修改其迭代目標函數(shù)，彌補了FLOC算法在處理多目標優(yōu)化問題時的不足。經(jīng)過實驗證明，本文提出的算法具有較好的全局尋優(yōu)能力，且聚類效果更佳。

參考文獻：

[1]Cheng Y，Church G M.Biclustering of expression data[C]//Proceeding of the 8th International Conference on Intelligent Systems for Molecular Biology.[S.l.]：ACM Press，2000：93-103.

[2]Yang Jiong，Wang Wei.Enhanced biclustering on gene expression data[C]//Proceedings of the 3rd IEEE Conference on Bioinformatics and Bioengineering，2003：321-327.

[3]Getz G，Levine E，Domany E.Coupled two-way clustering analysis of gene microarray data[C]//Proceedings of the Natural Academy of Sciences of the United States of America，2000，97（22）：12079-12084.

[4]Qu Jinbin，Michael N，Chen Luonan.Constrained subspace clustering for time series gene expression data[C]//4th InternationalConferenceon ComputationalSystems Biology，2010：9-11.

[5]Gu Jiajun，Lee J S.Bayesian biclustering of gene expression data[C]//International Conference on Bioinformatics and Computational Biology，2007：25-28.

[6]劉維，陳崚.基因表達數(shù)據(jù)的并行雙向聚類算法[J].小型微型計算機系統(tǒng)，2009，30（4）：683-689.

[7]Lazzeroni L，Qwen A.Plaid models for gene expression data[J].Statistica Sinica，2002（12）：61-86.

[8]Tanay A，Sharan R，Kupiec M，et al.Revealing modularity and organization in the yeast molecular network by integrated analysis of highly heterogeneous genomewide data[C]//Proceedings of the National Academy of Sciences of the United States of America，2004，101（9）：2981-2986.

[9]Bryan K，Cunningham P，Bolshakova N.Biclustering of expression data using simulated annealing[C]//18th IEEE Symposium on Computer-based Medical Systems，2005：383-388.

[10]Kennedy J，Eberhart R.Particle swarm optimization[C]//IEEE International Conference on Neural Networks，1995，4：1942-1948.

[11]周崸，孫俊，須文波.具有量子行為的協(xié)同粒子群優(yōu)化算法[J].控制與決策，2011，26（4）：582-586.

[12]Sun J，F(xiàn)eng B，Xu W B.Particle swarm optimization with particles having quantum behavior[C]//Proceedings of the Congress on Evolutionary Computation，2004：325-331.

[13]Ursem R K.Diversity-guided evolutionary algorithms[C]//Proceedings of the 7th International Conference on Parallel Problem Solving from Nature.Berlin Heidelberg：Springer，2002：462-471.

[14]Riget J，Vesterstr?m J S.A diversity-guided particle swarm optimizer-the ARPSO，Technical Report no.2002-02[R].EVALife，2002.

[15]高倩倩，須文波.量子行為粒子群算法在基因聚類中的應(yīng)用[J].計算機工程與應(yīng)用，2010，46（21）：152-155.

[16]孫俊.量子行為粒子群優(yōu)化算法研究[D].江蘇無錫：江南大學，2009.

[17]劉文華，梁永全，馮政.基于加權(quán)均方殘差的改進雙聚類算法[J].模式識別與人工智能，2016，29（6）：519-526.

[18]Das C，Maji P.Possibilistic biclustering algorithm for discovering value-coherent overlapping δ-biclusters[J].International Journal of Machine Learning&Cybernetics，2015，6（1）：95-107.

[19]郭紅，蔡莉.采用多目標微分進化算的基因表達數(shù)據(jù)雙向聚類算法[J].小型微型計算機系統(tǒng)，2010（10）：1997-2001.

[20]Cunningham S W，van der Lei T E.Decision-making for new technology：A multi-actor，multi-objective method[J].Technological Forecasting&Social Change，2009，76（1）：26-38.

[21]Golub T R，Slonim D K，Tamayo P，et al.Molecular classification of cancer：Class discovery and class prediction by gene expression monitoring[J].Science，1999，286（12）：205-214.

[22]Cho R J，Campbell M J，Winzeler E A，et al.A genomewide transcriptional analysis of the mitotic cell cycle[J].Molecular Cell，1998，2（1）：65-73