時(shí) 磊 劉 銳 虢 韜 楊 恒 王 偉 陳 玥 張 磊 曹小鈺 羅 飛

1(貴州電網(wǎng)有限責(zé)任公司輸電運(yùn)行檢修分公司 貴州 貴陽 550005)2(國網(wǎng)電力科學(xué)研究院武漢南瑞有限責(zé)任公司 湖北 武漢 430074)3(電網(wǎng)雷擊風(fēng)險(xiǎn)預(yù)防湖北省重點(diǎn)實(shí)驗(yàn)室 湖北 武漢 430074)4(武漢大學(xué)計(jì)算機(jī)學(xué)院 湖北 武漢 430079)

0 引 言

在圖理論中，兩組對(duì)象之間的連接關(guān)系可以建模為二分圖[1-2]。其中節(jié)點(diǎn)表示對(duì)象，邊表示兩組對(duì)象節(jié)點(diǎn)之間的關(guān)系。頻繁連接的對(duì)象集合，往往是應(yīng)用背景下有意義的社團(tuán)群體。邊的連接程度是找到這些有意義社團(tuán)群體的重要特征。二分圖中，全連接二分子圖稱為完全二分子圖(biclique)。稠密的、近乎完全連接的二分子圖稱為準(zhǔn)完全二分子圖(quasi-biclique)。因?yàn)樵跀?shù)據(jù)分析中常常存在缺失、假陽性等問題[3]，相對(duì)于完全二分子圖，準(zhǔn)完全二分子圖不要求完全連接，使其具備容忍誤差的優(yōu)勢(shì)，所以它更適合實(shí)際建模。

quasi-biclique可定義靈活的結(jié)構(gòu)特征來滿足不同應(yīng)用。例如，Yan[4]提出α-quasi-biclique，要求一組中每個(gè)節(jié)點(diǎn)至少連接另一組中α%數(shù)量的節(jié)點(diǎn)。在α-quasi-biclique的基礎(chǔ)上，進(jìn)一步產(chǎn)生了δ-quasi-biclique[5]，要求一組中的每個(gè)節(jié)點(diǎn)與另一組至少連接(1-δ)的節(jié)點(diǎn)。除了限制節(jié)點(diǎn)間連接數(shù)量，還有一類應(yīng)用考慮準(zhǔn)完全二分子圖邊的權(quán)重。例如Maulik[6]提出在加權(quán)二分圖中找最大邊值的準(zhǔn)完全二分子圖。

準(zhǔn)完全二分子圖的結(jié)構(gòu)特征定義了在二分圖中搜索的目標(biāo)，搜索算法用于完成發(fā)現(xiàn)目標(biāo)準(zhǔn)完全二分子圖。大多數(shù)準(zhǔn)完全二分子圖搜索是NP問題[5]，但是相應(yīng)工作均給出了近似優(yōu)化解。Bu等[7]使用相鄰矩陣來存儲(chǔ)作用網(wǎng)絡(luò)，并用鄰接矩陣的譜去尋找特定的網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)，其中具有負(fù)特征值的特征向量映射成準(zhǔn)完全二分子圖。Maulik等[6]提出了一種多目標(biāo)遺傳算法，同時(shí)優(yōu)化了三個(gè)目標(biāo)函數(shù)來獲得目標(biāo)準(zhǔn)完全二分子圖。Liu等[8]提出了一種基于Giraph用于發(fā)現(xiàn)最大頂點(diǎn)和最大平衡準(zhǔn)完全二分子圖的挖掘算法。文獻(xiàn)[9]提出了一種具有靈活模塊的BicNET算法，以發(fā)現(xiàn)具有參數(shù)化均勻性的類quasi-biclique模塊。

現(xiàn)有的準(zhǔn)完全二分子圖結(jié)構(gòu)定義和搜索方法依然存在不足。例如，α-quasi-biclique搜索時(shí)依然要求強(qiáng)條件：以完全二分子圖作為核來擴(kuò)展。δ-quasi-biclique雖然不需要以biclique作為核，但會(huì)低估所找到的quasi-biclique。另一方面，目前quasi-biclique 模型應(yīng)用廣泛，例如生物醫(yī)學(xué)[10]、商業(yè)市場(chǎng)行為分析[11]和社會(huì)學(xué)研究[12]等。但是現(xiàn)有大多數(shù)針對(duì)解決某一具體領(lǐng)域問題提出挖掘目標(biāo)quasi-biclique的算法，缺乏推廣性，很難應(yīng)用到其他情況。

本文從三個(gè)方面定義一種通用化的，基于全局-局部密度的準(zhǔn)完全二分子圖，在滿足預(yù)設(shè)的密度閥值前提下，要求目標(biāo)quasi-biclique規(guī)模盡可能大、內(nèi)部盡可能稠密和外部盡可能稀疏，同時(shí)給出相應(yīng)的啟發(fā)式算法完成目標(biāo)準(zhǔn)完全二分子圖的搜索。

1 方 法

1.1 基于全局-局部密度的準(zhǔn)完全二分子圖

給定兩個(gè)二分圖G=和B=，且B?G。N1、N2和B1、B2表示對(duì)象節(jié)點(diǎn)集。D和DB表示兩組點(diǎn)集之間的邊集。G和B采用0-1鄰接矩陣表示。N1和B1置于行上，N2和B2置于列上。dij是第i行、第j列位置的邊連接變量，1代表有邊，0代表無邊。

quasi-biclique與biclique最大的區(qū)別在于對(duì)邊完全連接的誤差容忍。因此，與其他quasi-biclique結(jié)構(gòu)特征目標(biāo)不同，本文關(guān)注quasi-biclique邊連接的稠密程度特征，通過定義quasi-biclique的內(nèi)部密度θ度量該特征。

(1)

當(dāng)B是G中基于全局-局部密度的quasi-biclique時(shí)，其密度θ必須大于給定閥值,為了使quasi-biclique密度特征更顯著，綜合文獻(xiàn)[13-14]，還要求其滿足三個(gè)目標(biāo)：

1) 最大規(guī)模：不存在滿足θ閾值的另一個(gè)quasi-bicliqueB′，使得B?B′；

2) 內(nèi)部稠密：對(duì)于有相同規(guī)模的quasi-biclique，其密度盡可能大；

3) 外部稀疏：對(duì)于具有相同規(guī)模的quasi-biclique，與G-B圖中節(jié)點(diǎn)的聯(lián)系邊盡量少。

為了在給定的G中搜索到基于全局-局部密度的quasi-biclique，定義以下與密度相關(guān)的參量引導(dǎo)啟發(fā)式算法搜索目標(biāo)。

σr和σc度量B中行r的行密度和列c的列密度：

(2)

(3)

類似的，τr和τc度量N1-B1中行r的行密度和N2-B2中列c的列密度：

(4)

(5)

φr和φc度量N1-B1中行r的行內(nèi)部一致性和N2-B2中列c的列內(nèi)部一致性：

(6)

(7)

(8)

T1={r|r∈N1-B1,τr>0}T2={c|c∈N2-B2,τc>0}

1.2 挖掘方法

由于在二分圖中搜索quasi-biclique已經(jīng)被證明是NP問題，因此本文提出一種啟發(fā)式算法來搜索基于全局-局部密度的目標(biāo)quasi-biclique。算法流程如下所示。

Input: a matrix X of size N1*N2, a threshold ρ for mini-mal inner density, a threshold s for minimal row/colume in-cluded in quasi-biclique.Output: quasi-bicliquesFlag=Truewhile Flag doFlag=False, B=DeleteZeroDensity(X)D[ ]=DivideMatrix(B)for the ith matrix from 1 to length of D M1[i]=DetectKernel (ρ, s, D[i])endforif exist a matrix m1 in M1 whose density θ >ρ then find the matrix m1 whose |r||c| is maximal in M1 Q=ExpandKernel(ρ, s, m1,X) if Q is not null then output Q X=X-Q flag=True endifendifendwhileFunction DetectKernel (ρ, s, M){Calculate θ of M.while θ<ρ and |R| >s and |C|> s do Find the rows in M with the lowest σr Find the columns in M with the lowest σc if σr <σc then Remove the rows with σr from M endif if σr >σc then Remove the columns with σc from M endif if σr==σc and |R|≥|C| then Remove the rows with σr from M endif if σr ==σc and |R|<|C| then Remove the columns with σc from M endif Calculate the inner density θ of shrink matrix Mendwhilereturn M}

Function ExpandKernel (ρ, s, M, X){K=Mwhileθ>ρ and |R|>s and |C|>s do Keep the old record M=K Find rows outside K but in X with largest φr Find columns outside K but in X with largest φcif φr >φc then the rows with largest φr on K is added into Kendifif φr <φc then add the columns with largest φc into Kendifif φr==φc and |R|≥|C| then the columns with largest φr on K is added into Kendifif φr ==φc and |R|<|C| then the rows with largest φc on K is added into K endif Calculate the inner density θ of matrix K endwhile return M}

以上啟發(fā)式搜索算法通過內(nèi)部密度，行和列密度，以及內(nèi)核引用的行和列內(nèi)部一致性決定哪些行或列被刪除(在內(nèi)核搜索中)或添加(在內(nèi)核擴(kuò)展中)進(jìn)行搜索目標(biāo)。整個(gè)過程包括四個(gè)主要部分。

1) 去空行/列 DeleteZeroDensity去除密度為零的行和列。如果輸入的矩陣是稀疏矩陣，DeleteZeroDensity可以降低數(shù)據(jù)處理的規(guī)模。

2) 對(duì)角化判斷與分割 DivideMatrix檢查待處理的二分圖鄰接矩陣是否存在對(duì)角化情形。如果存在，則分割,返回非對(duì)角化的所有子矩陣構(gòu)成的數(shù)組。鄰接矩陣存在對(duì)角化會(huì)影響內(nèi)核發(fā)現(xiàn)和擴(kuò)展。因此，要求內(nèi)核的搜索是在非對(duì)角化的矩陣中進(jìn)行。

圖1展示了二分圖G的鄰接矩陣是對(duì)角化(反對(duì)角化類似)的例子。G是對(duì)角化，當(dāng)且僅當(dāng)存在G=A∪B。其中A=，B=，其中A1∩B1=Φ，A2∩B2=Φ，D=DA+DB，G′=中的元素均為零。

圖1 對(duì)角化舉例

假設(shè)G的鄰接矩陣D是如圖1所示可以被分成兩個(gè)小的非對(duì)角矩陣A和B。在搜索quasi-biclique的內(nèi)核時(shí)，需要使內(nèi)核的內(nèi)部密度不斷提升，直至滿足預(yù)置的密度閥值。但是如果在對(duì)角化的不同子矩陣中尋找內(nèi)核，這樣會(huì)使密度減小，所以有如下定理1。

定理1將B中的任何行或列添加到A中，則新矩陣密度必有r2

證明：以添加一行為例。假設(shè)k1是矩陣A中的值為1的個(gè)數(shù)，那么按式(1)的計(jì)算方法，A的密度r1有：

設(shè)k2是新加入行中的1的個(gè)數(shù)，則加入行后的新矩陣密度r2為：

因?yàn)樗忻芏葹榱愕男泻土卸家褟腄中刪除，則有:

k1≥|A1|,k1≥|A2|

為了證明r1恒大于r2，可以轉(zhuǎn)換為證明r1-r2總大于零:

由于分母恒大于零，所以只需要考慮分子。展開分子后得到:

Numerator=k1|A1||A2|+k1|A1|k2+k1|A2|+

k1k2-|A1||A2|k1-|A1||A2|k2=

k1s2+k1k2+k2|A1|(k1-|A2|)

因?yàn)閗1≥|A1|，即分子總是大于零，所以r2

以上證明說明了需要對(duì)原始輸入矩陣進(jìn)行對(duì)角化檢查與分割。

3) 搜索quasi-biclique內(nèi)核 DetectKernel負(fù)責(zé)確定quasi-biclique的內(nèi)核。DetectKernel是在經(jīng)過分割后的非對(duì)角化子矩陣上通過移除低密度部分找到滿足密度要求的quasi-biclique無誤差子結(jié)構(gòu)。相比其他方法,本方法并不是搜索無誤差的biclique結(jié)構(gòu)，增加了方法適用的靈活性。具體搜索quasi-biclique內(nèi)核時(shí)，是以行、列密度作為指標(biāo)，移除低密度行或者列。當(dāng)行列密度相等時(shí)，算法試圖優(yōu)先去除行列維數(shù)較大的一維。這樣可以使得搜尋的quasi-biclique內(nèi)核趨向于方陣。在有些具體的應(yīng)用中，當(dāng)期望行、列位置上代表的對(duì)象有均等的機(jī)會(huì)出現(xiàn)在最終quasi-biclique中時(shí)，適合采用此策略。當(dāng)兩組對(duì)象之間有重要度差異，相應(yīng)地可以考慮優(yōu)先保留重要一組對(duì)象。每當(dāng)一輪去除低密度部分后，將計(jì)算當(dāng)前狀態(tài)下的內(nèi)部密度值，當(dāng)密度要求滿足閥值時(shí)，搜索內(nèi)核停止。

4) 擴(kuò)展內(nèi)核 以DetectKernel確定的內(nèi)核為基礎(chǔ)，ExpandKernel在原始輸入的二分圖鄰接矩陣上確定最終目標(biāo)quasi-biclique。由于DetectKernel確定內(nèi)核過程中可能去除了一些有用的行或者列，所以內(nèi)核擴(kuò)展是增加在內(nèi)核搜索時(shí)誤刪的行或列。在擴(kuò)展過程中，需要考慮內(nèi)核外部行列與內(nèi)核的一致程度，從而決定哪些外部的行或者列可以加入到內(nèi)核中。所以，此步驟操作中是以內(nèi)核外行或者內(nèi)核外列的內(nèi)部一致性指標(biāo)作為依據(jù)，選擇擴(kuò)展部分。擴(kuò)展過程中，伴隨著行、列的增加，quasi-biclique密度會(huì)隨之下降，當(dāng)不再滿足閥值要求時(shí)，對(duì)該內(nèi)核的擴(kuò)展搜索停止。

2 實(shí)驗(yàn)與討論

本文使用R語言實(shí)現(xiàn)搜索基于全局-局部密度quasi-biclique的算法。算法在個(gè)人PC上實(shí)現(xiàn)，硬件平臺(tái)CPU為Intel i5-4440，RAM為16 G。本算法的主要操作是比較矩陣的行、列之間的相似性與一致性，所以在R的程序?qū)崿F(xiàn)中，通過apply函數(shù)代替所有循環(huán)操作，從而極大地提高了程序的求解速度。

2.1 對(duì)比實(shí)驗(yàn)

如前所述，大部分quasi-biclique挖掘方法是整合在具體應(yīng)用中的，沒有獨(dú)立的源程序，這樣不能準(zhǔn)確地復(fù)現(xiàn)或者直接應(yīng)用到其他領(lǐng)域。為了客觀比較性能，選擇提供源碼的PPIExtend[6]與本方法進(jìn)行比較。比較實(shí)驗(yàn)使用的數(shù)據(jù)集是酵母蛋白質(zhì)相互作用關(guān)系。它包括了4 959個(gè)蛋白質(zhì)和它們之間10 640條作用邊。當(dāng)輸入?yún)?shù)為(0.1, 5, 5)時(shí)，PPIExtend耗時(shí)約500秒搜索到了59 124個(gè)quasi-biclique。在(0.9, 2, 2)參數(shù)下，本方法耗時(shí)約70秒搜索到8個(gè)quasi-biclique。為區(qū)別兩個(gè)方法的結(jié)果，PPIExtend搜索到的quasi-biclique記為L(zhǎng)QB，本方法記為GLDQB。GLDQBs和LQBs分別含有109和155個(gè)蛋白質(zhì)，其中71個(gè)蛋白質(zhì)相同。

首先，比較分析PPIExtend方法和本文方法挖掘結(jié)果相互覆蓋度。對(duì)于兩個(gè)quasi-biclique：X1=和X2=，如果A2∪B2中90%以上的蛋白質(zhì)在A1∪B1中被發(fā)現(xiàn)，則說X1覆蓋X2，反之亦然。對(duì)比兩個(gè)方法的結(jié)果，只有一個(gè)GLDQB(占整體的12.5%)被PPIxtend覆蓋，而對(duì)于PPIExtend，有21 286個(gè)LQB(占整體的36.0%)被本方法覆蓋。這一事實(shí)說明本方法可以使用更少的GLDQB來表示挖掘到目標(biāo)信息，比LQB更有效。PPIExtend的結(jié)果過于瑣碎，如此龐大的結(jié)果對(duì)于使用者直接理解其含義是困難的。

然后，考察本方法的條件約束效果。表1顯示GLDQB關(guān)鍵指標(biāo)項(xiàng)，依次為：內(nèi)部密度、外部密度、背景(以LQB為參照，具有與對(duì)應(yīng)GLDQB相同大小的LQB個(gè)數(shù))和外部個(gè)數(shù)(背景中外部密度小于GLDQB的LQB個(gè)數(shù))。

表1 GLDQBs主要指標(biāo)

表1中所有的GLDQB的內(nèi)部密度要遠(yuǎn)大于其外部密度。1、2、6、7和8號(hào) GLDQB沒有相應(yīng)的背景LQB,說明本方法發(fā)現(xiàn)了與PPIExtend不同的quasi-biclique?！巴獠總€(gè)數(shù)”的數(shù)目反映出本方法能恰好滿足密度閥值的情況下使得quasi-biclique規(guī)模最大，而且“外部個(gè)數(shù)”占“背景”比重小，反映LQB的外部密度稀疏性不如GLDQB，從而體現(xiàn)出GLDQB出的全局-局部密度的內(nèi)外部密度的顯著性特點(diǎn)。

2.2 藥物和目標(biāo)基因數(shù)據(jù)應(yīng)用

藥物可以對(duì)多個(gè)基因起作用，一個(gè)目標(biāo)基因可能受多種藥物的影響[15-16]。本文不考慮多肽類藥物，從DrugBank 數(shù)據(jù)庫中收集已批準(zhǔn)的針對(duì)人類疾病的1 186種小分子藥物與1 141個(gè)目標(biāo)基因，總計(jì)4 594種藥物-基因相互作用。通過利用準(zhǔn)完全二分子圖可以幫助了解藥物和目標(biāo)基因之間的頻繁關(guān)系，同時(shí)可推測(cè)未知的藥物與目標(biāo)基因作用關(guān)系。

準(zhǔn)完全二分子圖搜索算法中需要設(shè)定的主要參數(shù)是ρ，它指定最終quasi-biclique必須滿足的最小密度要求。該參數(shù)的大小會(huì)影響搜索時(shí)間。為了研究不同的ρ與運(yùn)行時(shí)間的關(guān)系，ρ按步長(zhǎng)0.1，從0.1到0.9依次執(zhí)行搜索。圖2是不同ρ下的運(yùn)行時(shí)間對(duì)比。通過編程優(yōu)化，quasi-bicluque搜索算法的實(shí)際運(yùn)行時(shí)間是可接受的。當(dāng)ρ為0.1時(shí)，用時(shí)最長(zhǎng)，需要37秒。同時(shí),結(jié)果表明算法的運(yùn)行時(shí)間并不是與ρ成完全反比。

圖2 不同ρ取值下的運(yùn)行時(shí)間

圖3是上述實(shí)驗(yàn)在對(duì)應(yīng)ρ下挖掘準(zhǔn)完全二分子圖個(gè)數(shù)。當(dāng)ρ增加時(shí)，密度閥值增加，所以quasi-biclique的數(shù)量首先降低，相應(yīng)地運(yùn)行時(shí)間也會(huì)降低。當(dāng)ρ超過0.7時(shí),矩陣中的中等密度區(qū)域被分解成許多小密度的quasi-biclique。因此，quasi-biclique的數(shù)量又會(huì)增加，需要的時(shí)間也會(huì)相應(yīng)增加。

圖3 對(duì)應(yīng)ρ取值下的準(zhǔn)完全二分子圖數(shù)

圖4是以對(duì)應(yīng)ρ值發(fā)現(xiàn)的每個(gè)quasi-biclique包含的藥物和目標(biāo)基因數(shù)。當(dāng)ρ較小時(shí)，發(fā)現(xiàn)了代表藥物數(shù)與目標(biāo)基因數(shù)具有不平衡現(xiàn)象的quasi-biclique，即當(dāng)ρ較小時(shí)，可發(fā)現(xiàn)一個(gè)藥物針對(duì)多個(gè)目標(biāo)基因，或者是一個(gè)基因可被多個(gè)藥物作用的情形。隨著ρ值的增加，發(fā)現(xiàn)的是代表藥物與目標(biāo)基因數(shù)趨向平衡的quasi-biclique，這與較小ρ的情形不同，表明ρ值與quasi-biclique的結(jié)點(diǎn)平衡特征也有關(guān)聯(lián)性。

圖4 不同ρ取值下,每個(gè)quasi-biclique包含的藥物數(shù) 與目標(biāo)基因數(shù)

生物意義上，取ρ=0.3的quasi-biclique分析生物功能特點(diǎn)。此參數(shù)下的quasi-biclique的基因包括許多典型的基因家族，如ATP結(jié)合盒(ABC)轉(zhuǎn)運(yùn)蛋白、?；o酶A脫氫酶家族、鋅金屬酶、鈣電壓門控通道亞基、G蛋白偶聯(lián)受體等。quasi-biclique起到橋梁的作用，將某些目標(biāo)基因家族與相應(yīng)的藥物治療作用聯(lián)系起來。例如，第16個(gè)quasi-biclique中包含了藥物DB00143、DB01224、DB00404、DB00364、DB00211、DB01068、DB00245、DB05630、DB00521、DB00892、DB00517、DB01165、DB01104，藥物功效包括了硫酸化蔗糖、治愈消化性潰瘍、營養(yǎng)補(bǔ)充劑、抗高血壓藥、抗焦慮藥、抗生素、消化藥物、抗精神病藥、抗驚厥藥和抗焦慮藥。針對(duì)基因有DRD4、KCNH2、RARA、RARB、RARG、RXRB、RXRG、HTR1B、HTR1E、HTR3A、HRH1、ADRA1A、ADRA1B。其相應(yīng)的目標(biāo)基因具有典型的精神病相關(guān)基因，如DRD4、HTR1B、ADRA1A。該quasi-biclique的內(nèi)部密度為0.302，外部密度是0.031，內(nèi)部密度相對(duì)外部的具有顯著密集特性。該quasi-biclique中還有很多藥物與目標(biāo)基因之間沒有報(bào)道有相關(guān)聯(lián)系(即該連接元素為0)，但是根據(jù)quasi-biclique中基因簇的功能相似，以及quasi-biclique內(nèi)部0.3的密集聯(lián)系的事實(shí)，該quasi-biclique內(nèi)這些還未報(bào)道的藥物與基因?qū)χg客觀上有很高的可能性是未知的藥物與目標(biāo)基因之前的作用關(guān)系，這可以為藥物與基因研究提供提示。

3 結(jié) 語

本文考慮到局部特征(內(nèi)部相關(guān))和全局特征(外部稀疏)，提出了基于局部密度(內(nèi)部密度)和全局密度(外部密度)的quasi-biclique。它可以在大型二分圖中捕獲重要的子結(jié)構(gòu)。通過分析和實(shí)驗(yàn)證明了基于全局-局部密度的準(zhǔn)完全二分子圖是發(fā)現(xiàn)雙向關(guān)系合適的模型，其搜索算法也是有效的。介于基于全局-局部密度的準(zhǔn)完全二分子圖是一個(gè)通用的模型，它可以推廣到社交網(wǎng)絡(luò)、電力系統(tǒng)等許多其他領(lǐng)域。