管澤雨，邱嘉迪，劉文碩，趙蘊杰

(華中師范大學(xué)物理科學(xué)與技術(shù)學(xué)院 生物物理研究所，武漢 430079)

蛋白質(zhì)由20種不同的氨基酸經(jīng)肽鍵聚合組成，通過形成特定的空間三級結(jié)構(gòu)以實現(xiàn)催化和調(diào)控等不同的生物學(xué)功能.不同的氨基酸排列順序會形成不同的蛋白質(zhì)結(jié)構(gòu)，蛋白質(zhì)三維空間結(jié)構(gòu)決定了其生物學(xué)功能[1-2].因此，蛋白質(zhì)三維結(jié)構(gòu)研究對理解其如何發(fā)揮生物學(xué)功能和設(shè)計相關(guān)的藥物具有重要的意義.近年來，蛋白質(zhì)殘基相互作用網(wǎng)絡(luò)普遍應(yīng)用于蛋白質(zhì)相關(guān)問題的研究.該方法中網(wǎng)絡(luò)的節(jié)點為組成蛋白質(zhì)的殘基，網(wǎng)絡(luò)的邊為非共價鍵殘基相互作用(如范德瓦爾斯和靜電相互作用等)[3].基于蛋白質(zhì)殘基相互作用網(wǎng)絡(luò)，可以進一步利用圖論的方法研究蛋白質(zhì)結(jié)構(gòu)穩(wěn)定性[4-5]，蛋白質(zhì)動力學(xué)[6-8]，酶活性和變構(gòu)調(diào)節(jié)[9]，信號轉(zhuǎn)導(dǎo)[10-11]等問題，為解決這些問題提供了一個嶄新的視角.例如，Vendruscolo等人通過蛋白質(zhì)殘基相互作用網(wǎng)絡(luò)的聚類系數(shù)(clustering coefficient)，平均路徑長度(average shortest path length)和中介中心度(betweenness centrality)的分析，提出了中介中心度極大值殘基是折疊過程的關(guān)鍵氨基酸[8].Amitai等人通過計算蛋白質(zhì)殘基相互作用網(wǎng)絡(luò)整體的接近中心度(closeness centrality)和相對溶劑可及性(relative solvent accessibility)，可以較為有效地識別蛋白質(zhì)活性位點氨基酸，在包含178個典型蛋白質(zhì)的數(shù)據(jù)集中正確率達到70%[9].

目前已有一些分析蛋白質(zhì)殘基相互作用網(wǎng)絡(luò)的工具.例如，AMINONET可計算疏水、親水或帶電氨基酸組成網(wǎng)絡(luò)的拓撲屬性[12]；RING可構(gòu)建殘基相互作用網(wǎng)絡(luò)，并通過Cytoscape[13]計算網(wǎng)絡(luò)的拓撲性質(zhì)[14-15]；NAPS[16]可分析蛋白質(zhì)網(wǎng)絡(luò)邊和節(jié)點的各種中心度性質(zhì)，查找兩個節(jié)點的最短路徑，k連通子網(wǎng)等.然而，現(xiàn)有的蛋白質(zhì)殘基相互作用網(wǎng)絡(luò)分析工具易用性較差，且網(wǎng)站訪問較不穩(wěn)定.例如，AMINONET需要在Java環(huán)境下運行，RING則需要將結(jié)果導(dǎo)出到Cytoscape進行分析，而Cytoscape本身也需要在Java環(huán)境下運行，需要依賴其它軟件和環(huán)境，使用不方便.另一方面，隨著人類基因組計劃的完成，蛋白組學(xué)迅猛發(fā)展，大量的疾病相關(guān)基因被發(fā)現(xiàn)，藥物作用的靶標(biāo)分子急劇增加，現(xiàn)有方法缺乏對蛋白質(zhì)結(jié)合口袋的分析，極大限制了殘基相互作用網(wǎng)絡(luò)方法的應(yīng)用范圍[17].因此，亟需易用性較強并可分析結(jié)合口袋的蛋白質(zhì)殘基相互作用網(wǎng)絡(luò)模型.

本文建立了基于D3[18]和NGL Viewer[19]的蛋白質(zhì)分析平臺.用戶僅需提交蛋白質(zhì)的PDB結(jié)構(gòu)信息，既可快速搭建蛋白質(zhì)殘基相互作用網(wǎng)絡(luò)，計算網(wǎng)絡(luò)的拓撲性質(zhì)，實現(xiàn)網(wǎng)絡(luò)的可視化，并結(jié)合DogSiteScorer計算的結(jié)合口袋信息計算它們的拓撲性質(zhì)，對理解蛋白質(zhì)的結(jié)構(gòu)、生物學(xué)功能和相關(guān)藥物設(shè)計有重要的意義.

1 算法

1.1 蛋白質(zhì)殘基相互作用網(wǎng)絡(luò)

本文用蛋白質(zhì)的三維結(jié)構(gòu)信息構(gòu)建蛋白質(zhì)殘基相互作用網(wǎng)絡(luò).網(wǎng)絡(luò)由節(jié)點和邊兩部分組成：蛋白質(zhì)結(jié)構(gòu)中的氨基酸為網(wǎng)絡(luò)中的節(jié)點.蛋白質(zhì)結(jié)構(gòu)統(tǒng)計分析表明，不相鄰的氨基酸之間主要有兩個距離分布的峰值.第一個距離峰值為0.3～0.5 nm，主要是氫鍵等短程相互作用；第二個距離峰值為0.7～0.8 nm，主要為靜電等長程相互作用[20].因此，蛋白質(zhì)相互作用研究大多數(shù)以0.8 nm作為距離截斷，如果兩氨基酸間任一對原子的距離小于0.8 nm則定義該氨基酸—氨基酸形成網(wǎng)絡(luò)的邊.

1.2 網(wǎng)絡(luò)特征計算

本文提供了3類網(wǎng)絡(luò)特征的計算，分別為度中心度，接近中心度和中介中心度，具體的計算公式如表1所示.度中心度直觀上反映了一個節(jié)點在網(wǎng)絡(luò)中的重要程度，定義為與該節(jié)點連接的邊的數(shù)目；接近中心度描述了網(wǎng)絡(luò)中的節(jié)點通過網(wǎng)絡(luò)到達其它節(jié)點的難易程度，數(shù)值為該節(jié)點到其它所有能到達節(jié)點的平均距離的倒數(shù)，反映了該節(jié)點對其它節(jié)點施加影響的能力；中介中心度描述了節(jié)點在整個網(wǎng)絡(luò)中的作用和影響力，定義為網(wǎng)絡(luò)中所有最短路徑經(jīng)過該節(jié)點的數(shù)量[21].

表1 3種centrality的定義Tab.1 The definitions of 3 kinds of centrality

1.3 網(wǎng)站服務(wù)的基本框架

蛋白質(zhì)殘基相互作用網(wǎng)絡(luò)服務(wù)用到了HTML、Javascript、PHP、MATLAB等語言以及D3.js、NGL.js等Javascript庫，基本框架如圖1所示.

主要流程為：

1) 服務(wù)器端調(diào)用FileRead.php，將臨時文件傳入到服務(wù)器的Inputs文件夾，對文件解壓.

2) 調(diào)用Matlab編寫的可執(zhí)行程序，讀取pdb文件中的數(shù)據(jù)，生成每個殘基的3種centrality數(shù)據(jù)，具體為(流程圖見圖1(b))：①將pdb文件中的蛋白質(zhì)三維結(jié)構(gòu)數(shù)據(jù)轉(zhuǎn)換為殘基網(wǎng)絡(luò)的鄰接矩陣，若其中兩個殘基不相鄰且存在一對原子距離小于0.8 nm則定義為相鄰.②按照1.2中centrality的定義計算出每個殘基的3種centrality值.

3) 清空Inputs文件夾，以便再次傳入數(shù)據(jù).將得到的每個殘基的centrality值及殘基之間的連接信息寫入Outputs文件夾中的centrality.json文件,生成3種centrality的直方圖，折線圖，散點圖.

4) 瀏覽器加載centrality.json文件，默認以closeness的倍數(shù)為節(jié)點的半徑繪制所上傳蛋白質(zhì)的力導(dǎo)向圖，并調(diào)用NGL Viewer顯示蛋白質(zhì)三維結(jié)構(gòu)，顯示3種centrality的直方圖，折線圖，散點圖.

1.4 D3和NGL Viewer的實現(xiàn)

D3是一個基于web標(biāo)準的JavaScript可視化庫，D3可以借助SVG,Canvas以及HTML將數(shù)據(jù)生動的展現(xiàn)出來[18].使用基于D3的力導(dǎo)向算法(Forced-Directed Algorithm)的相關(guān)API來實現(xiàn)可視化網(wǎng)絡(luò).力導(dǎo)向算法是Eades于1984年提出的一種布點作圖算法[22]，其基本思想為將網(wǎng)絡(luò)看成一個頂點為鋼環(huán)、邊為彈簧的物理系統(tǒng)，系統(tǒng)被賦予某個初始狀態(tài)后，彈簧彈力的作用使鋼環(huán)運動，直到系統(tǒng)總能量達到最小值時停止.通過改進力學(xué)模型可以得到不同的算法，并實現(xiàn)以下3點：1)節(jié)點分布均勻；2)邊交叉最小化；3)具有對稱性[23].此外，還通過添加tick事件來不斷更新圖形系統(tǒng)，實現(xiàn)動態(tài)的推拽效果.

NGL Viewer是實現(xiàn)蛋白質(zhì)結(jié)構(gòu)可視化的Web應(yīng)用程序，用戶可上傳并顯示蛋白質(zhì)三維結(jié)構(gòu)，不受第三方插件(如Flash和Java小程序)的影響.NGL Viewer支持常見的蛋白質(zhì)結(jié)構(gòu)文件格式(如pdb、mmcif)和各種分子表示(例如cartoon、spacefill、licorice、surface).此外，NGL Viewer可嵌入到網(wǎng)站中實現(xiàn)可視化.

2 結(jié)果分析

2.1 蛋白質(zhì)殘基相互作用網(wǎng)絡(luò)在線服務(wù)的主要功能

研究表明，蛋白質(zhì)殘基相互作用網(wǎng)絡(luò)呈現(xiàn)小世界網(wǎng)絡(luò)特征[8,24]，力導(dǎo)向圖在小世界網(wǎng)絡(luò)中有良好的可視化效果[25-26].力導(dǎo)向圖易于理解和實現(xiàn)，可以畫出相當(dāng)優(yōu)美的圖形布局，充分展現(xiàn)出圖的整體結(jié)構(gòu)及其自同構(gòu)特征[26].因此，基于力導(dǎo)向的可視化算法可以較好的描繪蛋白質(zhì)殘基相互作用網(wǎng)絡(luò).圖2為細胞周期蛋白依賴激酶2(CDK2，PDB code：1fin，A鏈)前100個氨基酸的蛋白質(zhì)殘基相互作用網(wǎng)絡(luò)，其中不同的5種顏色分別對應(yīng)5個不同的口袋，黑色為默認顏色，表示不參與口袋的形成，右側(cè)顯示為該蛋白質(zhì)的三維結(jié)構(gòu).

圖2 細胞周期蛋白依賴激酶A鏈前100個氨基酸的力導(dǎo)向圖Fig.2 Force-directed graph of top 100 amino acid of CDK2

網(wǎng)站的主要功能有以下幾個方面.

1) 當(dāng)光標(biāo)浮在節(jié)點上方時，顯示該節(jié)點的編號，closeness centrality，betweenness centrality，degree centrality，所在口袋編號(可能同時屬于多個口袋)；

2) 按照各centrality值高低搜索對應(yīng)殘基并染色；

3) 搜索某個口袋的所有殘基并染色，觀察其在網(wǎng)絡(luò)中的分布；

4) 列出各個口袋的各平均centrality值；

5) 利用NGL Viewer顯示蛋白質(zhì)三維結(jié)構(gòu)，并根據(jù)計算的centrality分布將三維結(jié)構(gòu)染色；

6) 由D3繪制的力導(dǎo)向圖具有拖拽功能，可以動態(tài)顯示蛋白質(zhì)殘基相互作用網(wǎng)絡(luò)；

7) 作出各centrality的區(qū)間分布圖，centrality隨氨基酸索引的折線圖，不同centrality分別作橫縱坐標(biāo)的散點圖.

2.2 蛋白質(zhì)殘基相互作用網(wǎng)絡(luò)在線服務(wù)的使用方法

在線服務(wù)分為任務(wù)區(qū)，演示區(qū)和教程區(qū)，使用的主要步驟如下：

1) 上傳蛋白質(zhì)PDB結(jié)構(gòu)文件，DoGSiteScorer(https://proteins.plus/#dogsite)計算得到的口袋結(jié)構(gòu)信息(壓縮為zip格式)，提交任務(wù).

2) 計算結(jié)果如圖3所示，主要分為兩個部分，區(qū)域I顯示蛋白質(zhì)的力導(dǎo)向圖，區(qū)域II顯示相應(yīng)的三維結(jié)構(gòu)圖，初始時默認染色第一個口袋.此時，有兩種染色功能可供選擇：①根據(jù)centrality值的大小染色，centrality有closeness centrality、betweenness centrality、degree centrality三種，可在區(qū)域I下拉列表中選擇.在文本框中輸入整數(shù)m(下標(biāo)從0開始)與n，點擊color后突出顯示排名在m與n之間的殘基，將鼠標(biāo)懸停在某個節(jié)點的上方，可以查看該節(jié)點對應(yīng)殘基的序號，各centrality值以及該氨基酸所在的口袋編號.②根據(jù)口袋染色，數(shù)據(jù)加載完成后，會自動生成與上傳口袋數(shù)相同的復(fù)選框按鈕，選擇要染色的口袋，點擊color，選中的口袋的所有氨基酸將依口袋的不同被染成不同的顏色.當(dāng)上傳蛋白質(zhì)的殘基數(shù)目較多時(例如大于300個)，可滾動鼠標(biāo)實現(xiàn)圖的縮放，縮放比例范圍為[0.5,2].區(qū)域II用于顯示蛋白質(zhì)三維結(jié)構(gòu)(提供Cartoon、Spacefill、Licorice、Surface 4種顯示方式)并根據(jù)計算的中心度值將其染色，將中心度值按照數(shù)值從大到小的順序進行三等分，分別為高中心度值(顯示為紅色)、中中心度值(顯示為綠色)以及低中心度值(顯示為藍色)，用戶可以在下拉列表中選擇不同類型的中心度進行染色.在區(qū)域I與區(qū)域II的下方，分別顯示3種中心度的區(qū)間分布圖(橫坐標(biāo)為相應(yīng)中心度的區(qū)間取值范圍，縱坐標(biāo)為氨基酸中心度取值在該范圍內(nèi)的數(shù)目)、中心度值隨殘基序號的折線圖，不同中心度分別為橫縱坐標(biāo)的散點圖.

3) 點擊Download區(qū)域的json、jpeg文件可以打開或者下載相應(yīng)文件.

利用力導(dǎo)向圖可以更加深入和直觀的了解蛋白質(zhì)結(jié)構(gòu)，挖掘有用的信息.

圖3 頁面主體功能區(qū)Fig.3 Main body of the resulting page

2.3 蛋白質(zhì)殘基相互作用網(wǎng)絡(luò)的網(wǎng)絡(luò)特征計算

以計算p38alpha c162s突變體(PDB code：1r3c)的closeness為例，上傳蛋白質(zhì)PDB結(jié)構(gòu)文件及DogSiteScorer網(wǎng)站計算的口袋結(jié)構(gòu)信息并提交任務(wù)[27].在結(jié)果頁面，選擇圖3所示區(qū)域I底部的ByCentrality染色方案，右側(cè)下拉列表中選擇closeness選項，在輸入框中輸入0，10(表示將closeness排名前10的氨基酸染成紅色)，點擊Color，將光標(biāo)移到紅色節(jié)點上方顯示出殘基名稱、序號、centrality值及所在口袋編號，所得結(jié)果如表2.其中P_0為ATP口袋，因而一半的殘基在ATP口袋，說明了ATP口袋有較強的保守性.betweenness的計算結(jié)果如表3所示.

可以看出closeness和degree近似服從正態(tài)分布，betweenness近似服從長尾分布.進一步，可下載centrality.json文件，找出betweenness最大的前10位，進行細致的分析(如表3所示).圖5為p38alpha c162s突變體結(jié)構(gòu)(PDB code：1r3c)的cartoon圖，紅色部分為betweenness值排名前10的殘基，主要分布在ATP口袋周圍，說明ATP口袋中的殘基對蛋白質(zhì)網(wǎng)絡(luò)拓撲結(jié)構(gòu)有較大的影響.值可預(yù)測蛋白質(zhì)的結(jié)合位點氨基酸.

表2 1r3c中closeness值Top 10對應(yīng)的氨基酸Tab.2 Top 10 amino acid by closeness of 1r3c

表3 1r3c中betweenness值Top 10對應(yīng)的氨基酸Tab.3 Top 10 amino acid by betweenness of 1r3c

網(wǎng)絡(luò)特征closeness centrality，betweenness centrality和degree centrality的分布規(guī)律如圖4所示.

圖4 centrality的分布規(guī)律Fig.4 The distribution of 3 kinds of centrality

圖5 p38alpha c162s突變體的結(jié)構(gòu)(PDB code：1r3c)Fig.5 The cartoon presentation of 1r3c

2.4 利用closeness識別關(guān)鍵氨基酸

研究表明，蛋白質(zhì)的活性位點，配體結(jié)合位點，進化保守殘基的closeness值大多較高[28].因此，根據(jù)蛋白質(zhì)殘基相互作用網(wǎng)絡(luò)中節(jié)點的closeness以枯草桿菌蛋白酶為例(PDB code:1bh6)，圖6為枯草桿菌蛋白酶三維結(jié)構(gòu)中網(wǎng)絡(luò)拓撲性質(zhì)的closeness分布.將closeness值按照數(shù)值從大到小的順序進行三等分，紅色表示closeness值較高的殘基，綠色表示closeness值中等的殘基，藍色表示closeness值較低的殘基.結(jié)果表明，枯草桿菌蛋白酶表面closeness值較高的殘基大多數(shù)分布在小分子的結(jié)合位點區(qū)域.

圖6 枯草桿菌蛋白酶(1bh6)closeness分布圖Fig.6 The distribution of closeness of subtilisin(PDB code:1bh6)

2.5 基于betweenness識別蛋白質(zhì)折疊過渡態(tài)系綜中的關(guān)鍵氨基酸[8]

研究表明，蛋白質(zhì)殘基相互作用網(wǎng)絡(luò)不同于隨機網(wǎng)絡(luò)與規(guī)則網(wǎng)絡(luò)，具有處于兩者之間的小世界網(wǎng)絡(luò)的特性，即具有較大的聚類系數(shù)和較小的平均最短路徑[5].因此，利用網(wǎng)絡(luò)節(jié)點的betweenness可以有效識別蛋白質(zhì)折疊過程中過渡態(tài)系綜的關(guān)鍵殘基.Vendruscolo等[29]通過實驗確定了酰磷酸酶(PDB code：1aps)折疊過程中的關(guān)鍵殘基，具體為TYR11、PRO54、PHE94.圖7為酰磷酸酶殘基相互作用網(wǎng)絡(luò)的betweenness數(shù)值分布，TYR11和PHE94的betweenness較高.

圖7 1aps蛋白的betweenness值隨殘基序列的分布Fig.7 The distribution of betweenness of protein 1aps(pdb code：1aps )

3 總結(jié)與展望

蛋白質(zhì)殘基相互作用網(wǎng)絡(luò)模型為研究蛋白質(zhì)的結(jié)構(gòu)與功能的關(guān)系提供了嶄新的視角.網(wǎng)絡(luò)拓撲結(jié)構(gòu)中的接近中心度、中介中心度等性質(zhì)，反映了單個氨基酸與蛋白質(zhì)整體結(jié)構(gòu)的關(guān)系[9].該方法在蛋白質(zhì)穩(wěn)定性，蛋白質(zhì)折疊動力學(xué)，控制酶活性和變構(gòu)調(diào)節(jié)，信號轉(zhuǎn)導(dǎo)等方面取得了一系列成果.

本文利用D3.js，HTML，PHP，Javascript，MATLAB和NGL Viewer等程序模塊搭建了蛋白質(zhì)殘基相互作用網(wǎng)絡(luò)在線服務(wù)及可視化分析平臺.用戶可以利用蛋白質(zhì)結(jié)構(gòu)信息搭建網(wǎng)絡(luò)結(jié)構(gòu)，并計算接近中心度等網(wǎng)絡(luò)的拓撲結(jié)構(gòu)信息.長度約300個氨基酸的中等大小蛋白質(zhì)，計算速度約為30 s.另外，該分析平臺還可分析蛋白質(zhì)的結(jié)合口袋特征，有較強的可擴展性，可添加網(wǎng)絡(luò)最短路徑，k個節(jié)點的極大完全子圖(k-clique)等模塊單元.結(jié)果表明，蛋白質(zhì)殘基相互作用網(wǎng)絡(luò)在線服務(wù)及可視化分析平臺對理解蛋白質(zhì)結(jié)構(gòu)，調(diào)控功能和藥物設(shè)計的相關(guān)研究有重要的幫助.