王 希, 潘 理, 蔣軍強(qiáng), 楊 勃, 李文彬
(湖南理工學(xué)院 信息科學(xué)與工程學(xué)院, 湖南 岳陽 414006)
蛋白質(zhì)是生物生存的基礎(chǔ), 也是構(gòu)成生物細(xì)胞和組織并維持生命活動(dòng)的基本物質(zhì).某些蛋白質(zhì)的缺失會(huì)導(dǎo)致生物體的疾病或影響生長, 嚴(yán)重的甚至?xí)?dǎo)致生物體的滅絕.根據(jù)蛋白質(zhì)對生物體生存繁衍的重要性, 蛋白質(zhì)可分為兩大類: 關(guān)鍵蛋白質(zhì)和非關(guān)鍵蛋白質(zhì).識別關(guān)鍵蛋白質(zhì)對于了解細(xì)胞的生存和發(fā)育的最低要求、分析疾病的原因并發(fā)現(xiàn)藥物作用的機(jī)制等都具有重要意義, 現(xiàn)已成為生物信息學(xué)研究的熱點(diǎn).
目前識別關(guān)鍵蛋白質(zhì)的方法主要分為兩類: 網(wǎng)絡(luò)拓?fù)浞椒╗1]和機(jī)器學(xué)習(xí)方法[2].網(wǎng)絡(luò)拓?fù)浞椒ㄍㄟ^蛋白質(zhì)相互作用進(jìn)行蛋白質(zhì)相互作用(Protein-Protein Interaction, PPI)網(wǎng)絡(luò)構(gòu)建, 隨后對PPI 網(wǎng)絡(luò)中各蛋白質(zhì)節(jié)點(diǎn)進(jìn)行拓?fù)渲行男耘判蛞詫?shí)現(xiàn)關(guān)鍵蛋白質(zhì)識別.機(jī)器學(xué)習(xí)方法以PPI 網(wǎng)絡(luò)中各蛋白質(zhì)節(jié)點(diǎn)的拓?fù)涮卣饕约芭c之相關(guān)的其它生物信息(如基因表達(dá)水平, 亞細(xì)胞位置等) 作為訓(xùn)練樣本的特征值, 隨后進(jìn)行模型訓(xùn)練以實(shí)現(xiàn)關(guān)鍵蛋白質(zhì)識別.
靜態(tài)PPI 網(wǎng)絡(luò)只包含了蛋白質(zhì)之間的相互作用, 但用于構(gòu)建靜態(tài)PPI 網(wǎng)絡(luò)的蛋白質(zhì)相互作用數(shù)據(jù)存在大量的假陰性與假陽性數(shù)據(jù).為了更為精準(zhǔn)地識別關(guān)鍵蛋白質(zhì), 研究人員提出了多信息融合的活躍PPI 網(wǎng)絡(luò)的構(gòu)建方法: 固定閾值方法[3]和3Sigma 閾值方法[4].兩種方法均采用蛋白質(zhì)相互作用信息與基因表達(dá)數(shù)據(jù)融合以提高PPI 網(wǎng)絡(luò)構(gòu)建數(shù)據(jù)的可靠性.固定閾值方法對靜態(tài)PPI 網(wǎng)絡(luò)中蛋白質(zhì)設(shè)置統(tǒng)一的活躍閾值, 若蛋白質(zhì)在某時(shí)點(diǎn)的基因表達(dá)水平高于閾值, 則蛋白質(zhì)在這個(gè)時(shí)點(diǎn)處于活躍狀態(tài), 可與其它在此時(shí)點(diǎn)同樣活躍的蛋白質(zhì)進(jìn)行相互作用; 否則, 該蛋白質(zhì)在此時(shí)點(diǎn)為非活躍狀態(tài), 由此對PPI 網(wǎng)絡(luò)構(gòu)建數(shù)據(jù)進(jìn)行篩選.但該方法由于活躍閾值設(shè)置過于絕對, 在蛋白質(zhì)活躍性判斷中, 部分蛋白質(zhì)在所有時(shí)點(diǎn)都處于非活躍狀態(tài), 與這部分蛋白質(zhì)相關(guān)的相互作用也就無法在PPI 網(wǎng)絡(luò)中得以體現(xiàn).3Sigma 閾值方法根據(jù)各蛋白質(zhì)的基因表達(dá)數(shù)據(jù)特性對各蛋白質(zhì)進(jìn)行單獨(dú)的活躍閾值設(shè)定以比較判別各時(shí)點(diǎn)的活躍狀態(tài), 隨后進(jìn)行PPI 網(wǎng)絡(luò)的數(shù)據(jù)篩選以及構(gòu)建.雖然活躍PPI 網(wǎng)絡(luò)對用于構(gòu)建網(wǎng)絡(luò)的蛋白質(zhì)相互作用進(jìn)行了有效篩選, 且關(guān)鍵蛋白質(zhì)的識別率也得到了提升, 但整個(gè)網(wǎng)絡(luò)仍是一個(gè)單層的PPI 網(wǎng)絡(luò), 沒有充分挖掘基因表達(dá)水平數(shù)據(jù)的時(shí)序特性, 蛋白質(zhì)之間的動(dòng)態(tài)特性并不能在網(wǎng)絡(luò)中體現(xiàn), 而現(xiàn)實(shí)中的蛋白質(zhì)活動(dòng)卻是一個(gè)隨著時(shí)間變化而變化的過程.
本文提出一種基于3Sigma 閾值的多層PPI 網(wǎng)絡(luò)構(gòu)建方法, 該方法將靜態(tài)PPI 網(wǎng)絡(luò)中的相互作用映射到多個(gè)觀測時(shí)點(diǎn)網(wǎng)絡(luò)的子網(wǎng)中, 不僅能對網(wǎng)絡(luò)構(gòu)建數(shù)據(jù)進(jìn)行有效篩選,提高了網(wǎng)絡(luò)的可靠性, 并且整個(gè)多層網(wǎng)絡(luò)使得蛋白質(zhì)相互作用之間的動(dòng)態(tài)變化特性也得以體現(xiàn),提高了網(wǎng)絡(luò)的合理性.在關(guān)鍵蛋白質(zhì)識別方面, 多層PPI 網(wǎng)絡(luò)的識別結(jié)果也優(yōu)于現(xiàn)有單層PPI 網(wǎng)絡(luò).
一個(gè)靜態(tài)PPI 網(wǎng)絡(luò)可以描述為一個(gè)無向圖.令無向圖G0=(V0,E0)為一個(gè)靜態(tài)PPI 網(wǎng)絡(luò), 其中
(1)V0={v1,v2,…,vn}是頂點(diǎn)集, 表示構(gòu)建PPI 網(wǎng)絡(luò)的蛋白質(zhì)集合;
(2)E0={(v i,v j)|v i,v j∈V0}是頂點(diǎn)之間邊的集合, 表示蛋白質(zhì)之間相互作用的集合.
在無向圖中, 若兩個(gè)頂點(diǎn)是一條邊的兩個(gè)端點(diǎn), 則稱這兩個(gè)頂點(diǎn)相鄰.即對任何vi,v j∈V0, 若存在ei∈E0, 使得ei=(vi,vj), 則vi和vj為一組相鄰頂點(diǎn).若圖G0中任意兩個(gè)不同的頂點(diǎn)之間都是相鄰的, 則稱G0為完全圖.對于圖G0=(V0,E0)和圖G0′=(V0′,E0′), 若V0′ ?V0且E0′ ?E0, 則稱G0′是G0的子圖.設(shè)G0=(V0,E0),V1?V0且V1≠?, 由V1構(gòu)建邊集E1={(vi,v j)|vi,v j∈V1}, 則稱圖G[V1]=(V1,E1)為關(guān)于V1的導(dǎo)出子圖.
鄰接矩陣表示圖中頂點(diǎn)之間相鄰關(guān)系的二維矩陣.靜態(tài)PPI 網(wǎng)絡(luò)G0的鄰接矩陣M0定義為
活躍PPI 網(wǎng)絡(luò)主要是利用基因表達(dá)水平時(shí)間序列數(shù)據(jù), 從靜態(tài)PPI 網(wǎng)絡(luò)中, 刪除不能在任意同一時(shí)點(diǎn)保持活躍的兩個(gè)蛋白質(zhì)的相互作用, 從而構(gòu)建出活躍狀態(tài)下的PPI 網(wǎng)絡(luò).
通常采用閾值方法判定蛋白質(zhì)的活躍性, 即如果一個(gè)蛋白質(zhì)在某一時(shí)點(diǎn)的基因表達(dá)水平值超過閾值, 則此蛋白質(zhì)在該時(shí)點(diǎn)處于活躍狀態(tài).如果靜態(tài)PPI 網(wǎng)絡(luò)中的某條邊(蛋白質(zhì)相互作用)關(guān)聯(lián)的兩個(gè)頂點(diǎn)(蛋白質(zhì))不存在相同的活躍時(shí)點(diǎn), 則該條相互作用不能保留在活躍PPI 網(wǎng)絡(luò)中.
活躍PPI 網(wǎng)絡(luò)構(gòu)建算法[5]如下:
?Algorithm 1: filterthe data of DIP Input: S-PPI, gene expression(EXP).Step 1: for each gene i do
calculate its active threshold Active_th(i) for t=1 to T do if expit>Active_th(i) apit=1 end if end for end for ?Step 2: for each edge(u,v)∈E in S-PPI do for t=1 to T do if there does not exist a time point i which satisfies apvt*aput=1 remove edge(u,v) from E end if end for end for Step 3: Output edge set E
該算法通過基因表達(dá)水平數(shù)據(jù)集以及活躍閾值的設(shè)定, 對靜態(tài)PPI 網(wǎng)絡(luò)中各蛋白質(zhì)在整個(gè)周期各時(shí)點(diǎn)的活躍性進(jìn)行判斷; 隨后對靜態(tài)PPI 網(wǎng)絡(luò)中每條相互作用進(jìn)行活躍性判斷, 對每條相互作用的兩個(gè)蛋白質(zhì)頂點(diǎn)進(jìn)行共同活躍時(shí)點(diǎn)匹配, 若兩個(gè)蛋白質(zhì)頂點(diǎn)在整個(gè)周期中存在相同的活躍時(shí)點(diǎn), 則將其放入活躍相互作用邊集中.由活躍相互作用邊集構(gòu)成的PPI 網(wǎng)絡(luò), 稱為活躍PPI 網(wǎng)絡(luò).
活躍PPI 網(wǎng)絡(luò)根據(jù)基因表達(dá)時(shí)間序列數(shù)據(jù), 從靜態(tài)PPI 網(wǎng)絡(luò)中導(dǎo)出一個(gè)活躍PPI 網(wǎng)絡(luò).與活躍PPI 網(wǎng)絡(luò)不同, 多層PPI 網(wǎng)絡(luò)將構(gòu)建多個(gè)時(shí)點(diǎn)的活躍PPI 網(wǎng)絡(luò).它對每個(gè)時(shí)點(diǎn)的所有蛋白質(zhì)進(jìn)行活躍性判斷, 并從靜態(tài)PPI 網(wǎng)絡(luò)導(dǎo)出每個(gè)時(shí)點(diǎn)的活躍PPI 子網(wǎng)絡(luò), 從而構(gòu)成多個(gè)時(shí)點(diǎn)的活躍PPI 網(wǎng)絡(luò).
假設(shè)G0=(V0,E0)為靜態(tài)PPI 網(wǎng)絡(luò).用EXP表示蛋白質(zhì)的基因表達(dá)水平矩陣, 其中expit表示蛋白質(zhì)i在觀測時(shí)點(diǎn)t的基因表達(dá)水平值, 這里t=1,2,…,T.用AP表示蛋白質(zhì)的活躍性矩陣, 其中apit表示蛋白質(zhì)i在觀測時(shí)點(diǎn)t的活躍性, 其計(jì)算公式為
其中TH為活躍閾值.用at={vi|apit=1}表示在t時(shí)點(diǎn)的所有活躍蛋白質(zhì)的頂點(diǎn)集, 則G0[at]為根據(jù)活躍頂點(diǎn)集at從靜態(tài)PPI 網(wǎng)絡(luò)G0中導(dǎo)出的t時(shí)點(diǎn)的活躍PPI 子網(wǎng)絡(luò).所有時(shí)點(diǎn)的活躍PPI 子網(wǎng)絡(luò)構(gòu)成多層PPI網(wǎng)絡(luò)
多層PPI 網(wǎng)絡(luò)構(gòu)建算法如下:
Algorithm2: build Multilayer PPIs Input: S-PPI, gene expression(EXP).Step 1: for each gene i do calculate its active threshold Active_th(i) for t=1 to T do if expit>Active_th(i) apit=1 end if end for end for Step 2: for t=1 to T do for each gene i do if apit==1 do node i add intoat end if end for end for
Step 3: for t=1 to T do for each edge(v,u) in S-PPI do if v and u in at edge(v,u) add into PPIt end if end for end for Step 4: Output Multilayer PPIs ?
下面通過例子說明多層PPI 網(wǎng)絡(luò)的構(gòu)建過程.
(1) 靜態(tài)PPI 網(wǎng)絡(luò)構(gòu)建以及基因表達(dá)水平文件獲取.靜態(tài)PPI 網(wǎng)絡(luò)G0=(V0,E0), 其中V0={p1,p2,p3,p4},E0={(p1,p2),(p1,p4),(p3,p4)}, 如圖1所示.
圖1 蛋白質(zhì)活躍性
表1為靜態(tài)PPI 網(wǎng)絡(luò)G0中4 個(gè)蛋白質(zhì)在3 個(gè)時(shí)點(diǎn)的基因表達(dá)水平數(shù)EXP.
(2) 蛋白質(zhì)活躍時(shí)點(diǎn)判斷及活躍蛋白質(zhì)集合獲取.依據(jù)3Sigma閾值方法[5], 計(jì)算得到各蛋白質(zhì)的基因表達(dá)水平的活躍閾值分別為: 0.6146, 0.3785, 1.7823, 0.7258.然后得到蛋白質(zhì)活躍性矩陣AP, 見表2.
根據(jù)表2, 得到各時(shí)點(diǎn)活躍蛋白質(zhì)頂點(diǎn)集a1={p3,p4},a2={p1,p4},a3={p1,p2}.
表1 基因表達(dá)水平矩陣EXP
表2 蛋白質(zhì)活躍性矩陣AP
(3) 導(dǎo)出各時(shí)點(diǎn)的活躍時(shí)點(diǎn)網(wǎng)絡(luò).根據(jù)各時(shí)點(diǎn)的活躍蛋白質(zhì)集合, 從靜態(tài)PPI 網(wǎng)絡(luò)導(dǎo)出各時(shí)點(diǎn)的活躍時(shí)點(diǎn)網(wǎng)絡(luò), 如圖2所示.最后組合為多層PPI 網(wǎng)絡(luò)G={G0[a1],G0[a2],G0[a3]}.
圖2 多層PPI 網(wǎng)絡(luò)
蛋白質(zhì)相互作用采用釀酒酵母(Saccharomyces cerevisiae)蛋白質(zhì)相互作用數(shù)據(jù)集(DIP 數(shù)據(jù)庫[10]), 包括5093 個(gè)蛋白質(zhì)及24743 條相互作用.基因表達(dá)水平數(shù)據(jù)來自GES3431(geneExpressionOmnibus[11]), 該文件包含了共6777 個(gè)蛋白質(zhì)36 個(gè)時(shí)點(diǎn)(3 個(gè)周期)的基因表達(dá)水平數(shù)據(jù).兩個(gè)數(shù)據(jù)集具有4846 個(gè)相同的蛋白質(zhì)(占DIP 數(shù)據(jù)集的95%).因此, 用于構(gòu)建活躍PPI 網(wǎng)絡(luò)和多層PPI 網(wǎng)絡(luò)的蛋白質(zhì)為4846 個(gè).在4846個(gè)蛋白質(zhì)中, 已知的關(guān)鍵蛋白質(zhì)有1167 個(gè)[5].
中心性度量方法是關(guān)鍵蛋白質(zhì)識別的一種常用方法[6~8].表3列舉了三種常用的中心性方法:DC(度中心性)、NC(邊聚類系數(shù)中心性)和LAC(局部平均連接中心性), 其中NC和LAC是網(wǎng)絡(luò)拓?fù)浞椒ㄖ袃煞N識別效果最好的中心性方法.
表3 中心性方法
表3中,deg(v)表示蛋白質(zhì)v的度;其中Zv,u表示包含節(jié)點(diǎn)v,u的這條邊的三角形的個(gè)數(shù);Nv是節(jié)點(diǎn)v的鄰居集合,Cv是靜態(tài)PPI 網(wǎng)絡(luò)中只包含節(jié)點(diǎn)v鄰居節(jié)點(diǎn)的導(dǎo)出子圖,deg Cv(w)為相應(yīng)導(dǎo)出子圖中節(jié)點(diǎn)w的度.
下面以DC中心性方法為例, 介紹加權(quán)多層中心性方法的計(jì)算步驟.假設(shè)多層PPI 網(wǎng)絡(luò)已構(gòu)建完成.
(1) 對每個(gè)時(shí)點(diǎn)子網(wǎng)絡(luò), 計(jì)算該網(wǎng)絡(luò)中各活躍蛋白質(zhì)的中心性度量值.用DC(i,t)表示在時(shí)點(diǎn)t子網(wǎng)絡(luò)中蛋白質(zhì)vi的DC中心性度量值.
(2) 對每個(gè)蛋白質(zhì)在各個(gè)子網(wǎng)絡(luò)中的中心性值進(jìn)行加權(quán)求和.蛋白質(zhì)vi的DC中心性值為其中wt是第t時(shí)點(diǎn)子網(wǎng)絡(luò)的權(quán)值系數(shù).
(3) 對所有蛋白質(zhì)的DC中心性值進(jìn)行排序, 取中心性值排名前100~600 的蛋白質(zhì), 計(jì)算其中包含已知關(guān)鍵蛋白質(zhì)的數(shù)量, 得到DC中心性方法的關(guān)鍵蛋白質(zhì)識別率.
NC和LAC均采用類似方法進(jìn)行計(jì)算.
用S-PPI、A-PPI、M-PPI 分別表示靜態(tài)PPI 網(wǎng)絡(luò)、活躍PPI 網(wǎng)絡(luò)和多層PPI 網(wǎng)絡(luò).下面通過實(shí)驗(yàn)比較三種PPI 網(wǎng)絡(luò)中關(guān)鍵蛋白質(zhì)的識別率.其中對A-PPI 和M-PPI, 蛋白質(zhì)活躍性判定均采用3Sigma 方法[5].
給定M-PPI 各時(shí)點(diǎn)層加權(quán)系數(shù)W={62, 103, 39, 3, 25, 3, 12, 8, 718, 49, 55, 9}.分別針對三種中心性DC、NC和LAC, 采用加權(quán)多層中心性方法, 計(jì)算閾值k從 1- 到4 變化范圍內(nèi)的關(guān)鍵蛋白質(zhì)識別數(shù)量.實(shí)驗(yàn)結(jié)果如圖3所示.
圖3 M-PPI 中三種中心性方法的關(guān)鍵蛋白質(zhì)識別情況
由圖3可看出, 針對中心性DC, 其識別效果隨著閾值系數(shù)k的增大而提高, 但當(dāng)閾值系數(shù)增長至3左右后, 隨著閾值系數(shù)的增大, 識別效果開始隨之降低; 對于中心性算法NC和LAC, 隨著閾值系數(shù)的提高, 其識別效果波動(dòng)較小, 變化幅度不大, 當(dāng)閾值系數(shù)增長至1.8 左右后, 其識別效果開始有了較為明顯的下滑趨勢, 當(dāng)閾值系數(shù)增長至3.6 左右時(shí), 識別效果大幅降低.
表4列出了三種中心性方法在M-PPI 中識別結(jié)果優(yōu)于A-PPI 最佳識別結(jié)果[5]的k值范圍.
表4 閾值k 參考范圍
由表4可看出, 對于DC, 當(dāng)k值為0.9~3.4 時(shí), 其M-PPI 的關(guān)鍵蛋白質(zhì)識別數(shù)量優(yōu)于A-PPI.對于NC, 當(dāng)k值為1.4~3.2 時(shí), 其M-PPI 的關(guān)鍵蛋白質(zhì)識別數(shù)量優(yōu)于A-PPI.對于LAC, 當(dāng)k值為1.2~1.9 和2.7~3.2時(shí), 其M-PPI 的關(guān)鍵蛋白質(zhì)識別數(shù)量優(yōu)于A-PPI.綜上可知, 當(dāng)k值為1.4~1.9 和2.7~3.2 時(shí), M-PPI 中三種中心性方法均可取得優(yōu)于A-PPI 的識別結(jié)果.圖4和表5列出了S-PPI、A-PPI(2.5)k=和M-PPI(1.9k=,W={62, 103, 39, 3, 25, 3, 12, 8, 718, 49, 55, 9})中DC、NC、LAC的關(guān)鍵蛋白質(zhì)識別結(jié)果.
圖4 關(guān)鍵蛋白質(zhì)識別對比圖
表5 關(guān)鍵蛋白質(zhì)識別數(shù)量
對于DC, 相比S-PPI、A-PPI, M-PPI 在Top600 中的識別率分別提升了29.8%和8.3%.對于NC, M-PPI相比S-PPI 和A-PPI, 分別提升了10.4%和5.2%.對于LAC, M-PPI 比S-PPI 和A-PPI 分別提升了10.1%和4.7%.表6給出了W和k取不同值時(shí)三種中心性方法在Top100~600 的識別結(jié)果, 均優(yōu)于目前A-PPI 中三種中心性的最優(yōu)識別結(jié)果.
表6 各加權(quán)系數(shù)下最優(yōu)閾值k 及對應(yīng)關(guān)鍵蛋白質(zhì)識別數(shù)量
由表6可以看出, 針對不同的加權(quán)系數(shù), 其對應(yīng)M-PPI 網(wǎng)絡(luò)的最優(yōu)閾值存在一定差異; 但針對不同加權(quán)系數(shù), 其閾值系數(shù)為2 或3 左右時(shí), M-PPI 網(wǎng)絡(luò)的識別結(jié)果均能優(yōu)A-PPI 網(wǎng)絡(luò)的識別結(jié)果.
本文在靜態(tài)PPI 網(wǎng)絡(luò)的基礎(chǔ)上, 利用基因表達(dá)水平時(shí)間序列數(shù)據(jù), 構(gòu)建了一種蛋白質(zhì)相互作用多層網(wǎng)絡(luò), 提出了基于多層網(wǎng)絡(luò)的中心性加權(quán)方法, 提高了關(guān)鍵蛋白質(zhì)的識別率.