李力 曹以誠 毛曉帆

(華南理工大學 生物科學與工程學院，廣東 廣州 510006)

對由節(jié)點和邊組成的圖進行劃分得到一系列的子圖，子圖內(nèi)的連接比較緊密而子圖間的連接相對稀疏，這一劃分過程稱為圖聚類，所得子圖稱為聚類［1］.2008年，Kawamura 等［2］首次將圖聚類方法應(yīng)用在基因芯片數(shù)據(jù)分類分析中，將DNA 芯片數(shù)據(jù)使用圖聚類處理后構(gòu)建分類器，分類精度較二分類法提高80%，并成功在由38 個對照組樣組成的分類簇中發(fā)現(xiàn)了胸腺癌相關(guān)基因，顯示了圖聚類方法在高維度數(shù)據(jù)學習和生物學對象相互關(guān)系分類方面的優(yōu)勢.基因相互關(guān)系網(wǎng)絡(luò)的數(shù)學模型是一個無向圖［3］，圖中的節(jié)點表示基因，邊線表示相應(yīng)基因表達譜的相似度.因此，基因表達數(shù)據(jù)的聚類問題可以簡化為圖論問題，這也是圖聚類方法可以較好應(yīng)用于基因芯片數(shù)據(jù)分析的基礎(chǔ).同時，基因網(wǎng)絡(luò)還是一個特殊的復雜網(wǎng)絡(luò)，節(jié)點數(shù)目龐大且節(jié)點之間的相互關(guān)系錯綜復雜［4］.在這個網(wǎng)絡(luò)中，模塊是實現(xiàn)各種生物功能的主體，各個模塊都承擔著不同的生物學功能，所以越來越多的研究者使用模塊識別的方法對基因芯片數(shù)據(jù)進行圖聚類分析.但是，一直以來這些模塊識別方法缺少一種定量的評價方法.2004年，Newman 等［5］基于協(xié)調(diào)混合(Associative Mixing)理論提出一個衡量網(wǎng)絡(luò)模塊劃分質(zhì)量的指標——模塊性(Modularity)，在這個指標下，較大的模塊性值代表著對網(wǎng)絡(luò)有一個較好的劃分.此外，由于模塊識別對單節(jié)點的分配過程會導致算法陷入局部最優(yōu)解，Mei等［6］在模塊性圖聚類過程中加入了一個擾亂概率因子.加入擾亂概率因子的優(yōu)點在于可以隨機破壞前一次聚類過程中輸出的子圖結(jié)構(gòu)，以搜索全局最優(yōu)解;但缺點也很明顯，即在多次迭代運算之后，對已輸出子圖結(jié)構(gòu)會產(chǎn)生無目的的破壞，大大影響了結(jié)果的穩(wěn)定性.

針對現(xiàn)有模塊識別的定量評價方法存在的不足，文中提出了一種基于模塊性指標和子圖平滑度(Smoothness)的全局圖聚類方法(Module Smoothness)，為了防止算法陷入局部最優(yōu)解，尋找最優(yōu)的全局模塊性，文中引入子圖平滑度的定義，對每一次聚類結(jié)果產(chǎn)生的異常值比平均水平多的子圖進行打散，再利用打散得到的單節(jié)點進行下一次聚類，多次迭代后得到具有較強圖結(jié)構(gòu)的基因網(wǎng)絡(luò)的劃分，即最優(yōu)聚類結(jié)果.文中還利用Module Smoothness 方法對一組人外周血淋巴細胞基因組表達數(shù)據(jù)進行了分析，并將聚類結(jié)果與其他常用基因芯片數(shù)據(jù)聚類方法進行了比較.

1 材料與方法

1.1 芯片數(shù)據(jù)與預(yù)處理

從GEO 數(shù)據(jù)庫［7］中選用放療毒性與DNA 轉(zhuǎn)錄異常研究實驗數(shù)據(jù)集GSE1725，共包含171 個人外周血淋巴細胞基因組表達數(shù)據(jù)樣本，來自57 個病人，芯片平臺Affymetrix Human Genome U95 V2，每個樣本包含6 808 個基因表達數(shù)據(jù)［8］.所有基因原始表達數(shù)據(jù)轉(zhuǎn)換為對數(shù)，使用Lowess 方法［9］進行歸一化處理.

1.2 圖聚類

圖聚類過程通常分為兩步:圖構(gòu)造和圖劃分.如圖1 所示，首先根據(jù)基因表達數(shù)據(jù)集計算每個表達譜之間的相關(guān)系數(shù)，即圖結(jié)構(gòu)的邊，得到一個抽象意義上的總的基因相互關(guān)系圖，圖中的節(jié)點和邊構(gòu)成的矩陣稱為模式特征矢量集;再根據(jù)合適的距離聚類算法找到基于聯(lián)接的大型圖中聯(lián)接良好的類似小集團的最佳簇，也即具有生物學意義的相互關(guān)聯(lián)的基因聚類.

圖1 圖聚類原理示意圖Fig.1 Schematic diagram of graph clustering

1.2.1 相關(guān)系數(shù)計算

在基因表達數(shù)據(jù)的聚類分析中，通常使用皮爾森公式［10］計算基因間的相關(guān)系數(shù).但是，該公式對正向反向共表達數(shù)據(jù)的度量是建立在統(tǒng)計基礎(chǔ)上的，無法反映基因表達數(shù)據(jù)中實際具有的時序性特征.為此，文中采用了王晗［11］提出的基于打分的相似度計算方法，具體算法簡要描述如下:

首先，將原始數(shù)據(jù)矩陣轉(zhuǎn)換為差分矩陣，對樣本總數(shù)N，整個表達譜的X 軸由樣本序列s(s=1，2，…，N)變成統(tǒng)一的單位時間序列t(t=1，2，…，S -1)，由此連接樣本間不同表達量的線段與X 軸形成S-1 個夾角θt，將θt的二倍角正弦值代入式(1)中得到任意兩個基因(i，j)在相鄰兩樣本間表達變化量的相似度Gij(t).

隨后，進行基因表達變化趨勢的相似度計算.如式(2)所示，函數(shù)fij(t)對基因i 和基因j 在所有時間段內(nèi)的變化趨勢進行打分，如果兩個單位時間內(nèi)兩基因的變化趨勢一致，則該時段的打分為前一段打分的值加上常數(shù)c，如不一致則重置為初始值.

最后，將前面計算得到的表達值相似度與變化趨勢相似度使用式(3)結(jié)合起來，將每個時段內(nèi)基因i和基因j 的變化量差異值和對變化趨勢的打分相乘，然后進行累加最后得到統(tǒng)一的相關(guān)系數(shù)rij(t).

式中，a、b 均為常數(shù).

式中，rij代表基因i 和基因j 之間統(tǒng)一的相關(guān)系數(shù).

1.2.2 模塊性的定義

通過上述方法將基因與基因之間相關(guān)系數(shù)映射得到加權(quán)構(gòu)造圖G=(V，E)，其中節(jié)點v(V)表示基因，e(E)表示基因表達譜間的距離.圖聚類的目標是將構(gòu)建的網(wǎng)絡(luò)圖劃分為若干相互關(guān)系更為緊密的子圖，最終使得子圖內(nèi)部節(jié)點相似度較高，而不同子圖的節(jié)點之間相似度較低.模塊性就是為了量化這種網(wǎng)絡(luò)結(jié)構(gòu)而提出的一個衡量網(wǎng)絡(luò)模塊劃分的標準［5］.對于一個給定了劃分、由n 個節(jié)點和m 條邊組成的圖，其模塊性(Q 值)計算公式如下:

式中，k 代表子圖的個數(shù)，eI是子圖I 內(nèi)部所有邊的數(shù)目，aI是子圖I 中所有節(jié)點的度之和.

1.2.3 基于平滑度的最優(yōu)模塊性值計算

模塊性圖聚類的通常做法是:先將每一個節(jié)點作為單獨的子圖，然后在這個初始劃分Cmax的基礎(chǔ)上，依次將每個節(jié)點移動到使Q 值增長最快的其他子圖中去，直至整個網(wǎng)絡(luò)的Q 值不能再增大為止［6］.至此，所有的節(jié)點均處于較合適的子圖中，但這只是相對于單節(jié)點移動的一個局部最優(yōu)解.為尋求最優(yōu)的全局模塊性，對每次聚類得到的所有子圖進行平滑度的計算.理論上關(guān)系相對集中的基因節(jié)點蔟，各節(jié)點之間的表達譜相關(guān)系數(shù)(即邊線距離)應(yīng)該服從正態(tài)分布的原則，其方差值越小，各節(jié)點間的關(guān)聯(lián)程度就越集中，該子圖的平滑度就越高，也意味著簇內(nèi)基因的表達情況有更為相似的規(guī)律;反之，蔟內(nèi)基因關(guān)系的異常值越多，方差值越大，子圖的平滑度就越低，簇內(nèi)基因的關(guān)聯(lián)性較差.子圖平滑度的計算方法如下.

式中，S 代表子圖平滑度，m 代表所有邊數(shù)，RI表示第i 個邊的邊長，ˉR 為所有邊長的算術(shù)平均值.

由此得到每個子圖的平滑度S 和所有子圖平滑度的平均值ˉS，將平滑度低于總體平均值的子圖打散成單個節(jié)點，再對所有單節(jié)點子圖重復聚類過程.整個算法流程可以看作是收縮和打散步驟迭代多次求全局最優(yōu)解的過程，其形式化定義如下:

Module Smoothness 算法:

1.將每一個節(jié)點作為一個單獨的子圖，記錄當前的子圖結(jié)構(gòu):Cbest和當前結(jié)構(gòu)模塊性值:Qbest

2.repeat /* 重復收縮和打散步驟* /

3.repeat /* 重復聚類過程* /

4.for each clustering x /* 遍歷所有子圖* /

6.end for

7.until Q 不再增加

8.if Q＞Qbest，更新Cbest和Qbest

9.for each clustering x /* 遍歷所有子圖* /

10.if S ＜ˉS，將子圖x 中每一個節(jié)點作為一個單獨的子圖

11.end for

12.until Qbest不再增加在上述算法中，表示將節(jié)點y 合并后子圖內(nèi)的所有邊長和增加值最大的子圖α.

1.3 對比分析

為了考察Module Smoothness 方法的聚類效能，將該方法與4 種常用的基因表達數(shù)據(jù)聚類方法相比較，考察這些方法在對同一數(shù)據(jù)聚類分析過程中的類間重疊度和品質(zhì)因子值.4 種常用聚類方法為:k均值聚類方法、自組織(SOM)聚類方法、Fuzzy 模糊聚類方法和經(jīng)典圖聚類方法(使用皮爾森公式和Ncut 算法［12］).

品質(zhì)因子FOM 是Yeung 等［13］提出的用來比較不同方法聚類效能的工具，它用考察類內(nèi)變異總和的方法來判斷類的質(zhì)量水平，從而找出最優(yōu)的聚類算法.FOM 方法也能夠幫助觀察實際數(shù)據(jù)中的最優(yōu)分類數(shù)，是評價聚類方法質(zhì)量的有效工具.但是，隨著聚類數(shù)的增加，F(xiàn)OM 值有自然增大的趨勢，因此，文中使用易東等［14］改進的FOM 計算方法來考察5種聚類方法的聚類效能，改進的計算方法如下:

式中，F(xiàn)OM'表示改進后的品質(zhì)因子，n 表示總節(jié)點數(shù)，此處代表總基因數(shù)，k 代表當前分類的簇數(shù).

2 結(jié)果與分析

2.1 聚類性能分析

類間重疊度常用來度量不同聚類間的分離程度，對于同一組數(shù)據(jù)而言，不同聚類方法的分類結(jié)果的類間重疊度低，就說明方法的聚類表現(xiàn)相對優(yōu)于其他方法［15］.5 種算法在放療毒性實驗數(shù)據(jù)集上的類間重疊度變化如圖2(a)所示，從圖中可知:Module Smoothness 方法的重疊度曲線整體上比其他4 種聚類方法更靠近X 軸，聚類表現(xiàn)優(yōu)于這些方法;隨著聚類數(shù)量的增加，該方法的重疊度值變化更加穩(wěn)定.

圖2(b)中顯示的是5 種算法在GSE1725 數(shù)據(jù)集上的FOM'變化曲線圖，通過對算法的FOM'變化曲線進行比較，可以判斷各聚類算法的整體效能.如圖2(b)所示，Module Smoothness 方法的聚類效能總體上要優(yōu)于其他4 種算法;在聚類數(shù)到達39 類之后，大多數(shù)算法的下降趨勢均趨于平穩(wěn)，可以認為該數(shù)據(jù)集最佳聚類數(shù)在39 左右，因此對5 種方法在聚類數(shù)為39 時的FOM'值進行比較，Module Smoothness 方法的FOM' 值比經(jīng)典圖聚類方法、k 均值、SOM 及Fuzzy 算法分別低28.41%、19.21%、9.84%和24.67%;此外，經(jīng)典圖聚類算法因為只是從統(tǒng)計學意義上對表達數(shù)據(jù)進行聚類，類內(nèi)變異度明顯大于其他方法，聚類質(zhì)量在5 種方法內(nèi)最差.

圖2 不同算法聚類分析外周血淋巴細胞基因表達數(shù)據(jù)的比較Fig.2 Comparison of the Peripheral blood lymphocytes genes expression data analyzed with various clustering algorithms

2.2 聚類結(jié)果

網(wǎng)絡(luò)智能檢索系統(tǒng)GeneCards 是魏茲曼科學研究院建立的提供以基因為中心的信息自動整合及挖掘工具［16］.為了驗證Module Smoothness 算法的聚類結(jié)果，從GeneCards 網(wǎng)站得到與ATF-2 轉(zhuǎn)錄因子通路相關(guān)的77 個基因，查找這些基因的聚類結(jié)果，其中12 個基因(SMAD2，TCF3，MYC，VDR，CEBFB，CDK2，AKT1，CDK4，TGIF2，RBBP4，HDAC1，SMAD4)被正確地聚類在一起.這12 個與ATF-2 轉(zhuǎn)錄因子通路相關(guān)基因的聚類結(jié)果見圖3，從圖3 中可以看出，這些聚在同一類中的基因在樣本中表達的時序變化規(guī)律非常相近.

圖3 12 個與ATF-2 轉(zhuǎn)錄因子通路相關(guān)基因的聚類結(jié)果Fig.3 Clustering results of 12 genes related to ATF-2 transcription factor pathway

GSE1725 試驗集作者Rieger 等［8］使用層次聚類的方法篩選出52 個相關(guān)基因，并按照基因功能將它們分為8 類.文中將其作為參考標準，使用性能分析結(jié)果中聚類效能較近的Module Smoothness 和自組織算法對這52 個基因進行聚類，所得3 種不同算法的聚類結(jié)果如表1 所示.

表1 3 種不同算法聚類結(jié)果1)Table 1 Clustering results of 3 different clustering algorithm

從表1 可知:自組織聚類結(jié)果與原文中的層次聚類結(jié)果相差不大，都將同屬于DNA 修復功能類的RAD23B 和PSF 基因分到不同的兩個類，而文中所建立的方法將RAD23B 和PSF 基因準確地聚在了同一個類中;此外，同屬細胞凋亡功能的TNFSF7 和SLC25A5 基因等也得到了正確的聚類.

以上的聚類結(jié)果分析說明，Module Smoothness方法從基因芯片數(shù)據(jù)中提取到了有效的生物學信息，將具有相近生物功能和表達譜具有相似時序變化規(guī)律的基因正確地聚類到了一起，其分類準確度優(yōu)于層次聚類和自組織聚類方法.

2.3 算法執(zhí)行效率

基因芯片數(shù)據(jù)集中的基因數(shù)量和樣本數(shù)量通常都非常龐大，分析過程中需要耗費大量的計算機資源，因而對聚類算法的執(zhí)行效率提出了較高的要求.使用Module Smoothness 與前述的經(jīng)典圖聚類方法、k 均值、SOM 及Fuzzy 算法，在相同配置的計算機上對GSE1725 數(shù)據(jù)集的全部6808 個基因進行聚類分析，結(jié)果見表2.

表2 不同算法的聚類執(zhí)行效率Table 2 Execution efficiency of various clustering algorithm

從表2 可知，除經(jīng)典圖聚類外，其他4 種算法聚類數(shù)都在44 左右;其中，k 均值憑借局部搜索尋找最優(yōu)解方法完成聚類，所耗時間最少，而Module Smoothness 在保證聚類結(jié)果的條件下，執(zhí)行效率和SOM 算法相近，比Fuzzy 算法高5.94%.

3 結(jié)語

文中提出了一種基于模塊性指標和子圖平滑度的全局圖聚類方法Module Smoothness，為了防止算法陷入局部最優(yōu)解，尋找最優(yōu)的全局模塊性，引入子圖平滑度的定義，對每一次聚類結(jié)果產(chǎn)生的異常值比平均水平多的子圖進行打散，再利用打散得到的單節(jié)點進行下一次聚類，多次迭代后得到具有較強圖結(jié)構(gòu)的基因網(wǎng)絡(luò)的劃分，即最優(yōu)的聚類結(jié)果.利用該方法與經(jīng)典圖聚類、k 均值、SOM 及Fuzzy 4 種常用的基因芯片數(shù)據(jù)聚類方法，對一組人外周血淋巴細胞基因組表達數(shù)據(jù)進行了分析比較.結(jié)果表明:該方法在聚類過程中的平均類間重疊度和FOM'值總體上優(yōu)于其他4 種算法，在將數(shù)據(jù)集分類到最佳聚類數(shù)39 時，其FOM' 值分別比上述4 種方法低28.41%、19.21%、9.84% 和24.67%;聚類結(jié)果的生物信息學分析顯示，該方法將具有相近生物功能和表達譜具有相似時序變化規(guī)律的基因準確地聚類到了一起，其分類準確度高于層次聚類和自組織聚類方法;在算法的執(zhí)行效率方面，Module Smoothness在保證聚類結(jié)果的條件下，執(zhí)行效率和SOM 算法相近，比Fuzzy 算法高5.94%.

［1］鄧健爽，鄭啟倫，彭宏，等.基于連通圖動態(tài)分裂的聚類算法［J］.華南理工大學學報:自然科學版，2007，35(1):118-122.Deng Jian-shuang，Zeng Qi-lun，Peng Hong，et al.Clustering algorithm based on dynamic division of connected graph［J］.Journal of South China University of Technology:Natural Science Edition，2007，35(1):118-122.

［2］Kawamura T，Mutoh H，Tomita Y，et al.Cancer DNA microarray analysis considering multi-subclass with graphbased clustering method ［J］.Journal of Bioscience and Bioengineering，2008，106(5):442-448.

［3］Rapaport F，Zinovyev A，Dutreix M，et al.Classification of microarray data using gene networks［J］.BMC Bioinformatics，2007，8(1):35.

［4］Guimerà R，Sales-Pardo M，Amaral L A N.Module identification in bipartite and directed networks［J］.Physical Review E，2007，76(3):036102.

［5］Newman M E J，Girvan M.Finding and evaluating community structure in networks ［J］.Physical Review E，2004，69(2):026113.

［6］Mei J，He S，Shi G，et al.Revealing network communities through modularity maximization by a contraction-dilation method［J］.New Journal of Physics，2009，11(4):1-14.

［7］Barrett T，Troup D B，Wilhite S E，et al.NCBI GEO:mining tens of millions of expression profiles:database and tools update［J］.Nucleic Acids Research，2007，35:760-765.

［8］Rieger K E，Hong W J，Tusher V G，et al.Toxicity from radiation therapy associated with abnormal transcriptional responses to DNA damage［J］.PNAS，2004，101(17):6635-6640.

［9］Yang Y H，Dudoit S，Luu P，et al.Normalization for cDNA microarray data:a robust composite method addressing single and multiple slide systematic variation［J］.Nucleic Acids Research，2002，30(4):15-e15.

［10］Guzzi P H，Veltri P，Cannataro M.Unraveling multiple miRNA-mRNA associations through a graph-based approach［J］.ACM-BCB，2012(10):649-654.

［11］王晗.整合變化量與變化趨勢的共調(diào)控基因相似性度量［D］.長春:吉林大學計算機科學與技術(shù)學院，2008.

［12］Shi X，Su S，Long J，et al.MicroRNA-191 targets Ndeacetylase/N-sulfotransferase 1 and promotes cell growth in human gastric carcinoma cell line MGC803［J］.Acta Biochim Biophys Sin (Shanghai)，2011，43(11):849-856.

［13］Yeung K Y，Haynor D R，Ruzzo W L.Validating clustering for gene expression data［J］.Bioinformatics，2001，17(4):309-318.

［14］易東，楊夢蘇，李輝智.基因表達數(shù)據(jù)聚類分析結(jié)果的評價方法研究［J］.中國衛(wèi)生統(tǒng)計，2002，19(6):332-335.Yi Dong，Yang Meng-su，Li Hui-zhi.An nuvel method for evaluation of clustering results for gene expression data［J］.Chinese Journal of Health Statistics，2002，19(6):332-335.

［15］Datta S，Datta S.Comparisons and validation of statistical clustering techniques for microarray gene expression data［J］.Bioinformatics，2003，19(4):459-466.

［16］Belinky F，Bahir I，Stelzer G，et al.Non-redundant compendium of human ncRNA genes in GeneCards ［J］.Bioinformatics，2013，29(2):255-261.