肖碧玉，李先彬，沈良忠，劉文斌

(溫州大學物理與電子信息工程學院，浙江溫州，325035)

系統(tǒng)生物研究表明:生命體的生長、發(fā)育、衰老過程以及疾病的發(fā)生與基因之間相互作用的變化密切相關。隨著以微陣列技術為代表的各種高通量技術的飛速發(fā)展，人們可以同時觀察到一個細胞中成千上萬個基因的活動狀況。如何利用這些數(shù)據(jù)挖掘出其中潛在變化，對于揭示衰老及疾病的發(fā)生發(fā)展，具有重要的生物學意義。因此，基于基因表達譜數(shù)據(jù)的差異分析成為系統(tǒng)生物學的一個研究熱點。基因差異表達的分析大致可以分為以下三個層次:

(1)生物學中，有些關鍵基因表達水平的高低往往是導致疾病產(chǎn)生的關鍵因素。因此，最早的差異分析主要集中在單個基因表達水平的變化，如Jacob等發(fā)現(xiàn)小鼠衰老期各組織的基因表達水平變化存在較大的差異，可歸類為神經(jīng)組織、血管組織、類固醇反應組織三種衰老類型［1］

(2)基因之間往往是通過相互作用來實現(xiàn)某種功能，因而挖掘基因作用關系的變化能夠揭示基因(簇)與功能之間的關系。Remondini通過比對基因時間序列網(wǎng)絡的度的分布，發(fā)現(xiàn)c-myc致癌基因的活性與基因的網(wǎng)絡結(jié)構(gòu)有密切關系［2］。Voy等通過研究網(wǎng)絡邊的變化，確定小鼠受輻射影響的基因簇［3］。Oldham等比對人類和黑猩猩差異網(wǎng)絡的拓撲重疊，挖掘與生物驅(qū)動進化有關的基因集合［4］。Omar Odiba等在差異網(wǎng)絡分析中將網(wǎng)絡中介性及中心性結(jié)合起來，挖掘出了生物標記(Biomarker)［5］。Zhang基于條件概率提出差異相關子網(wǎng)絡(DDN)，檢測兩個差異轉(zhuǎn)錄網(wǎng)絡的拓撲變化顯著性［6］。Southworth等構(gòu)建基因帶權差異網(wǎng)絡，挖掘與小鼠衰老密切相關的基因模塊［7］。

(3)雖然基因模塊往往與某種功能相關，但是這種功能往往需要其他基因模塊的協(xié)調(diào)與配合，基因模塊之間關系的變化也往往與功能的變化密切相關。Tesson等提出了DiffCoEx算法，可以同時挖掘基因模塊及基因模塊間的變化［8］。此外，在差異表達分析方面，F(xiàn)ang研究了二組樣本中表達水平高低的變化，提出了三種基因表達譜差異模式及其關系，這種模式特征的研究對于差異模式的挖掘具有重要的指導意義［9-10］在生物網(wǎng)絡中，節(jié)點往往通過與局部結(jié)構(gòu)作用實現(xiàn)功能，其與局部結(jié)構(gòu)的關系發(fā)生變化，可能會影響其功能。因此，如何理解和刻畫網(wǎng)絡節(jié)點的局部結(jié)構(gòu)特征及變化，對于認識生物系統(tǒng)的結(jié)構(gòu)和功能的關系具有非常重要的意義。點的聚類系數(shù)是比較經(jīng)典的網(wǎng)絡節(jié)點局部結(jié)構(gòu)描述測度，人們利用它做出了許多有意義的結(jié)果［16］。2003年，Przulj提出利用圖元及圖元向量刻畫網(wǎng)絡中節(jié)點鄰域關系，并且基于它們研究了生物網(wǎng)絡與隨機網(wǎng)絡的拓撲結(jié)構(gòu)差異［11-12］、癌癥相關基因［13］、生物網(wǎng)絡的進化樹［14］等。本文利用仿真實驗的方法分析和比較圖元向量和點的聚類系數(shù)兩個測度的性能，并且分別利用它們設計算法挖掘差異網(wǎng)絡中模塊化變化的節(jié)點簇，最后利用AGEMAP數(shù)據(jù)庫中小鼠的基因表達數(shù)據(jù)進行實驗，分析比較基于圖元向量和基于點的聚類系數(shù)挖掘小鼠基因時序差異網(wǎng)絡的差異模塊與小鼠衰老的關系。

1 基本概念

點的聚類系數(shù)定義為:

其中dv表示節(jié)點v的度，這dv個節(jié)點之間實際存在的邊數(shù)目ev表示節(jié)點v及其相鄰節(jié)點之間的邊數(shù)。圖元就是包含一定節(jié)點的非同構(gòu)子圖，其非同構(gòu)位置唯一標號構(gòu)成的向量稱為圖元向量。圖1列出了包含2、3、4個節(jié)點的圖元和其15維圖元向量。由于隨著圖元節(jié)點規(guī)模的增大，圖元向量的規(guī)模及計算的復雜度將指數(shù)級增加，本文采用圖1中的15維的圖元向量來研究差異網(wǎng)絡。此外，由于標號大的圖元向量可能包括標號小的維度信息(如標號14與3，7與2)，Milenkovi按照標號大的圖元向量包括標號小的圖元向量的個數(shù)對每一維圖元向量定義一個權值 wk(k=0，…，14)。

圖1 圖元結(jié)構(gòu)Fig.1 Graphlets and orbits

給定兩個差異網(wǎng)絡 G(V，E)和 G(V，E＇)，具有相同的節(jié)點集不同邊集，我們定義節(jié)點v∈V的圖元向量差異度為:

其中vk和 v′k分別表示節(jié)點 v在兩個差異網(wǎng)絡G(V，E)和 G(V，E＇)中的第 k維圖元向量。wk表示圖元向量第k維的圖元向量的權值，度是反映網(wǎng)絡節(jié)點鄰接結(jié)構(gòu)的一個非常重要的指標，結(jié)合圖元向量和點的聚類系數(shù)兩個測度，我們定義節(jié)點v∈V的兩種局部結(jié)構(gòu)差異度為:

其中dv，d′v分別表示節(jié)點v在兩個差異網(wǎng)絡G(V，E)和 G(V，E＇)中的度，cv，c′v分別表示節(jié)點 v在兩個差異網(wǎng)絡 G(V，E)和 G(V，E＇)中的點的聚類系數(shù)。直觀上，網(wǎng)絡中同一個簇的節(jié)點變化，應該具有一定的相似性，即對于一個簇中的節(jié)點u，v∈V＇，他們的局部結(jié)構(gòu)差異度可能基本相當。因此，我們給出衡量節(jié)點局部結(jié)構(gòu)變化差距的定義如下:

2 比較圖元向量和點的聚類系數(shù)的性能

為了定量證明圖元向量和點的聚類系數(shù)兩個測度的性能，本文采用計算機產(chǎn)生隨機網(wǎng)絡作仿真數(shù)據(jù)，該實驗方法已被廣泛采用［15-16］本仿真實驗產(chǎn)生的每一個隨機網(wǎng)絡，包含128個節(jié)點，每個節(jié)點的度p=16，這些點被劃分成4個社團，每個點隨機的與所在社團內(nèi)部的pin點相連，與其它的pout=p－pin個點相連。顯然，隨著pout的增大，網(wǎng)絡中的社團結(jié)構(gòu)越來越模糊;反之，社團結(jié)構(gòu)則越來越明顯。一共產(chǎn)生了3組隨機網(wǎng)絡，分別對應 pout=1、5、8，每組共100個網(wǎng)絡。為了驗證比較圖元向量和點的聚類系數(shù)的性能，我們對這每一組的100個網(wǎng)絡隨機加(減)10、20、30、40、50、60、70、80、90、100 條邊進行分析。

日前，由于實驗平臺、實驗技術等的限制，生物數(shù)據(jù)還存在一定的誤差，這對差異信息的研究帶來了一定的困難，所以，現(xiàn)在人們主要是對較明顯的差異信息進行研究，而往往把較小的變化當做噪聲處理。本實驗以加/減少量10(1%)的邊仿真噪聲，以加/減較多的邊100(10%)仿真差異變化，用噪聲魯棒性衡量測度容忍噪聲的能力，也就是計算在有噪聲干擾下測度度量時的變化量;用差異靈敏性度量測度發(fā)現(xiàn)差異變化的能力，也就是計算發(fā)生差異變化時測度發(fā)生的變化量。

為了分析和比較圖元向量、點的聚類系數(shù)的噪聲魯棒性和差異靈敏度，我們分別對3類網(wǎng)絡加(減)10條邊和100條邊，各節(jié)點的點的聚類系數(shù)和圖元向量的變化量做統(tǒng)計(見圖2)。圖中右側(cè)3個圖中，3類網(wǎng)絡中加(減)100條邊時圖元向量的變化量都高于加(減)10條邊時的變化，圖元向量測度能在容忍噪聲的同時，有效提取差異信息，而圖中左側(cè)3類網(wǎng)絡中點的聚類系數(shù)測度則不然?？梢姡鄬c的聚類系數(shù)，圖元向量的差異靈敏度和噪聲魯棒性比較好。

圖2 比較點的聚類系數(shù)和圖元向量的噪聲魯棒性、差異靈敏度Fig.2 The comparison of robustness and difference sensitivity of clustering coefficient and orbits under various noises

點的聚類系數(shù)(圖元向量)在各類網(wǎng)絡中的噪聲魯棒性、差異靈敏度是否相同呢?我們統(tǒng)計3類網(wǎng)絡加(減)10條邊、100條邊時，各節(jié)點的點聚類系數(shù)、圖元向量變化(見圖3)。圖3右側(cè)兩個圖中，在加(減)10條邊的噪聲魯棒性研究，和加(減)100條邊的差異靈敏性分析，圖元向量在3類網(wǎng)絡(Pout=1、5、8)中都相似。從圖3左側(cè)兩個圖中不難看出，在加(減)10條邊的噪聲魯棒性研究，和加(減)100條邊的差異靈敏性分析，各節(jié)點的點聚類系數(shù)的變化，在Pout=1的網(wǎng)絡的噪聲魯棒性最差、差靈敏度最好，在Pout=5的網(wǎng)絡噪聲魯棒性次之、差靈敏度也次之，在Pout=8的網(wǎng)絡噪聲魯棒性最好、差靈敏度最差，點聚類系數(shù)在3類網(wǎng)絡(Pout=1、5、8)中的噪聲魯棒性和差異靈敏度不同。相對點的聚類系數(shù)，圖元向量的適用性較強，應用范圍較廣。

我們統(tǒng)計 3 類網(wǎng)絡加(減)邊 10、20、30、40、50、60、70、80、90、100 條，局部結(jié)構(gòu)變化的節(jié)點的數(shù)量(見圖4)。分析圖4不難發(fā)現(xiàn)，在3類網(wǎng)絡中，各種加(減)邊時，圖元向量變化的節(jié)點幾乎是整個網(wǎng)絡128個節(jié)點，而點的聚類系數(shù)變化的節(jié)點則隨加(減)邊的數(shù)量的增加而增多，可見，相對點的聚類系數(shù)，圖元向量更能細致反應節(jié)點局部結(jié)構(gòu)的影響。

綜上可得，相對點的聚類系數(shù)，圖元向量具更好的噪聲魯棒性和差異靈敏度，有較強的適用能力，但是在時間復雜性方面較差。點的聚類系數(shù)和圖元向量各具優(yōu)缺點，我們可以根據(jù)需要酌情選擇。

圖3 比較點的聚類系數(shù)和圖元向量的適用性Fig.3 The comparison of clustering coefficient and orbits by using

圖4 比較點的聚類系數(shù)和圖元向量變化的節(jié)點數(shù)Fig.4 The comparison of clustering coefficient and orbits with the number of nodes Changed adding or cutting

3 實驗

3.1 材料

本文生物數(shù)據(jù)來自AGEMAP數(shù)據(jù)庫(http://cmgm.stanford.edu/～ kimlab/aging_mouse/)，其中包括 C57BL6小鼠16個組織的1、6、16、24個月的8 932個基因表達數(shù)據(jù)，由于其中4個組織數(shù)據(jù)存在較大噪聲和不完整性，本文使用Adrenal Glands(1)、Cerebellum(2)、Cerebrum(3)、Eye(4)、Gonads(5)、Heart(6)、Hippocampus(7)、kidney(8)、Lung(9)、Muscle(10)、Spinal Cord(11)、Thymus(12))等 12 個組織的數(shù)據(jù)。通常認為16月至24月為小鼠的衰老期，我們利用這兩個月基因表達數(shù)據(jù)研究差異簇與小鼠衰老的關系。共表達網(wǎng)絡是一個無向圖，每個節(jié)點表示一個基因，每條邊表示兩個基因的共表達關系。共表達關系通過pearson相關系數(shù)r確定

然后將Pearson系數(shù)r轉(zhuǎn)化另一變量

其中，n表示計算這兩個基因相關系數(shù)時所用的數(shù)據(jù)點的個數(shù)。轉(zhuǎn)化后的r＇值服從自由度為n－2的t分布。如果兩個基因之間的r＇值大于設定的p-value對應的t分布表，則就在這兩個基因之間加一條邊;否則在這兩個基因之間就不加邊。除Eye組織的p-value取值為1E-04，Spinal Cord組織為1E-07、其他各組織均取5E-06構(gòu)建網(wǎng)絡。David數(shù)據(jù)庫可以對小鼠基因簇的功能進行分析，如果一個基因簇有50%以上的基因顯著的共享一個或多個GO項，通常認為該基因簇具有生物學意義。本文GO項的P-value取0.05。Southworth等驗證了401條與小鼠衰老相關的GO項［8］，因此，同其分析一個基因簇顯著富集的GO項就可以判斷其是否與衰老有關。下面我們通過圖元向量和點的聚類系數(shù)二種測度，分別設計算法挖掘16和24月的基因表達數(shù)據(jù)中變化的基因簇，并對其功能進行分析。

3.2 挖掘差異基因簇

基于圖元向量挖掘差異基因簇，首先，提取局部結(jié)構(gòu)發(fā)生變化的差異基因;在16月(24月)基因共表達網(wǎng)絡中，取Dv1大于等于5提取16月(24月)的變化的基因。然后，根據(jù)提取出的基因多少和D1v分布，利用Luv分別在16月和24月基因共表達網(wǎng)絡中進行K-means聚類，最后，對聚類簇進行GO注釋和分析。基于點的聚類系數(shù)挖掘差異基因簇方法類似。

3.3 實驗結(jié)果分析

基于圖元向量和點的聚類系數(shù)挖掘差異基因簇的實驗結(jié)果如圖5。從圖中可以看出，基于圖元向量和點的聚類系數(shù)挖掘的差異基因簇，幾乎每個基因簇都顯著的富集一個或多個GO項，并且大部分的基因簇都與衰老相關。二種方法挖掘聚類簇顯著富集的GO項幾乎相似，主要相似GO項見表1，這些功能的差異與衰老具有密切關系。有些基因簇按照Southworth驗證的401條衰老GO項與衰老無關，這是由于目前有關衰老相關的GO項還不完善，這方面的研究工作還在不斷進行中，如Chunxiao Fu［17］、Jamie L.Barger［18］等也驗證其它一些與小鼠衰老有關的GO項。

表1 二種方法挖掘的相似GO項Table 1 Similar GO terms from two algorithms

圖5 基于圖元向量和點的聚類系數(shù)的實驗結(jié)果Fig.5 Experimental results of the two algorithms based on Graphlet and Clustering Coefficient

二種方法取相同的閾值提取差異基因，同組織中基于圖元向量提取的差異基因和基于點的聚類系數(shù)提取的差異基因大部分相同，且基于圖元向量提取的差異基因相對較多(見表2)。

表2 比較各組織基于圖元向量和點聚類系數(shù)提取的基因數(shù)Table 2 Compareing the number of genes mined by the two algorithms

這與前面圖元向量和點的聚類系數(shù)的仿真實驗結(jié)論一致，圖元向量更能反映局部結(jié)構(gòu)的信息。相對點的聚類系數(shù)，圖元向量的噪聲魯棒性和差異靈敏度較好，而基于圖元向量和點的聚類系數(shù)的實驗卻得到了相似較好的實驗結(jié)果，從而可知AGEMAP數(shù)據(jù)庫提供的數(shù)據(jù)質(zhì)量較高，而這點在其他的實驗中也得到證明［1，7］。此外，我們發(fā)現(xiàn)基于圖元向量和點的聚類系數(shù)挖掘的有些基因簇在David數(shù)據(jù)庫中不能注釋到與衰老相關的 GO項，但是在AGEMAP數(shù)據(jù)庫中發(fā)現(xiàn)它們具有這些功能。如Thymus組織利用基于圖元向量挖掘的一個基因簇(其中包括196個基因)，在David數(shù)據(jù)庫中注釋衰老GO:0005488～binding項的富集度為52.6%，顯著性為0.027。此簇中的Mm.20935基因不含GO:0005488～binding項功能，但在AGEMAP提供的數(shù)據(jù)中則顯示此組織基因 Mm.20935具有 GO:0005488～binding功能。另外，Thymus組織利用基于點的聚類系數(shù)挖掘的一個基因簇(其中包括90個基因)在David數(shù)據(jù)庫中注釋衰老GO:0005488～binding項的富集度為60%，顯著性為0.013。此簇中的 Mm.383175基因不含 GO:0005488～binding項功能，但在AGEMAP提供的數(shù)據(jù)中則顯示此組織基因Mm.383175具有GO:0005488～binding功能。因此，本文那些與衰老無關的基因(簇)存在潛在的小鼠衰老研究價值。

4 小結(jié)

生物分子往往是通過與局部結(jié)構(gòu)相互作用發(fā)揮功能，其與局部結(jié)構(gòu)關系發(fā)生變化可能引起其功能變化，所以基于局部結(jié)構(gòu)測度進行差異研究對于認識生命的進化、發(fā)育、衰老過程及疾病的產(chǎn)生等生物問題具有重要的意義。本文首先分析和比較了圖元向量和點的聚類系數(shù)兩個測度的性能，并基于圖元向量和點的聚類系數(shù)分別設計了挖掘模塊化變化簇的算法。利用AGEMAP數(shù)據(jù)庫中小鼠12個組織的基因表達數(shù)據(jù)進行實驗，本文設計的二個算法挖掘的差異簇都顯著的富集于一些GO條目，而且其中大部分都與衰老有關。

圖元向量和點的聚類系數(shù)作為刻畫網(wǎng)絡局部拓撲結(jié)構(gòu)的二種參數(shù)，各具優(yōu)缺點，在差異分析中具有廣闊的應用前景。本文是基于局部結(jié)構(gòu)變化挖掘模塊變化基因簇，而基于模塊基因簇識別變化模塊基因簇將是本文后續(xù)工作。

References)

［1］ Jacob M.Zahn，Suresh Poosala，Art.B Owen.AGEM AP:A Gene Expression Database for Aging in Mice［J］.PLOS Genetics，2007，3(11):e201.

［2］ Remondini D，O，Connell B，Intrator N，Sedivy JM，Neretti N，Castellani GC，Cooper LN.Targeting c-Myc-activated genes with a correlation method:Detection of global changes in large gene expression network dynamics［J］.Proc Natl Acad Sci U S A，2005，102(19):6902 －6906.

［3］ Voy BH，Scharff JA，Perkins AD，Saxton AM，Borate B，Chesler EJ，Branstetter LK，Langston MA.Extracting gene networks for low-dose radiation using graph theoretical algorithms［J］.PLoS computational biology，2006，2(7):757－768.

［4］ Oldham MichaelC，Horvath Steve，Geschwind DanielH.Conserva-tion and evolution of gene coexpression networks in human and chimpanzee brains［J］.Proceedings of the National Academy of Sciences，2006，103(47):17973 －17978.

［5］ Omar Odibat. RankingDifferentialGenesin Co-expression Networks［J］.Proceedings of the 2nd ACM Conference on Bioinformatics，Computational Biology and Biomedicine，2011，350－354.

［6］ Zhang B，Li H，Riggins RB，zhan M，Xuan J，Zhang Z，Hoffman EP，Charke R，Wang Y.Differential dependency network analysis to identify condition-specific topological changes in biological networks［J］.Bioinformatics，2009，25(4):526 － 532.

［7］ Lucinda K.Southworth，Art B.Owen，Stuart K.Kim.Aging Mice Show a Decreasing Correlation of Gene Expression within Genetic Modules［J］.PLOS Genetics，2009，(5):e1000776.

［8］ Bruno M Tesson，Rainer Breitling and Ritsert C Jansen.Diffcoex:a simple and sensitive method to find differentially coexpressed gene modules［J］.BMC Bioinformatics，2010，11(1):497.

［9］ Gang Fang，Gaurav Pandey.Mining Low-Support Discriminative Patterns from Dense and High-Dimensional Data.Knowledge and Data Engineering［J］.IEEE Transactions，2012，24(2):279 －294.

［10］ Gang Fang，Rui Kuang，Gauraw Pandey，Michael Steinbach，Chad L.Myers，Vipin Kumar.Subspace differential coexpression analysis:Problem Definition and a General Approach［J］.Pacific Symposium on Biocomputing，2010，15:145 －156.

［11］ Natasa Przulj.Biological network comparison using graphlet degree distribution［J］.Cancer Inform，2008，6:257 －273.

［12］ Przulj N，Corneil D G，Jurisica I.Modeling Interactome:Scale-Free or Geometric［J］.Bioinformatics，2004，20(18):3508 －3515.

［13］ TijanaMilenkovicand NatasaPrzulj. UncoveringBiological Network Function via Graphlet Degree Signatures［J］.Cancer Inform，2008，6:257－273.

［14］ Oleksii Kuchaiev，Tijana Milenkovic，Vesna Memisevic，Wayne Hayes，Natasa Przulj.Topological network alignment uncovers biological function and phylogeny［J］.Journal of the Royal Society Interface，2011，5(17):1341 －1354.

［15］ Yang B.，Liu D.Y.，Liu J.M..Complex Network Clustering Algorithms［J］.Journal of Software，2009，20(1):54 － 66.

［16］ Radicchi F.，Castellano C.，Cecconi F.Defining and identifying communities in networks［J］.PNAS，2004，101(9):2658 －2663.

［17］ Fu Chunxiao，Hickey Morgen，Morrison Melissa，McCarter Roger，Han Eun-Soo.Tissue specific and non-specific changes in gene expression by aging and by early stage CR［J］.Mechanisms of Ageing and Development，2006，(127):905 －916.

［18］ Barger Jamie L，Kayo Tsuyoshi，Vann James M，Arias Edward B，Wang Jelai，Hacker Timothy A，Wang Ying，Raederstorff Daniel，Morrow Jason D，Leeuwenburgh Christiaan，Allison David B，Saupe Kurt W，Cartee Gregory D，Weindruch Richard，Prolla Tomas A.A Low Dose ofDietary ResveratrolPartially Mimics Caloric Restriction and Retards Aging Parameters in Mice［J］.PLoS ONE，2008，(3):e2264.