距離和疏水模型輔助的蛋白質結構預測方法

王小奇,周曉根,胡 俊,張貴軍

1(浙江工業大學 信息工程學院,杭州 310023)2(密歇根大學計算醫學與生物信息學系，密歇根州安娜堡 45108，美國)

1 引 言

蛋白質是生命活動的物質基礎,對人類的健康起著決定性作用.準確預測蛋白質的結構對疾病研究、藥物設計等方面都有重要的科學意義[1].目前測定蛋白質結構的實驗方法主要有X射線晶體衍射、多維核磁共振和冷凍電鏡.但是這些實驗測定的方法代價昂貴[2-4],同時對于藥物設計的主要靶標—膜蛋白而言,實驗方法極難測定其結構.因此,根據Anfinsen法則[5],從氨基酸序列出發,利用計算機技術設計適當的算法從頭預測蛋白質結構已成為生物信息學研究中的熱點[6].

從頭預測方法近十年來迅速發展,取得了很大的進步,尤其是華盛頓大學Baker研究團隊開發的Rosetta[7,8]和密歇根大學Zhang研究小組開發的QUARK[9,10]在最近幾屆CASP[11-13]比賽中表現突出,但針對一些序列較長的蛋白(序列長度大于100),現有的方法預測精度依然不高[11,14].這主要是因為現有方法的構象空間采樣能力不足,搜索不到近天然態的構象,并且能量函數的不精確性導致每個構象難以量化,增加了該問題的復雜性,進而影響了從頭預測精度.

針對構象優化問題,國內外研究者相繼提出了大量的構象空間優化方法,如進化算法[15-18]、蒙特卡羅[19,20]、構象空間退火[21-24]、副本交換[25]等被廣泛應用于蛋白質結構預測中.研究表明,蛋白質結構先驗知識可以降低能量函數不精確帶來的誤差,并有效地減小構象的搜索空間,進而提高蛋白結構的預測精度[26-29].APL[30]利用二面角信息來提高蛋白質結構的預測精度;RMA[14]和AIMOES[31]分別利用預測的二級結構來增強構象空間的采樣;EVfold[32]和PconsFold[33]從預測的蛋白質接觸圖(contact map)中提取結構信息,并應用于蛋白質結構預測中;CONFOLD[34]和SCDE[35]利用二級結構和接觸圖來預測蛋白結構.DeepMind開發的AlphaFold系統[36]在CASP13的表現非常突出,主要是因為通過神經網絡預測氨基酸之間的距離和化學鍵角度,然后再根據兩種物理屬性對結構進行評分,最后通過梯度下降優化,這也進一步證實了先驗的結構知識對預測蛋白質結構的有效性.

盡管接觸圖可以提高蛋白質結構的預測精度,但是傳統的接觸圖僅僅估計了每對殘基間的距離是否大于某一閾值(8?).距離譜[9,10]可以反應殘基間的距離分布信息,基于距離譜作者前期在進化算法的框架下[37-39],提出了一種蛋白質構象空間優化方法[40]表明距離譜可以有效地提高預測的精度.然而,在蛋白質的折疊過程中,疏水作用是主要作用之一,對蛋白質結構的形成與維持有著特殊的作用[41],也被廣泛地用來預測蛋白質結構,但現階段的研究中氨基酸的疏水性主要針對于簡化的HP模型[42,43].

針對現有預測方法的不足,在進化算法的框架下提出一種距離和疏水模型輔助的蛋白質結構預測方法(DHMA).首先,根據親疏水氨基酸的回轉半徑來指導構象采樣;然后,利用距離譜構建距離分布模型和疏水概率模型來指導種群的更新,進而提高預測的精度.在10個測試蛋白的實驗結果表明,DHMA是一種有效的蛋白質結構預測方法.

2 DHMA算法設計

在DHMA中,根據氨基酸的親疏水性與距離譜,分別構建親疏水性氨基酸的回轉半徑、距離分布模型和疏水概率模型來指導構象的優化.

2.1 親疏水氨基酸的回轉半徑差

在變異操作中,利用親疏水氨基酸的回轉半徑差來引導目標個體的變異.

疏水氨基酸的回轉半徑定義如公式(1):

(1)

同理,計算得到親水性氨基酸的旋轉半徑RGP,并且根據如下公式計算親疏水氨基酸的旋轉半徑差:

DRG=RGP-RGH

(2)

在穩定的蛋白質結構中,疏水氨基酸殘隱藏在蛋白質分子內部.因此在變異操作中,首先判斷片段組裝后的新個體是否滿足DRG>0,如果DRG>0,則再采用能量函數進行評估;否則直接拒絕此次片段插入.

2.2 距離分布模型

在選擇過程中,首先利用基于距離譜的打分模型來評價構象的質量.距離分布打分模型計算如式(3):

(3)

其中U是距離譜中殘基對的集合,p和q是U中的任一殘基對,dpq為構象中第p個殘基與第q個殘基的Cα原子間的歐氏距離,Npq(dpq)是在殘基對(p,q)的距離譜中dpq對應的片段數.

距離分布估計模型的分值在很多程度上可以反應構象的質量,即Ds值越大表明構象的結構越好.因此,在選擇過程中,首先利用距離分布估計模型來指導種群的更新.

2.3 疏水概率模型

在穩定的蛋白質結構中,疏水氨基酸隱藏在蛋白質分子內部而彼此相互接觸.在距離譜中,只統計了歐式距離小于9?的殘基對間的距離分布信息,因此利用距離譜中疏水殘基對的疏水性設計了打分模型,具體定義如式(4):

(4)

其中f(p,q)計算如下:

(5)

其中cp和cq分別是距離譜中第p和q個殘基的疏水值[44].

疏水概率模型pH根據如公式(6)計算:

(6)

其中Hs與Hs′分別代表目標個體與變異個體的疏水得分,其中T根據如下公式計算:

(7)

其中u(p,q)定義如下:

(8)

同理可以計算得到變異個體的T′值.

在選擇操作中,當變異個體的距離分布值比目標個體低時,利用疏水概率模型來輔助指導種群的更新,進一步提高DHMA的預測精度.

2.4 算法描述

DHMA算法的流程描述如下:

Step 1.初始化:進行參數初始化和種群初始化,具體過程如下:

Step1.1設置種群規模NP、溫度常數kT、拒絕升溫常數m、片段長度frag、迭代常數λ∈(0,1)和最大進化代數Gmax;利用Rosetta的第一和第二階段產生NP個個體作為初始種群Pinit={x1,…,xNP};

Step1.2設置目標個體的索引i=1,i∈[1,NP],進化代數g=1;

Step 2.變異操作:利用親疏水氨基酸的回轉半徑差和能量函數指導目標個體xi的變異,具體過程如下:

Step2.1利用片段組裝技術對目標個體xi進行片段插入,并根據公式(2)計算組裝后個體的DRG;若DRG≤0,則拒絕此次片段插入,否則計算組裝前后個體的能量變化,并根據蒙特卡羅(Monte Carlo)機制判斷是否接收此次片段的插入;

鑒于東海2號機組于1978年投運，原子能電力已于2017年11月向規制委提交延壽申請。要想讓這份申請獲得批準，該公司必須在2018年11月底這一法定截止日之前獲得規制委對下述兩份計劃的批準：安全強化措施詳細執行計劃和延壽改造詳細計劃。

Step3.2若Ds>Ds′,則目標個體xi被保留,否則根據公式(6)計算疏水概率pΗ;

Step 4.判斷是否結束:具體步驟如下:

Step4.1若i=NP,則設i=1,否則i=i+1,并且轉至步驟Step 2;

Step4.2若g>λGmax,則frag=3,否則frag=9;

Step4.3若g≥Gmax,則結束迭代,否則g=g+1,且轉至Step 2.

3 實驗和結果

3.1 實驗材料

在實驗中,從PDB庫中選取序列長度從64到125的10個測試蛋白,這些測試蛋白的折疊類型包括α、β和α/β,測試蛋白的具體信息如表1中的第2-4列所示.選取這些蛋白質做測試是因為它們在生物體內有著重要的作用[45,46],并且在現有文獻中這些蛋白被廣泛地研究[14,47，48],因此,這些蛋白可以較好地測試DHMA算法的性能.

測試蛋白的片段庫和距離譜分別通過ROBETTA[49]和QUARK服務器獲得,并且在構建片段庫和距離譜的過程中均去除同源模板.實驗的參數設置如下:種群規模NP=100、溫度常數kT=2、拒絕升溫常數m=150、迭代常數λ=0.6、最大迭代次數Gmax=1.5L,其中L為序列長度.在初始化過程中能量函數采用Rosetta的能量函數score0和score1,在變異操作中采用Rosetta的能量函數score3[7，8].通過計算預測模型與天然態結構中Cα原子的均方根偏差(RMSD)和TM-score[50]來評價算法的預測精度.本次實驗中,獨立運行DHMA算法10次,保留中間采樣的構象,進而產生一個構象池,在同樣的代價下,獨立運行1.5L條蒙特卡羅軌跡并保留第三和第四階段采樣的所有構象.

3.2 結果分析

實驗中DHMA的預測結果與Rosetta和QUARK預測的結果進行比較,不僅因為Rosetta和QUARK是目前最成功的預測方法之一,主要是因為DHMA算法使用了QUARK預測的距離譜,同時DHMA算法的初始階段和變異過程中都采用了Rosetta的能量函數.表1和表2中的DH、Ro和QU分別是DHMA、Rosetta和QUARK的簡稱,其中QUARK的預測結果直接通過QUARK在線服務器直接得到.DHMA算法與Rosetta的預測結果主要分3步得到:

表1 第一個聚類結構的TM-core
Table 1 TM-core of first cluster structure

No.PDB IDSizeTypeDDROQU11ail70α0.720.390.3821bkrA108α0.530.340.3231bgf123α0.310.250.2841hz6A64α/β0.510.590.6151c8cA64α/β0.260.630.2861ctf125α/β0.450.490.3871vcc77α/β0.410.340.3081acf125α/β0.620.580.8591c9oA66β0.550.390.64101gvp87β0.340.290.33

1)在構象池中,根據能量函數選取能量最小的前10000個結構;

2)利用K-均值聚類算法進行聚類,選出聚類中心能量最小的一個類;

3)選出與類心最近的5個結構作為最后的預測結果.在聚類過程中,聚類半徑為能量值,聚類簇數為5.

表1是DHMA、Rosetta和QUARK聚類得到的第一個結構與天然態結構比對的TM-score值.如表1所示,蛋白1ail、1bkrA、1bgf、1vcc和1gvp的結果表明DHMA的預測精度均優于Rosetta和QUARK,并且其中1ail和1bkrA的預測結果要明顯優于Rosetta和QUARK.但對1acf和1c9oA蛋白來說,QUARK要明顯優于DHMA和Rosetta;對1c8cA和1ctf蛋白來說,Rosetta的預測精度優于DHMA和QUARK.總體而言,DHMA的預測精度優于Rosetta和QUARK,是一種有效的蛋白結構預測方法.

表2 DHMA、Rosetta和QUARK的實驗結果比較
Table 2 Comparison of results of DHMA,Rosetta and QUARK

No.PDBLowest-RMSD(?)Mean-RMSD(?)DHRoQUDHRoQU11ail5.247.395.268.248.227.1621bkrA7.4510.378.7810.1711.0410.7031bgf10.9512.8711.4612.5415.0513.5741hz6A3.372.862.843.693.364.2051c8cA7.844.409.258.744.5410.0761ctf5.574.103.447.037.817.6571vcc7.745.704.207.967.488.3181acf9.208.581.9911.5910.275.6091c9oA4.765.393.385.826.885.31101gvp8.348.5810.4212.4912.9312.74

表2中Lowest-RMSD代表預測得到的最低RMSD值,Mean-RMSD代表預測得到的平均RMSD值.對于測試蛋白1ail、1bkrA、1bgf和1gvp來說,DHMA預測模型的Lowest-RMSD都優于Rosetta和QUARK,尤其是1bkrA和1c8cA,DHMA算法的Lowest-RMSD比QUARK和Rosetta的精度高于1?.雖然,DHMA算法的Lowest-RMSD結果比QUARK預測的差,但是DHMA預測的Lowest-RMSD優于Rosetta預測的結果,其中對于1hz6A、1ctf和1vcc蛋白,DHMA預測的Lowest-RMSD比QUARK的差,但是在Mean-RMSD方面,DHMA算法的預測精度優于QUARK算白,DHMA預測的Lowest-RMSD比QUARK的差,但是在Mean-RMSD方面,DHMA算法的預測精度優于QUARK算法.同時,10個蛋白質中有4個蛋白,DHMA預測的Mean-RMSD優于QUARK和Rosetta預測的結果.總體而言,在Lowest-RMSD方面,DHMA比QUARK算法的預測結果差,但是在Mean-RMSD方面,DHMA算法的預測精度優于Rosetta和QUARK算法.

圖1是DHMA、Rosetta和QUARK預測的結構與天然態結構之間的比對圖.對蛋白1acf來說,DHMA和Rosetta結構在β區域均折疊錯誤,因此 QUARK的結構與天然態的結構更相似,但對1ail來說,DHMA得到的結構與天然態結構比Rosetta和QUARK的結構重疊區域更多,更相似.

3.3 采樣能力分析

在采樣能力分析中,本文將DHMA算法與Rosetta進行比較.部分蛋白的采樣結果如圖2所示,圖中橫坐標為構象與天然態結構之間的RMSD值;縱坐標為每個構象的能量值.

圖1 DHMA、Rosetta和QUARK預測結構與天然態結構的比對Fig.1 Comparison between predicted structure by DHMA、Rosetta and QUARK and native structure

如圖2所示,針對蛋白來1acf、 1bgf和1c9oA,DHMA算法增強了近天然態區域的采樣,尤其對1bgf來說DHMA顯著的增強了低RMSD區域的采樣能力.但是對圖2中的1hz6A蛋白來說,Rosetta比DHMA算法具有更好的采樣能力,這主要是因為現有的能量包含足夠多的信息,對一些序列較短的蛋白來說,現有的能量函數可以有效地指導構象的采樣.

圖2 構象空間采樣比較Fig.2 Comparsion of conformation space sampling

4個蛋白的RMSD分布如圖3所示,圖中橫坐標為構象與天然態結構之間的RMSD值;縱坐標為每個區域內構象數占所有構象總數的百分比.如圖3所示,針對大多數蛋白來說,DHMA算法增強了低RMSD區域的采樣.但對1hz6A蛋白來說,Rosetta比DHMA算法具有更好的RMSD分布,因此對1hz6A蛋白質來說,Rosetta的聚類結果比DHMA的精度高.

總體而言,DHMA算法能夠在構象空間進行更充分的采樣,采樣得到蛋白質近天然態構象.這主要是DHMA算法采用了距離分布和疏水模型指導構象的采樣,有效的彌補了能量函數不精確帶來的誤差.

圖3 RMSD分布比較Fig.3 Comparison of RMSD distribution

4 結 論

本文提出一種距離和疏水模型輔助的蛋白質結構預測方法.在進化算法的框架下,首先,根據親疏水氨基酸的旋轉半徑差來指導種群的變異;然后,構建距離分布模型和疏水概率模型來指導種群的更新,進而提高預測的精度.DHMA算法通過回轉半徑和疏水概率模型,有效的增強了疏水性在蛋白質的折疊過程中的驅動;同時通過構建距離分布模型,進一步加強了遠程作用力的引導.10個測試蛋白的實驗結果表明,DHMA算法具有較好的搜索性能和預測精度,是一種有效的構象空間搜索算法.

在下一步的研究中,我們將利用機器學習以及數據挖掘技術來預測高質量的距離譜,并且將其與疏水特性融合來改進現有的能量函數,以期得到更好的預測結果.同時,在以后的研究中,將會采用其他優化方法來進一步提高預測精度,如多目標優化技術和多模態優化方法.