李 亭，劉 俊，周曉根，張貴軍

1(浙江工業(yè)大學 信息工程學院，杭州 310023)2(密歇根大學 計算醫(yī)學與生物信息學系，安娜堡 48109)

1 引 言

蛋白質(zhì)是生命活動的主要承擔者[1].蛋白質(zhì)的三維結(jié)構(gòu)則決定了其所具有的特定功能[2].獲取蛋白質(zhì)的三維結(jié)構(gòu)是研究其生物機理的基礎(chǔ)，對蛋白質(zhì)功能注釋、疾病研究、藥物設計等方面有著重要的意義[3，4].目前主要通過生物實驗測定方法獲取蛋白質(zhì)的三維結(jié)構(gòu)，包括X-射線晶體衍射、核磁共振和冷凍電鏡技術(shù).這些實驗測定方法耗時長且成本高昂[5-6].同時隨著基因測序技術(shù)的快速發(fā)展，已測定蛋白質(zhì)序列數(shù)目和已測定蛋白質(zhì)結(jié)構(gòu)數(shù)目之間存在巨大差距，并且這種差距正在逐年擴大.在理論探索和應用需求的推動下，根據(jù)Anfinsen法則[7]，從氨基酸序列出發(fā)，利用計算機技術(shù)結(jié)合優(yōu)化算法的從頭蛋白質(zhì)三維結(jié)構(gòu)預測成為了生物信息學領(lǐng)域的研究熱點[8].

從頭蛋白質(zhì)結(jié)構(gòu)預測主要有兩大難點：1)構(gòu)建合適的能量函數(shù)，引導構(gòu)象向天然態(tài)結(jié)構(gòu)(即能量最低點狀態(tài))折疊；2)設計高效的構(gòu)象搜索方法，在復雜的能量景觀中探索低能量區(qū)域[9].傳統(tǒng)的能量函數(shù)[10，11]是基于理化模型構(gòu)建的，沒有考慮到不同蛋白質(zhì)的特殊性，因此很難準確地引導構(gòu)象向天然態(tài)結(jié)構(gòu)折疊.近幾年，得益于已測定蛋白質(zhì)序列的快速增長和深度學習技術(shù)的不斷完善，根據(jù)目標序列的進化信息提取先驗知識輔助結(jié)構(gòu)預測取得了巨大進展[12-17].尤其是殘基-殘基接觸預測，為蛋白質(zhì)折疊提供了重要的約束信息，很大程度上彌補了能量函數(shù)不精確造成的影響.許多利用殘基-殘基接觸預測蛋白質(zhì)結(jié)構(gòu)的方法被提出[18-21]，其中CONFOLD[22]將殘基-殘基接觸和二級結(jié)構(gòu)信息轉(zhuǎn)化為距離、二面角和氫鍵約束，進行蛋白質(zhì)結(jié)構(gòu)預測.自CASP13以來，基于深度學習的殘基間距離預測迅速成為新的研究熱點[23]，殘基間距離預測描述了殘基對在不同距離區(qū)間的概率，相比殘基-殘基接觸包含了更多的幾何約束信息，更有利于結(jié)構(gòu)預測[24].AlphaFold[25]根據(jù)預測的殘基間距離分布構(gòu)建蛋白質(zhì)特定的勢能函數(shù)，并利用隨機梯度下降尋找最優(yōu)結(jié)構(gòu)，最終，AlphaFold在CASP13中脫穎而出.數(shù)十年來，研究學者們提出了大量方法進行蛋白質(zhì)構(gòu)象空間搜索，其中蒙特卡洛[26-28]和副本交換蒙特卡洛[29]取得了巨大成功，并成功的應用于國際知名的從頭預測服務器Rosetta[30，31]和QUARK[32，33]，而進化算法[34-40]、多模態(tài)優(yōu)化[41]、分子動力學[42-44]模擬也展現(xiàn)出各自的優(yōu)勢.在差分進化算法框架下，SCDE[45]利用接觸圖和二級結(jié)構(gòu)提出策略來預測蛋白質(zhì)結(jié)構(gòu)；ItFix[46]建立粗粒度模型，給定適當?shù)乃阉鞑呗院湍芰亢瘮?shù)進行結(jié)構(gòu)預測.片段組裝[47-49]策略被廣泛應用于從頭蛋白質(zhì)結(jié)構(gòu)預測方法，能夠有效地減小構(gòu)象搜索空間，提高采樣效率.

針對上述問題，本文提出了一種基于距離約束和二面角優(yōu)化的蛋白質(zhì)結(jié)構(gòu)預測方法(DCDA).首先，對預測的殘基對距離分布信息進行篩選，構(gòu)建基于殘基對距離分布的構(gòu)象評估模型，結(jié)合片段組裝技術(shù)大范圍搜索構(gòu)象空間；進而，構(gòu)建Loop區(qū)域特定的局部構(gòu)象評估模型，并結(jié)合基于二面角的差分進化采樣策略，增強結(jié)構(gòu)靈活的Loop區(qū)域采樣，突破片段庫的約束，進一步提高預測模型的精度.15個測試蛋白質(zhì)的實驗結(jié)果表明，DCDA在預測精度上優(yōu)于Rosetta、QUARK和CONFOLD，是一種有效的蛋白質(zhì)結(jié)構(gòu)預測方法.

2 DCDA算法設計

2.1 殘基對距離特征提取

預測的殘基對距離分布相比預測的殘基對接觸包含了更詳細的幾何約束信息，可以有效地減少構(gòu)象搜索空間，更準確、快速地引導蛋白質(zhì)折疊.為了避免冗余信息的影響，本文根據(jù)距離預測概率峰值對距離矩陣進行過濾.首先，去除掉殘基對序列分離小于6的殘基對；然后，對于預測的距離分布中任意殘基對之間距離小于等于2的兩個殘基對，僅保留預測的距離概率峰值最大的殘基對，兩個殘基對之間的距離計算公式如下：

(1)

其中，(i1，j1)表示兩個殘基對中的第1個殘基對；(i2，j2)表示兩個殘基對中的第2個殘基對.

圖1為蛋白質(zhì)1E2A_A過濾前后的殘基對距離分布圖.

2.2 基于殘基對距離分布的構(gòu)象評估模型

預測的殘基對距離分布包含的是離散信息，即殘基對距離分別在[2，2.5)?、[2.5，3)?、…、[19.5，20)?區(qū)間內(nèi)的概率.然而，結(jié)構(gòu)預測過程中構(gòu)象殘基間的真實距離是連續(xù)的，為了更精確的評估構(gòu)象，本文將離散的距離分布信息轉(zhuǎn)換成連續(xù)的分布函數(shù)，確保對于每一個給定的距離都有與之對應的概率值.

本文采用三次樣條插值為每一個殘基對構(gòu)造連續(xù)的距離分布函數(shù).令每個預測距離區(qū)間的中值為三次樣條插值的節(jié)點，分別在兩端添加一個邊界點，即節(jié)點分別為x0=0，x1=2.25，x2=2.75，…，x36=19.75，x37=∞；以38個節(jié)點為依據(jù)，構(gòu)造37個三次函數(shù)fk(x)：

fk(x)=ak+bkx+ckx2+dkx3，k={0，1，2，…，35，36}

(2)

其中，ak，bk，ck，dk為待定系數(shù)，共有148個待定系數(shù).則殘基對的距離分布函數(shù)可以表示為：

(3)

(4)

為了確定待定系數(shù)值，需構(gòu)造148個方程進行求解.

圖1 過濾前后的殘基對距離分布圖Fig.1 Comparison distance map before and after filtering

首先，要求每一個分段函數(shù)經(jīng)過對應的節(jié)點：

(5)

其中，pk表示節(jié)點xk對應的距離區(qū)間的概率值，令p0=0、p37=0；由公式(5)可構(gòu)造74個方程.

(6)

可構(gòu)造72個方程.最后，通過邊界條件再構(gòu)造2個方程.通過求解148個方程得到系數(shù)值，得到方程表達式.

定義殘基對(i，j)對應的距離分布函數(shù)為f(i，j)(x)，x表示殘基對(i，j)在構(gòu)象中的真實距離.根據(jù)三次樣條得到的方程來計算采樣過程中構(gòu)象的勢能，構(gòu)建構(gòu)象評估模型，以引導構(gòu)象向天然態(tài)折疊.

根據(jù)所有殘基對的距離分布函數(shù)，構(gòu)建構(gòu)象評估模型：

(7)

其中Q是過濾后的殘基對集合，dij為待評估構(gòu)象的第i號殘基和第j號殘基的Cβ原子(甘氨酸為Cα原子)間的歐氏距離，f(i，j)(dij)是x=dij時f(i，j)(x)的函數(shù)值，表示殘基對(i，j)的距離為dij的概率.

2.3 片段組裝

由于蛋白質(zhì)構(gòu)象空間的高維特性，在巨大的構(gòu)象空間中進行采樣是不合適的.片段組裝技術(shù)利用已知結(jié)構(gòu)的局部信息，將每一個殘基的二面角約束在一組離散值內(nèi)，從而極大地縮小了構(gòu)象搜索空間.對于給定目標構(gòu)象，片段組裝過程如下：首先，隨機選擇一個滑動窗口，在該窗口對應的片段庫中隨機選擇一個片段替換窗口內(nèi)原有的片段，如圖2所示，L為序列長度；然后，利用公式(7)計算目標構(gòu)象和組裝得到的構(gòu)象的得分，并根據(jù)boltzmann準則判斷是否組裝成功；如果成功，則用組裝后的構(gòu)象替換目標構(gòu)象；否則，再次隨機選擇窗口進行片段組裝，直至組裝成功或組裝次數(shù)達到200次.

圖2 片段組裝示意圖Fig.2 Schematic diagram of fragment assembly

2.4 基于二面角的差分進化采樣策略

本文提出了基于差分進化算法的Loop區(qū)域二面角采樣策略，在利用片段庫信息的同時，避免結(jié)構(gòu)靈活的Loop區(qū)域受到片段庫的約束，探索更多結(jié)構(gòu)合理的構(gòu)象.首先，利用DSSP算法[50]計算目標構(gòu)象的二級結(jié)構(gòu)，隨機選擇一段二級結(jié)構(gòu)為Loop的局部區(qū)域，并構(gòu)建Loop區(qū)域特定的構(gòu)象評估模型；然后，以窗口寬度為3個殘基，步長為1個殘基，從選定Loop區(qū)域的起始殘基位滑動到結(jié)束殘基位，形成K個滑動窗口.對每個窗口內(nèi)的片段執(zhí)行變異、交叉操作生成候選構(gòu)象，添加到候選構(gòu)象池；最終，通過選擇操作，從候選構(gòu)象池中選擇最優(yōu)構(gòu)象替換目標構(gòu)象.具體過程如下：

a)構(gòu)建局部構(gòu)象評估模型.針對選定的Loop區(qū)域，從距離分布圖中選取殘基編號分別位于選定區(qū)域左右兩端的殘基對構(gòu)建局部評估模型.計算公式如下：

(8)

其中，M是選定區(qū)域左右兩端的殘基對集合，wij表示殘基對(i，j)的權(quán)重，由殘基i和殘基j所在二級結(jié)構(gòu)與選定Loop區(qū)域的相對位置關(guān)系決定的.計算公式如下：

(9)

其中，indexi表示殘基對中殘基i所在的二級結(jié)構(gòu)片段與選定Loop結(jié)構(gòu)之間非Loop結(jié)構(gòu)的數(shù)量；indexj表示殘基對中殘基j所在的二級結(jié)構(gòu)片段與選定Loop結(jié)構(gòu)之間非Loop結(jié)構(gòu)的數(shù)量.

(10)

其中，F(xiàn)是縮放因子，一般在[0，2]之間選擇，本文取0.5.

(11)

其中，CR表示交叉概率，l表示向量中角度的索引值，lrand表示一個隨機整數(shù)，它的取值范圍是[1，6].

d)選擇操作.根據(jù)公式(8)計算目標構(gòu)象和候選構(gòu)象池中每個構(gòu)象的得分，按得分由低至高選取候選構(gòu)象池中的構(gòu)象替換目標構(gòu)象，直至替換成功或遍歷完候選構(gòu)象池中所有構(gòu)象.候選構(gòu)象替換目標構(gòu)象的公式如下：

(12)

(13)

其中，KT表示溫度常數(shù).

2.5 算法描述

DCDA算法的流程描述如下：

輸入：目標蛋白序列，片段庫，殘基預測距離

輸出：蛋白質(zhì)三維結(jié)構(gòu)

1. 初始化種群T←Rosetta stage1，stage2

2.forg←1 to G do

3.fori←1 toNP

4.Ti←片段組裝

5. 隨機選取一段Loop區(qū)域，構(gòu)建局部構(gòu)象評估模型

6. 創(chuàng)建滑動窗口{W1，W2，…，WK}

7.fork←1 toK

11.endfor

13.fork←1 to K

16.break；

17.endif

18.endfor

19.endfor

20.endfor

21.return種群中得分最低的構(gòu)象

3 實驗結(jié)果與分析

3.1 測試蛋白和實驗設置

為了驗證算法的有效性，從PDB庫中選取序列長度范圍從72到147的15個蛋白質(zhì)測試DCDA的性能.這15個測試蛋白質(zhì)的折疊類型包括α、β和α/β，詳細信息如表1中的第2-4列所示.這些蛋白質(zhì)的三維結(jié)構(gòu)已經(jīng)由生物實驗方法測定，且被廣泛用于蛋白質(zhì)測試[51，52]，可以較好地測試DCDA的性能.

表1 測試蛋白的TM-score
Table 1 TM-score of test protein

No.PDBLenTypeDCDARosQUACONF11G8Q_A90α0.620.30.360.6221H4L_D147α0.590.270.420.6531W53_A84α0.70.420.370.6642BL7_A79α0.840.460.60.7951B4B_A72α/β0.820.330.410.6261CF7_B82α/β0.750.440.510.6571CQA_A123α/β0.680.270.350.5581DJ7_A109α/β0.550.270.380.3791IPI_A114α/β0.730.310.550.51101MWP_A96α/β0.510.260.340.3111PZW_A80α/β0.570.290.370.45122H8E_A120α/β0.660.30.340.57133CHB_D103α/β0.330.290.30.32141OK0_A74β0.630.260.220.22152BT9_A90β0.720.340.410.53Average0.6470.3210.3950.521

測試蛋白質(zhì)的片段庫和殘基間距離分布信息分別從ROBETTA[53]和trRosetta[54]服務器獲取，并且在構(gòu)建片段庫和距離分布信息的過程中均已去除同源模板.實驗參數(shù)設置如下：種群規(guī)模NP=100、種群更新次數(shù)G=1000、溫度常數(shù)KT=2、交叉概率CR=0.5、縮放因子F=0.5.在種群初始化過程中采用Rosetta score0和score1能量函數(shù)引導構(gòu)象搜索.在片段組裝和基于二面角的差分進化過程中分別采用構(gòu)象評估模型Sd和局部評估模型Sloop_d引導構(gòu)象采樣.通過計算預測模型和天然態(tài)結(jié)構(gòu)的Cα原子的均方根偏差(RMSD)和TM-score[55]來評估算法的預測精度.

3.2 實驗結(jié)果分析

在測試蛋白質(zhì)上，將DCDA算法的預測結(jié)果與Rosetta、QUARK和CONFOLD預測的結(jié)果進行了比較.Rosetta、QUARK和CONFOLD是蛋白質(zhì)結(jié)構(gòu)預測領(lǐng)域內(nèi)認可度最高的方法之一，與其進行比較能更客觀的反映DCDA的性能.Rosetta和QUARK均采用了片段組裝技術(shù)，是經(jīng)典的兩個蛋白質(zhì)結(jié)構(gòu)預測服務器；CONFOLD則利用了殘基接觸、二級結(jié)構(gòu)等先驗信息構(gòu)建幾何約束，代表著新興的預測算法.

表1和表2分別列出了DCDA、Rosetta、QUARK和CONFOLD在15個測試蛋白質(zhì)上預測模型的TM-score和RMSD，其中Ros、QUA和CONF分別代表Rosetta、QUARK和CONFOLD.其中Rosetta的預測結(jié)果是通過安裝本地版并使用默認參數(shù)預測得到的，QUARK和CONFOLD的預測結(jié)果是直接通過其在線服務器預測直接得到的.表中加粗的數(shù)據(jù)表明在對應指標下，相應算法的預測精度最優(yōu).

如表1所示，DCDA在15個測試蛋白質(zhì)中的14個上獲得了TM-score最高的預測模型，尤其對于2BL7_A和1B4B_A，DCDA預測模型的TM-score均在0.8以上；DCDA預測出了14個測試蛋白質(zhì)的正確折疊(TM-score≥0.5)，占全部測試蛋白的93.3%，是所有對比方法中最高的；在平均TM-score上，相比于其余3種方法中預測精度最高的CONFOLD提高了24.2%.由此可見，DCDA預測模型的TM-score明顯優(yōu)于Rosetta、QUARK和CONFOLD.表2列出了4種方法預測模型與天然態(tài)結(jié)構(gòu)的RMSD.DCDA在14個測試蛋白質(zhì)上預測出了RMSD最小的結(jié)構(gòu)模型，平均RMSD為4.288?，相比較于CONFOLD的平均RMSD降低了46.1%.尤其是蛋白1OK0_A，DCDA的RMSD相比QUARK降低了7.01?.由此可見，在RMSD方面，DCDA算法的預測精度優(yōu)于Rosetta、QUARK和CONFOLD算法.整體而言，DCDA在TM-score和RMSD評價指標上均優(yōu)于對比算法，是一種有效的蛋白質(zhì)結(jié)構(gòu)預測算法.

表2 測試蛋白的RMSD
Table 2 RMSD of test protein

No.PDBDCDARosQUACONF11G8Q_A4.0611.4310.654.521H4L_D5.215.9210.894.8731W53_A3.2511.1312.075.3642BL7_A2.037.53.372.4651B4B_A2.859.328.163.9861CF7_B2.797.386.363.2171CQA_A4.2613.558.026.0881DJ7_A6.3413.4210.6612.2691IPI_A3.3611.834.8514.52101MWP_A6.1513.5414.8813.57111PZW_A5.3911.137.519.03122H8E_A3.7814.3510.348.87133CHB_D8.9713.3512.3214.11141OK0_A3.3510.7210.3611.56152BT9_A2.5410.126.824.96Average4.28811.6469.1517.956

為了驗證基于二面角的差分進化采樣策略的有效性，本文設計了未使用基于二面角的差分進化采樣策略的對比實驗DC，對比結(jié)果如表3所示.與DC算法相比，DCDA在15個測試蛋白質(zhì)中的11個蛋白質(zhì)上獲得了RMSD最小的結(jié)構(gòu)模型，在13個蛋白質(zhì)上獲得了TM-score最高的預測模型.對比表2和表3的結(jié)果可以看出，僅使用殘基間距離分布信息的DC算法，在15個測試蛋白質(zhì)中，14個測試蛋白質(zhì)的TM-score比Rosetta高，其平均TM-score值比Rosetta高0.275，其平均RMSD比Rosetta低6.267?.

表3 DCDA組件比較
Table 3 Comparison of the DCDA component

No.PDBDCDARMSDTM-scoreDCRMSDTM-score11G8Q_A4.060.626.040.5121H4L_D5.20.594.560.6331W53_A3.250.72.40.7742BL7_A2.030.841.80.8251B4B_A2.850.823.280.761CF7_B2.790.753.890.6771CQA_A4.260.685.040.6181DJ7_A6.340.5511.750.4691IPI_A3.360.733.990.68101MWP_A6.150.518.690.5111PZW_A5.390.573.790.62122H8E_A3.780.665.150.55133CHB_D8.970.3311.420.29141OK0_A3.350.635.460.43152BT9_A2.540.723.420.7Average3.9950.6615.3790.596

實驗數(shù)據(jù)表明基于二面角的差分進化采樣策略可以有效探索片段庫約束空間之外的構(gòu)象，提高結(jié)構(gòu)靈活的Loop區(qū)域結(jié)構(gòu)的多樣性，進而提高預測模型的精度.此外，對比表3中DC與表1和表2中Rosetta、QUARK、CONFOLD的結(jié)果可以發(fā)現(xiàn)，僅使用殘基間距離分布信息的方法預測的結(jié)果也明顯優(yōu)于對比的方法，進一步表明預測的殘基間距離分布信息能夠更好的引導構(gòu)象采樣；并且，結(jié)合基于二面角的差分進化采樣策略能夠進一步提升預測精度.

圖3 DCDA預測結(jié)構(gòu)與天然態(tài)結(jié)構(gòu)的比對Fig.3 Comparison between predicted structure by DCDA and native structure

圖3是DCDA預測的模型與天然態(tài)結(jié)構(gòu)之間的比對圖，其中包含α和α/β折疊類型.從圖中可以發(fā)現(xiàn)DCDA預測的蛋白質(zhì)結(jié)構(gòu)模型與天然態(tài)結(jié)構(gòu)之間有著較好的重疊，和天然態(tài)結(jié)構(gòu)相似，具有很好的匹配度.

3.3 算法采樣能力分析

本節(jié)對DCDA和Rosetta的采樣能力進行了分析和對比.部分蛋白的采樣分布如圖4所示，圖中橫坐標為搜索過程中產(chǎn)生的構(gòu)象與天然態(tài)結(jié)構(gòu)之間的RMSD，縱坐標為每個區(qū)域內(nèi)構(gòu)象數(shù)占所有構(gòu)象總數(shù)的百分比，實線是DCDA的采樣分布，虛線是Rosetta的采樣分布.可以發(fā)現(xiàn)，DCDA在4個測試蛋白上的近天然態(tài)采樣比例均高于Rosetta，DCDA采樣概率峰值的RMSD比Rosetta更接近零，這表明DCDA的近天然態(tài)采樣能力強于Rosetta.對于蛋白質(zhì)1W53_A和2BL7_A，DCDA分別有89.95%和91.3%的過程構(gòu)象的RMSD小于5?，而Rosetta在1W53_A上的過程構(gòu)象的RMSD均大于5?，在2BL7_A上的過程構(gòu)象的RMSD僅有5.51%小于5?，對于1CQA_A，盡管DCDA和Rosetta的近天然態(tài)采樣能力都不強，但是DCDA的采樣分布比Rosetta更接近天然態(tài)構(gòu)象.

圖4 RMSD分布比較Fig.4 Comparison of RMSD distributions

總體而言，DCDA在搜索過程中能夠采樣到更多近天然態(tài)構(gòu)象.主要是DCDA算法采用了距離分布信息和基于二面角的差分進化采樣策略，有效的提升了預測精度.

4 結(jié) 論

本文提出一種距離約束和二面角優(yōu)化的蛋白質(zhì)結(jié)構(gòu)預測方法.首先，對預測的距離分布圖進行篩選；然后，利用三次樣條插值根據(jù)預測的蛋白質(zhì)殘基間距離分布構(gòu)建構(gòu)象評估模型，結(jié)合片段組裝技術(shù)，大規(guī)模采樣構(gòu)象空間；進而，在進化算法的框架下，針對結(jié)構(gòu)靈活的Loop區(qū)域，設計基于二面角的差分進化采樣策略，構(gòu)建局部構(gòu)象評估模型，利用基于滑動窗口的二面角交叉變異生成多個候選構(gòu)象，并根據(jù)局部構(gòu)象評估模型實現(xiàn)構(gòu)象更新.基于二面角的差分進化采樣策略能夠突破片段庫的約束，采樣到更多近天然態(tài)構(gòu)象，提高算法的預測精度.在15個測試蛋白質(zhì)上的實驗結(jié)果表明，DCDA算法具有較強的搜索性能和較高的預測精度，是一種有效的構(gòu)象空間搜索算法.在下一步研究中，我們將利用深度學習技術(shù)構(gòu)建深度殘差神經(jīng)網(wǎng)絡預測蛋白質(zhì)殘基二面角，根據(jù)預測的二面角和距離信息構(gòu)建約束函數(shù)，進一步提高蛋白質(zhì)結(jié)構(gòu)預測精度.