薛 云 ，傅俊橦，李杰進，王杜齊，鄺秋華，張美珍，肖 化

(華南師范大學(xué)物理與電信工程學(xué)院，廣州510006)

隨著DNA 微陣列檢測技術(shù)的發(fā)展，產(chǎn)生的基因表達數(shù)據(jù)量呈現(xiàn)爆炸性增長態(tài)勢［1］. 如何對基因表達數(shù)據(jù)進行有效的分析，挖掘有用信息已經(jīng)成為后基因組時代的研究熱點.

作為針對基因分析的數(shù)據(jù)挖掘模型，雙聚類算法應(yīng)運而生［2－7］.雙聚類和傳統(tǒng)聚類算法的本質(zhì)區(qū)別在于，雙聚類可以對基因表達數(shù)據(jù)矩陣的行和列同時進行聚類，挖掘基因表達數(shù)據(jù)中的局部信息，從而確定在1個實驗條件子集下表現(xiàn)出相似行為的基因.這種無監(jiān)督學(xué)習(xí)技術(shù)是發(fā)現(xiàn)基因之間共表達或共調(diào)控關(guān)系、預(yù)測基因功能、分析轉(zhuǎn)錄調(diào)控、闡釋生物學(xué)通路和提供疾病機理等研究的重要手段［2－6］.

因為多個基因在關(guān)聯(lián)表達時，其各自的表達數(shù)值不需要完全相同［6］，所以基因表達水平的升降模式比基因的精確表達水平更有意義. 因此，2003年Ben-Dor 等［6］提出OPSM (Order preserving submatrices)雙聚類模型，其根據(jù)是基因表達水平在不同時間點或不同實驗條件下表現(xiàn)出相同的上調(diào)或下調(diào)趨勢. Fang［8－9］、Hochbaum［10］等分別提出了幾種近似的OPSM 模型，諸如POPSM、AOPSM，但是這些模型不能保證OPSM 的精確度. 目前眾多針對OPSM 問題的精確算法在處理巨型基因數(shù)據(jù)矩陣的問題上都存在嚴(yán)重缺陷:時間復(fù)雜度和空間復(fù)雜度過高.因此在執(zhí)行時，在保證精確度的前提下必須限制所搜索雙聚類的行列數(shù)量，這影響了部分具有重大生物意義的基因聚類(Deep OPSM 模式)的發(fā)現(xiàn).本文提出一種高效的精確性算法，用于挖掘出基因數(shù)據(jù)矩陣中的OPSM，針對Deep OPSM 問題，分析并驗證了算法解決該問題的優(yōu)點.

1 基本概念和定義

1.1 基因表達數(shù)據(jù)集

基因表達數(shù)據(jù)可用矩陣D=(dij)m×n表示:

其中包含m個基因?qū)ο?，n個樣本，元素dij表示第i個基因在第j個條件下的表達水平值.

1.2 公共子序列

序列的子序列是指從給定字符序列中任意地(不一定連續(xù))去掉若干個字符后所形成的序列.如果1個序列既是c1序列的子序列又是c2的子序列，則稱這個序列為序列c1和c2的公共子序列.

1.3 OPSM 與Deep OPSM 的定義和意義

假設(shè)矩陣D 存在1個子矩陣，對應(yīng)的行列集合為(R，C)，其中R 包含于矩陣D 的行集合中，C 包含于矩陣D 的列集合中，當(dāng)子矩陣(·)中每一行元素都呈現(xiàn)一致性變化趨勢時，此時子矩陣(R，C)就稱為OPSM 聚類［7］.圖1A 為基因g1～g6在a ～g 這7 種實驗條件下的表達值(簡稱“七因子值”)折線圖，表達值的變化趨勢雜亂無章，無規(guī)律可循. 但如果將基因g1、g2、g6在a、c、d、f、g 這5 種實驗條件下的表達值從矩陣中抽選值(簡稱“五因子值”)畫出對應(yīng)折線圖(圖1B)，這時可以清晰地看出同升同降的共性，即所挖掘的OPSM［6］.進一步對各行按元素值從小到大排序，則可以得出一致上升的模式(圖1C).

圖1 基因數(shù)據(jù)OPSM 挖掘前后的表達值圖示Figure 1 Expression data of OPSM before and after mining

行數(shù)少、列數(shù)多的OPSM 稱為Deep OPSM，這是OPSM 的一種特殊類型.從生物學(xué)的意義來說，Deep的意義在于即使基因數(shù)較少(2個基因)，只要它們在眾多的實驗環(huán)境中，都具有相同變化趨勢的表達值時，則判斷這2個基因存在緊密聯(lián)系，具有重要的生物學(xué)意義［11］.過去研究大都集中在尋找大群體基因中的OPSM 聚類，但目前生物學(xué)家更感興趣的是“小群體基因在很多實驗環(huán)境下緊密地互調(diào)節(jié)關(guān)系”［12］，因而更應(yīng)該去關(guān)注和研究一個大基因群體里面數(shù)量很少的基因，它們極有可能在某些生物進程或疾病機理上充當(dāng)著重要的角色［13］，因此Deep OPSM 在基因表達數(shù)據(jù)的挖掘中理應(yīng)受到更多研究.但大部分傳統(tǒng)挖掘算法不能很好地解決Deep OPSM 問題，當(dāng)挖掘大群體基因中OPSM 聚類時，如果設(shè)定較低的行閾值，將大量消耗計算機內(nèi)存，甚至導(dǎo)致內(nèi)存溢出.而如果為了減少計算機消耗而設(shè)定較高的行閾值，將導(dǎo)致大量的Deep OPSM 被忽略.本文提出一種新的精確算法，用于挖掘基因表達數(shù)據(jù)集，以解決這一問題，首先提出一種高效方法尋找2個序列的所有公共子序列，并結(jié)合數(shù)據(jù)結(jié)構(gòu)得到全部OPSM，也包含全部Deep OPSM.

2 針對OPSM 的精確算法

首先針對基因數(shù)據(jù)矩陣中的每2 行尋找其公共子序列，然后在整個基因數(shù)據(jù)矩陣的范圍內(nèi)，對找到的公共子序列進行支持度統(tǒng)計. 當(dāng)公共子序列的支持度超過設(shè)定閾值時，保存該子序列并輸出為目標(biāo)OPSM.

在解決時間復(fù)雜度高的問題時，利用每2個序列構(gòu)造出來的匹配矩陣，快速地找到2個序列的所有公共子序列.在解決空間復(fù)雜度高的問題時，利用了C++編程環(huán)境中的STL map，存儲產(chǎn)生的序列并統(tǒng)計相應(yīng)的支持度.算法流程如圖2 所示.

2.1 數(shù)據(jù)預(yù)處理

原始的基因表達數(shù)據(jù)矩陣中存在著部分缺失的元素，所以在數(shù)據(jù)預(yù)處理過程中必須對這部分缺失的元素進行填充.數(shù)據(jù)預(yù)處理具體操作包括:(1)將數(shù)據(jù)矩陣中缺失的元素填充為0;(2)根據(jù)基因譜中基因表達水平的大小，對各個實驗條件按照表達水平遞增進行排列;(3)用列標(biāo)識置換數(shù)據(jù)矩陣中對應(yīng)的元素.

通過以上的數(shù)據(jù)預(yù)處理操作，基因數(shù)據(jù)矩陣可以轉(zhuǎn)換成1個序列數(shù)據(jù)集(表1)

此時，OPSM 的挖掘問題轉(zhuǎn)化為序列模式挖掘問題，對基因數(shù)據(jù)矩陣的局部聚類即轉(zhuǎn)化為尋找頻繁公共子序列的問題.

2.2 尋找全部公共子序列ACS 的算法

本文提出一種高效的方法尋找2個序列的全部公共子序列即ACS (all common subsequences).

定義1 假設(shè)有2個序列c1={A B C D E F }和c2={F A B C D E }，將2個序列的其中一個序列作為行標(biāo)識，另一個序列作為列標(biāo)識.當(dāng)且僅當(dāng)元素dij對應(yīng)的行標(biāo)識與列標(biāo)識相同時，將元素dij置為相應(yīng)的行(列)標(biāo)識，其他元素置“0”，即組成1個矩陣M，稱之為匹配矩陣(圖3).

圖2 算法流程圖Figure 2 Flowchart of the algorithm

推論1 由序列c1以及序列c2生成的匹配矩陣M 中，當(dāng)且僅當(dāng)元素B 的位置在元素A 的位置的右下方(即B 行列號都大于A)時，元素B 可接在元素A 的后面構(gòu)成序列c1和c2的公共子序列.

2.3 查找和存儲ACS

根據(jù)匹配矩陣的定義和敘述，從匹配矩陣中尋找公共子序列的方法如下:

表1 數(shù)據(jù)預(yù)處理結(jié)果Table 1 Results of data preprocessing

圖3 匹配矩陣Figure 3 Matching matrix

掃描匹配矩陣中第一列，查找非零元素如A，此元素A 為第一個公共子序列的起點.

掃描下一列的元素，查找非零元素如B，判斷該元素B 的行號:

如果元素B 的行號大于上一個元素A 的行號，則元素B 可以接在原有公共子序列的后面，形成新的公共子序列.

如果元素B 的行號小于上一個元素A 的行號，則元素B 不可以接在原有公共子序列的后面，應(yīng)該生成新的公共子序列的起點.

以此類推，直到判斷完匹配矩陣的最后一列，并以樹的形式存儲查找結(jié)果. 由匹配矩陣生成的結(jié)果存儲如圖4 所示.

圖4 存儲樹結(jié)構(gòu)Figure 4 Structure of storage tree

遍歷該樹，得到2個序列的所有公共子序列.在找出每2個序列的所有公共子序列后，對結(jié)果進行存儲.對于未出現(xiàn)的新子序列，創(chuàng)建新的存儲單元;對于此前已經(jīng)出現(xiàn)過的子序列，要在其行支持度記錄上加1.因為子序列的個數(shù)多，所以必須使用一種適合快速檢索的數(shù)據(jù)結(jié)構(gòu)進行存儲.

2.4 頻繁項集統(tǒng)計方法

利用C++標(biāo)準(zhǔn)程序庫中STL map 對產(chǎn)生的公共子序列進行存儲并統(tǒng)計其支持度. 對比支持度與設(shè)定的最小支持度，將頻繁項集作為結(jié)果輸出.STL map 的底層數(shù)據(jù)結(jié)構(gòu)是紅黑樹. 紅黑樹是一種自平衡二叉查找樹，可以在O(lg n)時間內(nèi)做查找、插入和刪除，這里的n 是樹中元素的數(shù)目.運用STL map對產(chǎn)生的公共子序列進行存儲并統(tǒng)計支持度，能夠降低算法的空間和時間復(fù)雜度.

采用STL map 對生成的子序列進行存儲具有以下優(yōu)點:(1)用紅黑樹的數(shù)據(jù)結(jié)構(gòu)能快速檢索元素，可以在O(lg n)時間內(nèi)做查找和插入，而雙聚類問題本身是NP 問題，因此該數(shù)據(jù)結(jié)構(gòu)非常適用于雙聚類的精確算法;(2)可以根據(jù)生成的子序列的長度進行判斷，從而對子序列進行條件存儲.當(dāng)程序只對較長的生成子序列進行存儲時，能夠快速地找到長的頻繁項集，這對于挖掘出Deep OPSM 非常有利;(3)根據(jù)條件存儲，可以同時在多臺計算機分析同1個數(shù)據(jù)集，實現(xiàn)并行計算，提高分析處理速度，適應(yīng)大型數(shù)據(jù)矩陣的分析需求.

3 案例分析

為了驗證本文方法的有效性，利用真實的基因數(shù)據(jù)矩陣來進行實驗. 使用半乳糖酵母基因數(shù)據(jù)集中的時間序列，其中包含205個半乳糖酵母基因在20個時間點(大于2個完整的細(xì)胞周期)的表達數(shù)據(jù)［14］.它是一組在真實的實驗條件下獲得的基因表達數(shù)據(jù)，是Yee 等［15］從生物數(shù)據(jù)庫中搜索出的205條基因，它們分別歸屬于4個功能類別(表2)，這4個功能類別也被廣泛用于衡量聚類質(zhì)量的外部標(biāo)準(zhǔn).

系統(tǒng)開發(fā)平臺參數(shù)如下:

處理器:Intel(R)Core(TM)i3 CPU M 380 @2.53 GHz;內(nèi)存:4G;運行速度2 533.3 MHz;計算機系統(tǒng)為Windows7;運行軟件為VS2010;實驗編程語言為C++.

表2 數(shù)據(jù)集的4個功能類別Table 2 Four functions of the genes in this dataset

3.1 OPSMs 可視化分析

使用本文提出的算法挖掘出全部OPSMs，并把4個典型OPSMs 進行可視化操作，繪制成折線圖(圖5)，直觀地發(fā)現(xiàn)其共同特點. 橫坐標(biāo)為對應(yīng)OPSMs 聚類的列號，縱坐標(biāo)為對應(yīng)的聚類元素在基因表達數(shù)據(jù)矩陣上的表達值. 每一條折線所對應(yīng)的聚類基因，其表達值在原矩陣中呈現(xiàn)同增同減的一致變化趨勢，即OPSM.

圖5 隨機挑選實驗中找到的4個OPSMs 折線圖Figure 5 Broken line charts of OPSMs found in randomly picked experiments

3.2 算法性能比較

對文獻［7］算法和本文算法所挖掘出的OPSMs進行可視化操作，分別繪制成折線圖進行比較.行閾值均設(shè)置為20，本文算法的列閾值設(shè)置為6 時，由于剪掉的序列模式較多，所以2 種算法均能在內(nèi)存負(fù)擔(dān)能承受的范圍內(nèi)挖掘出相應(yīng)的結(jié)果(圖6). 均比較第33 790個OPSMs.橫坐標(biāo)為對應(yīng)OPSMs 聚類的列號，縱坐標(biāo)為對應(yīng)的聚類元素在基因表達數(shù)據(jù)矩陣上的表達值.2 種算法對應(yīng)結(jié)果中每一條折線對應(yīng)的聚類基因，其表達值在原矩陣中呈現(xiàn)一致的同增同減變化趨勢，即為OPSM.

為了挖掘出Deep OPSM 而把2個算法的行閾值降低為2 時，文獻［7］算法在運行過程中內(nèi)存占用量極大，導(dǎo)致結(jié)果無法挖掘出來，普通的服務(wù)器根本無法承受其內(nèi)存負(fù)擔(dān). 本文算法在同等行閾值條件下，運行過程中占用的內(nèi)存比文獻［7］算法占用的內(nèi)存小，說明本文算法更適合挖掘Deep OPSM.如圖7 所示，統(tǒng)計行閾值α 分別為2、3、4 時，可挖掘出列數(shù)分別為10 ～15 的Deep OPSM 的數(shù)量. 在相同行閾值下，OPSM 數(shù)量隨著列數(shù)的增加而減少;在列數(shù)相同時，OPSM 數(shù)量隨著行閾值的增加而急劇減少. 而行數(shù)較少、列數(shù)較多的OPSM 往往更具有生物學(xué)意義，本算法能有效地將Deep OPSMs 挖掘出來.

3.3 GO 分析

GO(Gene Ontology)［16］注釋被用來檢驗本文算法生成的OPSM 聚類的生物學(xué)意義，對產(chǎn)生的聚類結(jié)果作真實性驗證. GO 是基因本體聯(lián)合會(Gene Onotology Consortium)所建立的數(shù)據(jù)庫，旨在建立一個適用于各種物種的、對基因和蛋白質(zhì)功能進行限定和描述的、并能隨著研究不斷深入而更新的語言詞匯標(biāo)準(zhǔn).

圖6 算法性能比較Figure 6 Performance comparisons of two algorithms

圖7 低閾值Deep OPSM 統(tǒng)計圖Figure 7 Statistical chart of Deep OPSMs based on low threshold

從本文算法產(chǎn)生的OPSM 聚類中挑選出部分聚類，進行GO 注釋實驗.結(jié)果證明，實驗產(chǎn)生的OPSM聚類，能夠?qū)?yīng)某些基因關(guān)聯(lián)調(diào)控的生物特征.列出部分GO 注釋的查找結(jié)果如表3 所示. 其中當(dāng)Pvalue 值較低時，說明了該簇基因在1個或者多個GO 項目中具有某些特別的生物學(xué)意義. 部分P-value 值小于1.0 ×10－25的GO 項目，表明本文算法能夠找到包含關(guān)聯(lián)調(diào)控信息的基因簇，從生物意義上證明了本文方法的可行性.

表3 GO 分析部分結(jié)果Table 3 Part of GO analysis results

4 結(jié)論

OPSM 模型作為一種基于模式的雙聚類方法，在分析基因數(shù)據(jù)矩陣等方面被廣泛的應(yīng)用. 在一個OPSM 聚類中，形成聚類的若干基因在特定的條件子集下有一致的表達模式. 這種關(guān)聯(lián)的共同表達隱含著基因的關(guān)聯(lián)調(diào)控，所以在基因數(shù)據(jù)矩陣上進行的雙聚類分析有極大的生物學(xué)意義. 本文根據(jù)OPSM 模型，建立了一種快速有效的精確性尋找方法，來挖掘分散在基因數(shù)據(jù)矩陣中的OPSM 聚類.結(jié)果證明，該方法能夠快速地找到符合條件的OPSM聚類，并且能夠通過條件存儲，針對長頻繁模式進行尋找分析，挖掘出更具生物學(xué)意義的Deep OPSM 聚類.此外，通過條件存儲，可以在多臺計算機上實現(xiàn)并行計算，提高分析處理速度，適應(yīng)大型數(shù)據(jù)矩陣的分析需求.最后從生物學(xué)(以半乳糖酵母基因數(shù)據(jù)中的時間序列)的角度，驗證了該方法的可行性.

［1］蔡莉，郭紅. 一種改進的基因表達數(shù)據(jù)雙聚類算法［J］.福州大學(xué)學(xué)報:自然科學(xué)版，2010，38(1):41－47.Cai L，Guo H. An improved biclustering a lgorithm for gene expression data［J］. Journal of Fuzhou University:Natural Science，2010，38(1):41－47.

［2］Lazzeronil L，Owen A. Plaid models for gene expression data［J］. Statistica Sinica，2000，12(1):61－86.

［3］Andrew C T，Prokopyev O A. Solving the order-preserving submatrix problem via integer programming［J］. Informs Journal on Computing，2010，22(3):387－400.

［4］Liu J Z，Wei W. OP-cluster:Clustering by tendency in high dimensional space［C］∥3rd IEEE International Conference on Data Mining. Melbourne，USA，2003:187－194.

［5］Yun X，Li Y T，Deng W J，et al. Mining order-preserving submatrices based on frequent sequential pattern mining［C］∥Health Information Science，HIS 2014.Shenzhen，China，2014:184－193.

［6］Ben-Dor A，Chor B，Karp R，et al. Discovering local structure in gene expression data:The order-preserving submatrix problem［J］. Journal of Computational Biology，2003，10(3/4):373－384.

［7］Cheung L，Kevin Y Y，Cheund D W，et al. On mining micro-array data by order-preserving submatrix［J］. International Journal of Bioinformatics Research A，2007，3(1):42－64.

［8］Fang Q. Effective discovery of order-preserving submatrices［D］. Hong Kong:Hong Kong University of Science and Technology，2012.

［9］Fang Q，Wilfred N G，F(xiàn)eng J L，et al. Mining order-preserving submatrices from probabilistic matrices［J］. ACM Transactions on Database Systems，2014，39(1):1－44.

［10］Hochbaum D S，Levin A. Approximation algorithms for a minimization variant of the order-preserving submatrices and for biclustering problems［J］. ACM Transactions on Algorithms，2013，9(2):1－12

［11］Byron J G，Grifiith O L，Ester M，et al. On the deep order-preserving submatrix problem:A best effort approach［J］. IEEE Transactions on Knowledge and Data Engineering，2012，24(2):309－325.

［12］Williamson M P，Sutcliffe M J. Protein-protein interactions［J］. Biochemical Society Transactions，2010，38(4):875－878.

［13］Albert R. Scale-free networks in cell biology［J］. Journal of Cell Science，2005，118(21):4947－4957.

［14］Ideker T，Thorsson V，Ranish J A，et al. Integrated genomic and proteomic analyses of a systematically perturbed metabolic network［J］. Science，2001，292(5518):929－934.

［15］Yee Y K，Medvedovic M，Bumgarner R E. Clustering gene-expression data with repeated measurements［J］.Genome Biology，2003，4(5):1－17.

［16］David M，Brun C，Remy E，et al. GOToolBox:Functional analysis of gene datasets based on gene ontology［J］. Genome Biology，2005，5(12):1－8.