基于全局-局部密度的準(zhǔn)完全二分子圖建模與挖掘方法

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

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

0 引 言

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

quasi-biclique可定義靈活的結(jié)構(gòu)特征來滿足不同應(yīng)用。例如，Yan[4]提出α-quasi-biclique，要求一組中每個節(jié)點至少連接另一組中α%數(shù)量的節(jié)點。在α-quasi-biclique的基礎(chǔ)上，進一步產(chǎn)生了δ-quasi-biclique[5]，要求一組中的每個節(jié)點與另一組至少連接(1-δ)的節(jié)點。除了限制節(jié)點間連接數(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]使用相鄰矩陣來存儲作用網(wǎng)絡(luò)，并用鄰接矩陣的譜去尋找特定的網(wǎng)絡(luò)拓撲結(jié)構(gòu)，其中具有負特征值的特征向量映射成準(zhǔn)完全二分子圖。Maulik等[6]提出了一種多目標(biāo)遺傳算法，同時優(yōu)化了三個目標(biāo)函數(shù)來獲得目標(biāo)準(zhǔn)完全二分子圖。Liu等[8]提出了一種基于Giraph用于發(fā)現(xiàn)最大頂點和最大平衡準(zhǔn)完全二分子圖的挖掘算法。文獻[9]提出了一種具有靈活模塊的BicNET算法，以發(fā)現(xiàn)具有參數(shù)化均勻性的類quasi-biclique模塊。

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

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

1 方 法

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

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

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

(1)

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

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

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

3) 外部稀疏：對于具有相同規(guī)模的quasi-biclique，與G-B圖中節(jié)點的聯(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)核擴展中)進行搜索目標(biāo)。整個過程包括四個主要部分。

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

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

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

圖1 對角化舉例

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

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

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

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

因為所有密度為零的行和列都已從D中刪除，則有:

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|)

因為k1≥|A1|，即分子總是大于零，所以r2

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

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

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

2 實驗與討論

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

2.1 對比實驗

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

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

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

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

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

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

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

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

圖2 不同ρ取值下的運行時間

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

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

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

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

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

3 結(jié) 語

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