崔衍，薛源

（南京郵電大學物聯(lián)網(wǎng)學院，南京 210003）

0 引言

各種生物的生命活動都是由基因調(diào)控的，近年來DNA微陣列技術(shù)憑借其優(yōu)點成為研究基因的主要技術(shù)［1］。DNA微陣列技術(shù)產(chǎn)生的基因表達數(shù)據(jù)是一個矩陣，反映的是基因轉(zhuǎn)錄產(chǎn)物mRNA在細胞中的豐度，通過分析這些數(shù)據(jù)能夠得到當前細胞的生理活動狀態(tài)。該矩陣的行代表的是基因，列代表的是條件（實驗環(huán)境）。對基因表達數(shù)據(jù)進行分析的傳統(tǒng)方法是聚類，但是傳統(tǒng)的聚類只能對基因進行聚類或者對條件進行聚類，這樣得到的是部分基因在所有條件下相似的表達模式或者所有基因在部分條件下相似的表達模式［2］，并且得到的聚類元素集合沒有相交，即一個基因或者條件只能在一個聚類里，不能屬于多個聚類［3］，這并不能很好的應用于對基因表達數(shù)據(jù)進行分析，因為一部分基因可能在多個條件下調(diào)節(jié)同一生物功能［4］。Cheng和Church［5］在2000年首次將雙聚類應用于基因表達數(shù)據(jù)分析。雙聚類顧名思義，就是同時對行（基因）和列（條件）進行聚類，得到一個具有相似表達模式的子矩陣。雙聚類已經(jīng)被證明是一個NP-hard問題，智能優(yōu)化算法在解決NP-hard問題方面有著獨特的優(yōu)勢，因此將智能優(yōu)化算法應用于雙聚類算法是一個很好的研究方向。

Yang［6］在2010年根據(jù)蝙蝠回聲定位提出了蝙蝠算法（bat algorithm，BA）。BA算法在解決復雜的優(yōu)化問題方面有著很好的性能，引起了學者的極大關(guān)注，現(xiàn)在已經(jīng)應用于多個領(lǐng)域來解決問題。本文提出了一個基于改進蝙蝠算法的雙聚類算法。改進了原始蝙蝠算法的位置和速度更新公式，并引入變異算子來提高局部搜索能力。最后在酵母菌數(shù)據(jù)集上進行實驗，與多個雙聚類算法的實驗結(jié)果進行比較。結(jié)果表明本文提出的基于蝙蝠算法設計雙聚類算法是一個可行的方法，能夠挖掘出更優(yōu)的雙聚類。

1 雙聚類定義

Cheng和Church提出CC算法的時候給出了雙聚類的定義［5］，定義如下：

定義1設M行N列的基因表達矩陣為A，G為包含M個基因的基因集，C為包含N個條件的條件合，a i j為基因表達矩陣A的一個元素。雙聚類定義為B=(I,J)，其中I為G的一個子集，J為C的一個子集，b ij為子矩陣B的元素，子矩陣B的平均平方殘差為：

其中：a iJ、a Ij、a IJ分別為子矩陣B的第i行平均值、子矩陣B的第j列平均值和子矩陣B(I,J)的平均值。存在一個δ≥0，如果子矩陣B(I,J)滿足H(I,J)≤δ，則稱該子矩陣為一個δ-bicluster。

雙聚類的目的就是在基因表達數(shù)據(jù)矩陣中尋找滿足條件的子矩陣，使得子矩陣中基因集在對應的條件集上具有相似的表達模式。

2 蝙蝠算法

蝙蝠算法是基于蝙蝠的回聲定位行為提出來的［7］。蝙蝠的回聲定位能夠幫助蝙蝠探測到目標的方向和位置，還可以幫助蝙蝠躲避障礙物。蝙蝠算法中將蝙蝠的一些行為理想化，理想化的規(guī)則如下：

（1）所有的蝙蝠都利用回聲定位來感知自身與目標的距離，并且以某種特殊的方式區(qū)別目標與障礙物。

（2）每只蝙蝠都擁有相同的脈沖頻率f、不同的波長λ和不同的響度A，在位置x以速度v隨機飛行。它們根據(jù)目標的接近程度調(diào)整發(fā)射脈沖的波長λ（或頻率f），并調(diào)整脈沖發(fā)射頻率r。

（3）假設響度從最大值A0變化到Amin。

該算法的主要思想是基于以上理想化規(guī)則，改變脈沖發(fā)射頻率、脈沖頻率、聲音響度，來尋找最優(yōu)解［8］。定義蝙蝠的頻率、速度、位置更新公式：

其中f i為第i只蝙蝠的脈沖頻率，fmax和fmin分別為脈沖頻率的最大值和最小值，β是一個服從均勻分布的隨機值且β∈[0,1]，和分別是第i只蝙蝠在第t代和t-1代的飛行速度；和分別為第t代和第t-1代的位置；x*為當前群體中蝙蝠的最優(yōu)位置（最優(yōu)解）。

為了改善算法的局部搜索能力，使用如下公式更新進行局部搜索：

其中，At為第t次迭代所有蝙蝠的平均響度，ε∈[-1,1]是一個隨機數(shù)。隨著迭代的進行，響度A i和發(fā)射脈沖率r i都會進行更新，更新公式如下：

其中，α和γ是常數(shù)，是第i只蝙蝠的初始脈沖發(fā)射率。

蝙蝠算法的偽代碼如下所示：

算法1原始蝙蝠算法

目標函數(shù)：f(x),x=(x1,x2,…,x d)T;

1.初始化蝙蝠種群x i，速度v i，蝙蝠x i的頻率f i，脈沖發(fā)射率r i，響度A0，i=1,2,…,n

2.初始化參數(shù)蝙蝠種群數(shù)n,α,γ和最大迭代數(shù)ge nmax

3.根據(jù)目標函數(shù)找到當前全局最優(yōu)解

4.fori=1 togenmax:

5.根據(jù)公式（4），（5），（6）生成新解

6.ifrand>r i（rand是一個隨機數(shù)，且rand∈[0,1]）

7.根據(jù)公式（7）生成一個新解

8.end if

9.通過目標函數(shù)評價新解

10.ifrand

11.將7生成的解作為新解；

12.根據(jù)公式（8）,（9）更新脈沖發(fā)射率r i和響度A i

13.end if

14.更新全局最優(yōu)解x*

15.end for

3 基于改進蝙蝠算法的雙聚類算法

在這一部分我們提出基于改進蝙蝠算法的雙聚類算法。原始的蝙蝠算法是適用于解決連續(xù)問題，不適用于解決離散問題，顯然原始蝙蝠算法并不適用于雙聚類算法，因此對蝙蝠算法做出如下改進：

（1）采用V型傳遞函數(shù)對位置更新公式作出改變，使蝙蝠算法適用于解決離散問題［9］。

（2）采用動態(tài)遞減慣性權(quán)重和速度調(diào)整因子，來增強全局搜索和局部搜索［10］。

（3）擴展遺傳算法，采用一種變異算子進行局部搜索［11］。

3.1 雙聚類編碼方式

對雙聚類進行編碼，使之能夠表示蝙蝠算法的種群中的一個個體。編碼方式采用文獻［12］中的編碼方式。規(guī)則如下：用長度為N+M的二進制字符串進行編碼，其中N和M分別代表基因表達式矩陣的行（基因）數(shù)和列（條件）數(shù)。二進制字符串的前N位表示基因，后M位表示條件。若某位為1，則說明該位對應的基因或者條件屬于雙聚類，否則，不屬于雙聚類。如一個基因表達矩陣大小為7×6，包含7個基因，6個條件，該表達矩陣中的一個雙聚類編碼為1001001010010，表明該雙聚類含有3個基因和2個條件，大小為6，基因分別為基因表達矩陣中的第1、4、7個基因，條件分別為基因表達矩陣中的第2、5個條件。

3.2 初始化種群

蝙蝠算法是一種基于種群的進化算法，因此需要對種群進行初始化。采用文獻［13］的方法對種群進行初始化。首先用K-means算法分別對基因和條件進行聚類。共生成a個基因簇和b個條件簇，再將基因簇和條件簇進行組合，得到a×b個雙聚類，對這些雙聚類進行編碼，以此作為初始種群。

3.3 速度更新

我們使用動態(tài)遞減慣性權(quán)重［10］來更新速度。在原始BA算法的速度更新公式中，上一代速度v t-1i的系數(shù)為1，不利于進行全局搜索，容易陷入局部最優(yōu)。因此引入動態(tài)遞減慣性權(quán)重，首先蝙蝠擁有較大的慣性權(quán)重，使蝙蝠的速度能夠在較大的范圍變化，有助于進行全局搜索。隨著迭代的進行，權(quán)重逐漸變小，使蝙蝠的速度在較小的范圍變化，有助于進行局部搜索。公式如下：

其中，ωmax是最大權(quán)重，ωmin是最小權(quán)重，t代表當前迭代次數(shù)，genmax代表最大迭代次數(shù)。arctan函數(shù)是單調(diào)遞增函數(shù)，所以ωi(t)是單調(diào)遞減函數(shù)，因此前期權(quán)重大，后期權(quán)重小。使算法得到前期速度快后期速度慢。前期快速全局搜索，后期慢速局部搜索，有效提高算法的速度和精度。

3.4 位置更新

在原始蝙蝠算法中，蝙蝠的速度負責更行蝙蝠的坐標，且兩者都是連續(xù)值。由上文可知將蝙蝠算法應用于雙聚類算法中時，每個蝙蝠的坐標是由一串二進制數(shù)字表示的，因此需要引入傳遞函數(shù)來負責由速度值變換到坐標值。本文采用V型傳遞函數(shù)［9］來更新坐標，如方程。

其中x ti(k)表示第i個蝙蝠在第t次迭代的第k維的坐標值，v ti(k)表示第i個蝙蝠在第t次迭代的第k維的速度值，(x ti(k))-1是x ti(k)取反得到的。

3.5 局部搜索

對于局部搜索部分，原始蝙蝠算法使用隨機游走公式（4）為每個蝙蝠生成一個新的解。然而，本文提出的算法的搜索空間是離散空間，因此需要使用另一種方法來進行局部搜索。本文使用位翻轉(zhuǎn)變異算子［11］。變異算子如下：

其中rand是一個隨機值且rand∈[0,1]，Pm是一個變異概率，x*全局最優(yōu)解。

3.6 適應度函數(shù)

本文提出的算法的主要目標是找到具有較低均方殘差的同時還要兼顧具有較大體積的雙聚類，適應度函數(shù)如下：

其中ω1和ω2分別為雙聚類的均方殘差和雙聚類的體積的權(quán)值。均方殘差residue(B)的計算如公式（1），s ize(A)為基因表達矩陣的大小，計算公式如下：

size(B)為雙聚類B的體積，計算公式如下：

3.7 基于蝙蝠算法的雙聚類算法實現(xiàn)

基于蝙蝠算法的雙聚類算法的實現(xiàn)過程如下：首先，用3.2節(jié)描述的種群初始化方法初始化種群，再依據(jù)3.1節(jié)描述的編碼規(guī)則對每個蝙蝠個體進行編碼；其次，以3.4節(jié)、3.5節(jié)和3.6節(jié)描述的速度更新、位置更新、局部搜索綜合而成的搜索方法對解空間進行搜索尋找質(zhì)量更高的雙聚類；最后重復前兩個步驟，直到算法完成規(guī)定的迭代次數(shù)，輸出雙聚類。該算法的偽代碼如算法2所示。

算法2 BIBA

輸入：基因表達矩陣，蝙蝠種群數(shù)n，MSR閾值δ，最大迭代次數(shù)genmax，雙聚類數(shù)量BiNum,變異概率Pm

輸出：雙聚類集合

1.fornum=1 toBi Num:

2.初始化響度A i，脈沖發(fā)射率r i，頻率f i，速度v i，i=1,2,…,n

3.通過K-means算法分別對矩陣行和列進行聚類，得到行聚類和列聚類

4.將行聚類和列聚類相乘，得到一定個數(shù)的雙聚類，選取M SR值最小的前n個雙聚類進行編碼，將其作為初始種群x i

5.for種群迭代次數(shù)iter=1 to max_gen

6.fori=1 ton

7.根據(jù)公式（4），（11），（10），（12），（13）分別更新頻率，速度權(quán)重，速度，坐標

8.ifrand>r i（rand是一個隨機數(shù)，且rand∈[0,1]）

10.end if

12.ifrand

13.接受9產(chǎn)生的新解，然后按照公式（8）,（9）更新響度A i和脈沖發(fā)射率r i

14.end if

15.end for

16.更新全局最優(yōu)解x*

17.end for

18.將迭代得到的雙聚類存入雙聚類結(jié)果集中

19.end for

20.輸出結(jié)果

4 結(jié)果與分析

4.1 數(shù)據(jù)集

實驗所使用的數(shù)據(jù)集是釀酒酵母菌數(shù)據(jù)集［14］。數(shù)據(jù)集由2884個基因和17個條件組成，所有基因在所有條件下的表達值都位于0~600之間，其中有34個缺失值，用0~800之間的隨機數(shù)進行填充。

4.2 參數(shù)設置

初始化參數(shù)設置如下：蝙蝠數(shù)量n=100，響度A=0.6，脈沖發(fā)射率r=0.5,α=γ=0.9，頻率fmin=0,fmax=1，最大種群迭代次數(shù)t=2500，變異算子，遞減慣性權(quán)重ωmax=0.9,ωmin=0.42,ω1=0.95,ω2=0.05。

4.3 實驗結(jié)果

應用改進蝙蝠算法共產(chǎn)生100個雙聚類，隨機挑選其中四個雙聚類來繪制它們的表達值，每個雙聚類選擇了50個基因進行繪制，如圖1所示，橫軸代表的是條件，縱軸代表的是基因表達值，從這些圖中看出每個雙聚類的基因表達值變化趨勢具有一定的相似性。在表1中分別列出了這四個雙聚類的MSR值和體積SIZE。圖2顯示了該算法的收斂曲線，橫軸代表的是迭代次數(shù)，縱軸代表的是適應度函數(shù)值，可以發(fā)現(xiàn)該算法前期全局搜索收斂較快，后期局部搜索較慢，能夠找到更優(yōu)的結(jié)果。

圖1 4個雙聚類的基因表達值

圖2 算法收斂曲線

表1 4個雙聚類的M SR、SIZ E和C I

4.4 蝙蝠算法改進前后對比

表2是基于原始二進制蝙蝠算法的雙聚類算法（bi-clustering based on original binary bat algo?rithm，BOBBA）和基于改進蝙蝠算法的雙聚類算法（BIBA）比較結(jié)果。

表2 BOBBA與BIBA的比較

其中，AMSR是所有雙聚類的平均均方殘差，MMSR是所有雙聚類中最小MSR值，ASIZE是所有雙聚類的平均體積，ACON是所有雙聚類的條件數(shù)平均值，AGENE是所有雙聚類基因平均值。并采用CI［15］對雙聚類算法結(jié)果進行評估：

CI越小，說明雙聚類的體積越大，均方殘差越小，則質(zhì)量越好。

由表2可以看出，BIBA算法無論是在A M S R方面還是在A S IZ E方面，結(jié)果都是要優(yōu)于BOBBA的，CI指標也是低于BOBBA的。

4.5 不同算法的對比分析

表3列出了七種算法在酵母菌數(shù)據(jù)集下的運行結(jié)果。評價標準和4.4的評價標準一致。

表3 7種算法結(jié)果比較

由表3可以看出，Single-object GA算法得到的AM S R最小，但是ASI ZE也是最小，因此得到的雙聚類綜合質(zhì)量并不是最好，其它算法同樣存在著相同的問題。而BIBA算法得到的雙聚類AMSR較小的同時，ASIZE也較大，C I為0.21，與其他算法相比值更低，因此得到的雙聚類的綜合質(zhì)量也更優(yōu)。

5 結(jié)語

基因表達數(shù)據(jù)隱藏著豐富的生物信息，通過雙聚類對基因表達數(shù)據(jù)矩陣進行挖掘可以發(fā)現(xiàn)這些生物信息。本文提出了一種新的基于改進蝙蝠算法的雙聚類算法。它將較新的、性能較好的智能優(yōu)化算法——蝙蝠算法應用到尋找雙聚類里。通過在酵母菌數(shù)據(jù)集上進行實驗，并與原始的蝙蝠算法和其他6種雙聚類算法的實驗結(jié)果進行比較，該方法表現(xiàn)出了較為理想的效果。但是，該算法也存在許多不足之處：計算消耗時間較長，尋優(yōu)仍需繼續(xù)改進，評價函數(shù)也需要進一步調(diào)整。對于未來的工作，將繼續(xù)修改該算法，提升算法的性能，尋找質(zhì)量更高的雙聚類。