基于多目標(biāo)優(yōu)化的蛋白質(zhì)三維結(jié)構(gòu)預(yù)測(cè)

王雨林，張 彪，沈紅斌

(上海交通大學(xué) 圖像處理與模式識(shí)別研究所, 上海 200240)

蛋白質(zhì)在生命細(xì)胞中發(fā)揮著各種各樣的功能，是生命活動(dòng)的主要承擔(dān)者之一[1-2]．自然界的蛋白質(zhì)主要由20種基本氨基酸構(gòu)成，存在著一條或多條肽鏈，長(zhǎng)度從幾十到幾千不等．蛋白質(zhì)的結(jié)構(gòu)基本決定了其特別的生物化學(xué)特性，自20世紀(jì)50年代起圍繞求解蛋白質(zhì)三維結(jié)構(gòu)的研究獲得了廣泛的關(guān)注[3-5]．但由于實(shí)驗(yàn)的困難性[6]，從蛋白質(zhì)序列出發(fā)利用計(jì)算機(jī)算法預(yù)測(cè)每個(gè)氨基酸的空間坐標(biāo)并進(jìn)而最終實(shí)現(xiàn)三維結(jié)構(gòu)的預(yù)測(cè)獲得廣泛重視．總體而言，當(dāng)前基于氨基酸序列來預(yù)測(cè)蛋白質(zhì)結(jié)構(gòu)的方法大致分為3類[7-8]：同源性建模法、穿線法以及從頭預(yù)測(cè)法．

同源建模法是一種十分成功的蛋白質(zhì)結(jié)構(gòu)預(yù)測(cè)方法，其基本的假設(shè)是結(jié)構(gòu)的變化率遠(yuǎn)低于序列的變化率，氨基酸序列包含了預(yù)測(cè)蛋白質(zhì)結(jié)構(gòu)的充足信息[9]．基于以上思想，同源建模法通過比對(duì)序列與PDB(protein data bank)數(shù)據(jù)庫(kù)中已經(jīng)實(shí)驗(yàn)求解獲得結(jié)構(gòu)的模板蛋白質(zhì)序列的相似性[10-13]，當(dāng)相似性高于一定的閾值時(shí)，說明它們可能具有相似的空間結(jié)構(gòu)，并用模板的結(jié)構(gòu)作為待預(yù)測(cè)的結(jié)構(gòu)基礎(chǔ)，在施加部分相關(guān)約束的同時(shí)優(yōu)化缺失的重原子與支鏈，最終優(yōu)化獲得預(yù)測(cè)目標(biāo)結(jié)構(gòu)．當(dāng)能找到高可靠模板時(shí)，這類同源建模法能取得很高的準(zhǔn)確性，但當(dāng)無法尋找到同源性蛋白時(shí)，此方法難以取得好的效果．

穿線法在兩個(gè)蛋白質(zhì)即使序列相似性很低，也可能具有相同的結(jié)構(gòu)方面取得了較大的研究進(jìn)展[14]．自然界中存在著數(shù)量很大的氨基酸序列，但是它們對(duì)應(yīng)的不同三維結(jié)構(gòu)卻比較有限．在進(jìn)行蛋白質(zhì)結(jié)構(gòu)預(yù)測(cè)時(shí)，相比于直接進(jìn)行序列與序列的比對(duì)，序列與三維結(jié)構(gòu)的比對(duì)往往能取得更好的效果[15]．穿線法大致可以分為兩類[16-18]：一類提取蛋白質(zhì)的特征，如蛋白質(zhì)的二級(jí)結(jié)構(gòu)、氨基酸殘基的可溶性等，來構(gòu)建特征向量作為蛋白質(zhì)的特征表示，然后用動(dòng)態(tài)規(guī)劃等相似算法將序列與當(dāng)前的特征進(jìn)行比對(duì)得到序列對(duì)應(yīng)的三維結(jié)構(gòu)；另一類方法基于統(tǒng)計(jì)勢(shì)能的方式，如統(tǒng)計(jì)在特定距離上兩個(gè)氨基酸殘基對(duì)出現(xiàn)的頻率，來評(píng)估序列與蛋白質(zhì)結(jié)構(gòu)的相似性，從而預(yù)測(cè)獲得三維結(jié)構(gòu)．

盡管基于同源建模法與穿線法的蛋白質(zhì)三維結(jié)構(gòu)預(yù)測(cè)取得了很大的成功，但是仍然有許多重要的蛋白質(zhì)因?yàn)闊o法找到合適的模板結(jié)構(gòu)而難以取得好的效果，此時(shí)需要采取從頭預(yù)測(cè)法的方法．基于從頭預(yù)測(cè)法的研究一般分為兩類[19]：① 通過模擬蛋白質(zhì)的折疊過程，從隨機(jī)構(gòu)象出發(fā)優(yōu)化預(yù)先給定的能量函數(shù)至最低狀態(tài)；② 基于碎片組裝的方式進(jìn)行結(jié)構(gòu)預(yù)測(cè)[20]．兩種方法的的核心都在于獲得能準(zhǔn)確評(píng)估結(jié)構(gòu)的能量函數(shù)[21]和高效的構(gòu)象空間搜索方法．第一類方法目前廣泛采用分子動(dòng)力學(xué)模擬的方式，通過解析天然蛋白質(zhì)折疊時(shí)的機(jī)理與相互作用來指導(dǎo)預(yù)測(cè)序列的折疊與構(gòu)象生成．然而一方面模擬蛋白質(zhì)折疊的過程十分耗時(shí)[22]，并且即使投入大量的算力，對(duì)長(zhǎng)度較長(zhǎng)的蛋白也難以達(dá)到高的精確度；另一方面，實(shí)驗(yàn)?zāi)M蛋白質(zhì)折疊的過程本身與自然界的天然折疊方式存在偏差，生成的結(jié)構(gòu)與天然結(jié)構(gòu)可能會(huì)有諸多不同．

為了解決上述問題，以ROSETTA[23]為代表的第二類方法通過引入構(gòu)象的約束來減少計(jì)算的復(fù)雜性．ROSETTA基于以下一個(gè)重要假設(shè)[23]：短序列片段具有很強(qiáng)的局部結(jié)構(gòu)偏好，并且這種偏好的強(qiáng)度與多樣性高度依賴于序列本身．ROSETTA首先計(jì)算序列的二級(jí)結(jié)構(gòu)，蛋白質(zhì)的二級(jí)結(jié)構(gòu)整體規(guī)劃了多肽鏈沿主鏈骨架方向的空間排列和走向；然后以二級(jí)結(jié)構(gòu)為約束條件，將序列以滑動(dòng)窗口的方式與蛋白質(zhì)數(shù)據(jù)庫(kù)中的結(jié)構(gòu)進(jìn)行相似性比對(duì)，數(shù)據(jù)庫(kù)中的結(jié)構(gòu)具有坐標(biāo)信息，與序列有高相似度的結(jié)構(gòu)坐標(biāo)可以轉(zhuǎn)化后作為序列的坐標(biāo)．算法通過生成具有坐標(biāo)的碎片的方式，隨機(jī)在序列上選取位點(diǎn)替換為備選的碎片結(jié)構(gòu)，來改善構(gòu)象的結(jié)構(gòu)，并使用單個(gè)能量函數(shù)進(jìn)行評(píng)估，朝全局自由能降低的方向優(yōu)化構(gòu)象，最終得到能量最低的結(jié)構(gòu)．ROSETTA 算法取得了很好的效果[23]，但是在部分蛋白的預(yù)測(cè)上不能找到全局最優(yōu)解，其中一個(gè)重要的潛在原因是僅使用單個(gè)能量函數(shù)進(jìn)行優(yōu)化．由于能量函數(shù)并非完全精確，對(duì)構(gòu)象的評(píng)價(jià)和構(gòu)象本身與天然結(jié)構(gòu)的相似度不是完全的正相關(guān)，在優(yōu)化過程中存在著舍棄掉最優(yōu)解的情況．

正是考慮到很難用統(tǒng)一的單能量函數(shù)衡量所有構(gòu)象的質(zhì)量，文中采用多個(gè)能量函數(shù)優(yōu)化構(gòu)象的TRIOFOLD 方法，多個(gè)能量函數(shù)的引入能進(jìn)一步擴(kuò)大搜索空間，對(duì)于構(gòu)象的評(píng)價(jià)將更加全面，能幫助優(yōu)化搜索的方向，也有助于克服能量函數(shù)本身的精度問題．實(shí)驗(yàn)結(jié)果驗(yàn)證了TRIOFOLD方法的可行性，該方法有望成為領(lǐng)域已有算法的有益補(bǔ)充．

1 TRIOFOLD多目標(biāo)優(yōu)化算法

文中采用從頭預(yù)測(cè)法的方法，通過碎片替換的方式對(duì)生成的構(gòu)象進(jìn)行更新，并利用多目標(biāo)優(yōu)化的方法對(duì)生成的構(gòu)象進(jìn)行優(yōu)化，最終在迭代終止時(shí)生成蛋白質(zhì)結(jié)構(gòu)的非支配集合，算法流程如圖1．

圖1 TRIOFOLD多目標(biāo)優(yōu)化算法流程Fig.1 Flowchart of TRIOFOLD multi-objective optimization algorithm

1.1 能量函數(shù)

預(yù)測(cè)過程中，蛋白質(zhì)可以通過能量函數(shù)得到一個(gè)表征自身結(jié)構(gòu)質(zhì)量的數(shù)值，分別引入3個(gè)能量函數(shù)ROSETTA、CHARMM[24]、RWPLUS[25]對(duì)程序生成的臨時(shí)構(gòu)象進(jìn)行評(píng)價(jià)篩選．CHARMM是基于物理學(xué)意義的能量函數(shù)，其定義了一套廣泛應(yīng)用于分子動(dòng)力學(xué)的力場(chǎng)，可以對(duì)蛋白質(zhì)結(jié)構(gòu)的鍵長(zhǎng)、鍵角、原子之間的相互作用等進(jìn)行定量化計(jì)算，得到一系列的數(shù)值來評(píng)估蛋白質(zhì)結(jié)構(gòu)的真實(shí)性；RWPLUS是基于統(tǒng)計(jì)學(xué)意義的能量函數(shù)，其使用隨機(jī)游走鏈作為參考狀態(tài)對(duì)蛋白質(zhì)結(jié)構(gòu)進(jìn)行評(píng)價(jià)；ROSETTA是以上兩種類型的綜合，它的一部分能量項(xiàng)計(jì)算力的作用，另一部分能量項(xiàng)基于統(tǒng)計(jì)得到，最后將所有能量項(xiàng)數(shù)值加權(quán)求和作為蛋白質(zhì)評(píng)估的分?jǐn)?shù)．

1.2 基于單目標(biāo)函數(shù)的構(gòu)象優(yōu)化

ROSETTA算法首先計(jì)算序列的二級(jí)結(jié)構(gòu)信息和序列與蛋白質(zhì)數(shù)據(jù)庫(kù)中數(shù)據(jù)的相似信息作為約束，在已知結(jié)構(gòu)的序列碎片數(shù)據(jù)庫(kù)中為序列的所有子序列選定具有高相似性的短碎片結(jié)構(gòu)．然后根據(jù)序列信息生成初始構(gòu)象，選取長(zhǎng)度固定的滑動(dòng)窗口隨機(jī)在序列上選擇一個(gè)局部子序列，使用當(dāng)前子序列的備選碎片結(jié)構(gòu)替換原始結(jié)構(gòu)，生成新的臨時(shí)構(gòu)象．在模擬退火算法下，使用ROSETTA的能量函數(shù)對(duì)構(gòu)象進(jìn)行評(píng)估，基于Metropolis準(zhǔn)則[26]決定臨時(shí)構(gòu)象的保留或是舍棄．不斷迭代以上過程，直到達(dá)到規(guī)定的接受次數(shù)，輸出最終的結(jié)構(gòu)．

具體而言，對(duì)于新生成的構(gòu)象采用單能量函數(shù)進(jìn)行評(píng)估，比較當(dāng)前構(gòu)象與上一次構(gòu)象的能量函數(shù)值，計(jì)算出新構(gòu)象的可行性概率，按一定概率接受新的構(gòu)象：

(1)

式中：x(n)為當(dāng)前構(gòu)象；x(n+1)為生成的新構(gòu)象；E(n)和E(n+1)分別為對(duì)應(yīng)的能量評(píng)估；參數(shù)T為控制模擬退火快慢的調(diào)節(jié)因子．若當(dāng)前構(gòu)象的能量小于上一步構(gòu)象的能量，當(dāng)前的替換操作將會(huì)被接受；若當(dāng)前構(gòu)象的能量大于上一步構(gòu)象的能量，算法將產(chǎn)生一個(gè)在區(qū)間[0,1]的隨機(jī)數(shù)ε，對(duì)于式(1)計(jì)算的概率p，如果ε

ROSETTA的重點(diǎn)之一是最小化序列中殘基局部和全局的交互作用來優(yōu)化生成的構(gòu)象結(jié)構(gòu)，使用單個(gè)能量函數(shù)作為約束提高搜索效率，通過大量替換碎片結(jié)構(gòu)的方式來遍歷構(gòu)象的解空間．但是能量函數(shù)本身是通過對(duì)自然界中的真實(shí)蛋白質(zhì)采用基于統(tǒng)計(jì)或是基于物理力學(xué)建模的方法擬合而成，存在著誤差并非完全精確．對(duì)于單個(gè)能量函數(shù)，以上的優(yōu)化過程容易在單一勢(shì)能方向上過度擬合，使構(gòu)象陷入局部最優(yōu)解中難以跳出．同時(shí)由于單一能量函數(shù)勢(shì)能面狹窄，整體上生成的蛋白質(zhì)相似程度很高，難以得到多樣化的結(jié)構(gòu)．

考慮到這個(gè)問題，領(lǐng)域算法常用的一個(gè)方案是采用加權(quán)求和的方式，把多個(gè)目標(biāo)優(yōu)化函數(shù)轉(zhuǎn)換為單目標(biāo)優(yōu)化問題：

(2)

式中：x為蛋白質(zhì)的構(gòu)象；f1,f2,f3為能量函數(shù)；E1,E2,E3為構(gòu)象的評(píng)估能量值；α1,α2,α3為加權(quán)系數(shù)．能量函數(shù)本身是許多能量項(xiàng)的加權(quán)求和．但是往往權(quán)重因子α1,α2,α3難以確定，對(duì)于權(quán)重較大的能量項(xiàng)，其可能出現(xiàn)的噪聲也將對(duì)整體目標(biāo)函數(shù)產(chǎn)生很大影響．同時(shí)，利用對(duì)能量項(xiàng)進(jìn)行線性加權(quán)的單目標(biāo)優(yōu)化方式，在非凸的解空間尋優(yōu)問題還存在很大的挑戰(zhàn)．

文中將多個(gè)能量作為獨(dú)立的指標(biāo)進(jìn)行使用，實(shí)現(xiàn)多目標(biāo)的平行優(yōu)化流程框架，由于不同能量函數(shù)基于不同的原理得到，且對(duì)于蛋白質(zhì)的三維結(jié)構(gòu)有著各自獨(dú)特的傾向性，使用多個(gè)并行能量函數(shù)同時(shí)對(duì)生成構(gòu)象進(jìn)行評(píng)估能夠有效的緩解單個(gè)能量函數(shù)的精度問題，同時(shí)擴(kuò)大了勢(shì)能面，增加了解的多樣性．多目標(biāo)優(yōu)化問題定義為：

minf(x)=

min[fROSETTA(x),fCHARMM(x),fRWPLUS(x)]

x∈X

(3)

式中：fROSETTA(x)為ROSETTA的能量函數(shù)；fCHARMM(x)為CHARMM的能量函數(shù)；fRWPLUS(x)為RWPLUS的能量函數(shù)；x為蛋白質(zhì)構(gòu)象；X為構(gòu)象空間．多個(gè)能量函數(shù)分別作用于構(gòu)象上，將涉及到處理非支配集合的生成與最優(yōu)解的篩選問題．

1.3 非支配集合

僅使用單個(gè)能量函數(shù)對(duì)值進(jìn)行比較時(shí)，簡(jiǎn)單的大小判斷便能決定出優(yōu)劣．但是當(dāng)一個(gè)構(gòu)象由多個(gè)目標(biāo)函數(shù)共同評(píng)估時(shí)，將會(huì)存在支配與非支配關(guān)系[27]．

如圖2，在二維空間中進(jìn)行闡述．圖上橫縱坐標(biāo)分別為兩個(gè)不同的能量函數(shù)值，所有點(diǎn)代表解構(gòu)象在能量空間中的投影．解的關(guān)系分為兩類：支配與非支配．如果兩個(gè)解滿足

(4)

即解xA的第一個(gè)能量值與第二個(gè)能量值均比解xB對(duì)應(yīng)的能量值大，則稱解xB支配于解xA；如果兩個(gè)解滿足:

(5)

即解xA的第一個(gè)能量值比解xB大，解xA的第二個(gè)能量值比解xB小，則解xA與解xB為非支配關(guān)系．圖2中，解xA支配解xD、xE，解xD和解xE為非支配關(guān)系．

1.4 基于多目標(biāo)函數(shù)的構(gòu)象優(yōu)化

針對(duì)蛋白質(zhì)結(jié)構(gòu)的預(yù)測(cè)問題，結(jié)合非支配排序遺傳算法(non-dominated sorting genetic algorithm II，NSGA-Ⅱ)[28]思想提出TRIOFOLD方法．NSGA-II算法是進(jìn)化算法的一種，解決了在多個(gè)目標(biāo)優(yōu)化下解的排序的問題，具有較低的算法復(fù)雜度，TRIOFOLD流程如圖1．算法通過大量的碎片替換迭代，盡可能多地生成不同的構(gòu)象，在其中篩選出最優(yōu)解．由于迭代次數(shù)龐大，傳統(tǒng)的排序算法、整體時(shí)間復(fù)雜度很高．當(dāng)前的一般多目標(biāo)算法在選取非支配解時(shí)，由于只維護(hù)了非支配解，被支配的解完全被舍棄掉，使一些被支配的解無法發(fā)揮作用，從而損失了部分信息．而文中所提出的TRIOFOLD方法利用NSGA-II算法再次在被支配的解集中選取可利用的解來對(duì)其進(jìn)行進(jìn)一步優(yōu)化并將優(yōu)化后的結(jié)果重新評(píng)定，有助于從被支配的解集合中選擇保留部分可靠的解．圖3為多目標(biāo)函數(shù)優(yōu)化的流程．

圖3 多目標(biāo)函數(shù)優(yōu)化流程Fig.3 Multi-objective function optimization process

1.4.1 非支配排序

為了從構(gòu)象集合中篩選出最優(yōu)構(gòu)象集合，算法根據(jù)構(gòu)象的多個(gè)能量函數(shù)指標(biāo)進(jìn)行分層篩選，分別得到不同層級(jí)非支配前沿的構(gòu)象．

圖4記錄了每一個(gè)構(gòu)象被支配數(shù)與自身支配的構(gòu)象集合，每一個(gè)圓點(diǎn)代表一個(gè)構(gòu)象，箭頭代表支配關(guān)系．

圖4 非支配排序示例Fig.4 Example of non-dominatedsorting

1.4.2 精英保留策略與擁擠距離

(6)

1.4.3 迭代初始解的選取

由于更新操作由單一構(gòu)象進(jìn)行，在多輪迭代時(shí)，需要在當(dāng)前的最優(yōu)構(gòu)象集合P中選取用于下一輪迭代的起始構(gòu)象．基于單目標(biāo)能量函數(shù)構(gòu)象選取的理念，首先分層計(jì)算集合P中構(gòu)象的擁擠距離，選取每層具有最大擁擠距離的構(gòu)象作為迭代初始構(gòu)象，能夠生成更多樣化的結(jié)果．然后采用歸一化的方式計(jì)算選取每一層的最優(yōu)解的概率為：

(7)

式中：j為當(dāng)前層rank等級(jí)，其取值范圍為rank等級(jí)分配數(shù)量；u為最大rank等級(jí)；α為rank等級(jí)對(duì)應(yīng)權(quán)重．最后基于伯努利分布生成一個(gè)隨機(jī)變量，選擇隨機(jī)變量值對(duì)應(yīng)的區(qū)間構(gòu)象作為下一輪的迭代初始構(gòu)象．為了保證漸近收斂，算法在模擬退火的框架下進(jìn)行，參數(shù)α是模擬退火溫度系數(shù)T的負(fù)相關(guān)表達(dá)，隨著迭代次數(shù)的增加，溫度系數(shù)T越低，rank等級(jí)越高的α系數(shù)將越大，選擇到的概率越低，最終的構(gòu)象選取將逐漸向rank=1收斂，算法總是選擇在能量函數(shù)評(píng)價(jià)下最好的構(gòu)象．

在所有迭代結(jié)束后，對(duì)最終的最優(yōu)解集P，選出rank=1的非支配集合輸出，采用KNEE[29]算法從中挑選出最好的5個(gè)解作為最終的蛋白質(zhì)預(yù)測(cè)結(jié)構(gòu)．

1.4.4 算法整體流程

算法1 TRIOFOLD多目標(biāo)優(yōu)化流程輸入:初始線性構(gòu)象x'1,迭代次數(shù)s、單輪迭代碎片替換次數(shù)k;輸出:第s輪迭代結(jié)束后最優(yōu)構(gòu)象集合Ps中rank=1的構(gòu)象;1:初始化當(dāng)前迭代輪數(shù)i=1,構(gòu)象集合Xi、Pi-1、Pi、Ui均為空集;2:連續(xù)進(jìn)行k次碎片替換操作,得到構(gòu)象集合Xi={xi1,xi2,xi3,…,xik};3:合并構(gòu)象集合Xi與最優(yōu)構(gòu)象集合Pi-1,得到臨時(shí)構(gòu)象集合Ui;4:按照rank等級(jí)對(duì)集合Ui進(jìn)行非支配排序;5:構(gòu)象按rank從小到大依次填入Pi中直至到達(dá)數(shù)量上限,記集合Ui中剩余構(gòu)象最低rank等級(jí)為n;6:若n>4,舍棄Pi中所有rank>4的構(gòu)象;7:若n=4,取出Pi與Ui中rank=4的所有構(gòu)象,基于式(6)計(jì)算構(gòu)象擁擠距離,從小到大填入Pi中,舍棄多余構(gòu)象;8:基于式(6)計(jì)算出Pi中每一個(gè)rank等級(jí)的最優(yōu)構(gòu)象xbest;9:基于式(7)計(jì)算xbest出的選擇概率Pbest,按概率選出i+1輪迭代起始構(gòu)象xi+11;10:當(dāng)i?s時(shí),i=i+1,重復(fù)步驟2～9;當(dāng)i>s時(shí),進(jìn)入步驟11;11:輸出最優(yōu)構(gòu)象集合中所有rank=1的構(gòu)象.

2 實(shí)驗(yàn)

為評(píng)估TRIOFOLD算法的效果，實(shí)驗(yàn)基于ROSETTA的碎片拼接框架，所有代碼使用C++實(shí)現(xiàn)，在40個(gè)核256G內(nèi)存的Linux服務(wù)器上完成實(shí)驗(yàn)．實(shí)驗(yàn)數(shù)據(jù)選取第13屆蛋白質(zhì)結(jié)構(gòu)預(yù)測(cè)大賽(CASP13)上的10條序列作為初始的數(shù)據(jù)集，它們的序列長(zhǎng)度小于150，更便于集中驗(yàn)證方法的可行性．當(dāng)序列長(zhǎng)度大于150以后，隨著長(zhǎng)度的增加，構(gòu)象空間的大小將呈指數(shù)級(jí)增長(zhǎng)，基于短碎片的方法難以在短時(shí)間內(nèi)有效遍歷所有的可行解．

在所提出的TRIOFOLD方法實(shí)驗(yàn)中，一次運(yùn)行的迭代次數(shù)設(shè)置為1 250 次，每一次迭代將生成的構(gòu)象集合X與上一輪構(gòu)象最大數(shù)量為100的最優(yōu)集合P合并，按rank順序依次遞增排序后生成新的最優(yōu)集合，并挑選出最優(yōu)解作為下一輪迭代初始解．單輪程序執(zhí)行總計(jì)構(gòu)象更新次數(shù)為250 000次．實(shí)驗(yàn)結(jié)果為最優(yōu)集合rank=1的非支配前沿構(gòu)象，考慮到計(jì)算量大且耗時(shí)，實(shí)驗(yàn)采用多線程計(jì)算的方式，最終的結(jié)果以每個(gè)線程的rank=1的非支配前沿構(gòu)象直接相加得到．

實(shí)驗(yàn)結(jié)果與ROSETTA的原始方法進(jìn)行比較，兩種方法均基于傳統(tǒng)的碎片拼接思想進(jìn)行模擬，著重比較多個(gè)能量函數(shù)優(yōu)化與單個(gè)能量函數(shù)優(yōu)化間結(jié)果的差異．兩種方法分別運(yùn)行150次．大量的計(jì)算迭代保證了最終預(yù)測(cè)結(jié)構(gòu)盡可能接近真實(shí)蛋白結(jié)構(gòu)．同時(shí)相同的計(jì)算量保證算法的可比對(duì)性．由于ROSETTA僅有單個(gè)能量函數(shù)，ROSETTA方法每一輪程序計(jì)算只有1個(gè)結(jié)構(gòu)生成，最終輸出150個(gè)結(jié)構(gòu)，TRIOFOLD 方法每一輪輸出rank=1的所有結(jié)構(gòu)，最后對(duì)150 輪rank=1的集合合并，計(jì)算最終的非支配解集．根據(jù)具體蛋白質(zhì)的不同，最終生成的非支配解結(jié)構(gòu)數(shù)量在20～60不等，并用KNEE算法進(jìn)行解的排序．

表1展示了TRIOFOLD和ROSETTA在10個(gè)目標(biāo)序列上的預(yù)測(cè)結(jié)果，給出了預(yù)測(cè)結(jié)構(gòu)和天然結(jié)構(gòu)之間的全局距離測(cè)試分?jǐn)?shù)(global distance-total score，GDT-TS)．GDT-TS描述了預(yù)測(cè)結(jié)構(gòu)與真實(shí)結(jié)構(gòu)的相近程度，值越大越接近于真實(shí)結(jié)構(gòu)，選出每一個(gè)目標(biāo)蛋白GDT-TS評(píng)估前5的結(jié)構(gòu)作為比較．表1中可以看出TRIOFOLD整體性能優(yōu)于ROSETTA，以最優(yōu)結(jié)構(gòu)為例，在7個(gè)蛋白質(zhì)上，TRIOFOLD取得更優(yōu)的結(jié)果．

表1 CASP13數(shù)據(jù)集上兩種方法的預(yù)測(cè)結(jié)果

圖5為蛋白質(zhì)T0970進(jìn)行一輪快速非支配排序后前4個(gè)rank的構(gòu)象在能量空間的分布情況，非支配排序能夠有效地將生成的構(gòu)象進(jìn)行分層，并且在空間中均勻分布．

圖5 非支配排序后構(gòu)象在能量空間中的分布Fig.5 Conformation distribution in energy space after non-dominated sorting

3個(gè)能量函數(shù)的優(yōu)化下構(gòu)象之間間距較大，體現(xiàn)了構(gòu)象的多樣化．當(dāng)前迭代的構(gòu)象集合將按照rank從低到高的順序填入最優(yōu)集合．

圖6選取蛋白質(zhì)T0970，展示了兩種方法生成的最終構(gòu)象集合在3個(gè)能量函數(shù)中的空間分布圖，圖中不同的點(diǎn)代表了構(gòu)象在能量空間的坐標(biāo)．z軸的CHARMM能量函數(shù)由于結(jié)構(gòu)間差異較大，在比較時(shí)做了標(biāo)準(zhǔn)化處理，減小離群點(diǎn)影響的同時(shí)保證能量的相對(duì)大小關(guān)系．在相同的計(jì)算量下，TRIOFOLD方法為篩選出最優(yōu)解提供的備選解集更小，更容易獲得最優(yōu)解．對(duì)于RWPLUS和CHARMM能量函數(shù)，ROSETTA方法未進(jìn)行優(yōu)化，整體的能量值相比TRIOFOLD更高，生成的結(jié)構(gòu)是最優(yōu)解的概率更低．

圖6 TRIOFOLD和ROSETTA生成構(gòu)象分布Fig.6 Conformational distribution of TRIOFOLD and ROSETTA

圖7為TRIOFOLD和ROSETTA兩種方法在蛋白質(zhì)T0958上最好的5個(gè)模型的預(yù)測(cè)結(jié)果．圖7中上方為TRIOFOLD的預(yù)測(cè)結(jié)構(gòu)，下方為ROSETTA預(yù)測(cè)結(jié)構(gòu)與原始結(jié)構(gòu)的對(duì)比．TRIOFOLD預(yù)測(cè)結(jié)構(gòu)的平均GDT-TS得分為54.16，均方根偏差(root-mean-square deviation,RMSD)為6.983，對(duì)應(yīng)的ROSETTA能量為-169．ROSETTA 方法預(yù)測(cè)結(jié)構(gòu)的平均GDT-TS得分為45.85，RMSD為7.301，對(duì)應(yīng)的ROSETTA能量值為-179．其原因?yàn)橛捎诓煌哪芰亢瘮?shù)對(duì)不同的結(jié)構(gòu)傾向性不同，多目標(biāo)函數(shù)優(yōu)化使得TRIOFOLD在二級(jí)結(jié)構(gòu)的空間預(yù)測(cè)上更為合理，其二級(jí)結(jié)構(gòu)的預(yù)測(cè)均與原始結(jié)構(gòu)更加相近．ROSETTA的預(yù)測(cè)在α螺旋上存在著更大的扭轉(zhuǎn)角偏差，部分結(jié)構(gòu)上長(zhǎng)度差異也較大．β折疊與無規(guī)則卷曲部分與原始結(jié)構(gòu)在空間位置上有一定差異．表2為對(duì)應(yīng)的統(tǒng)計(jì)結(jié)果，其中RMSD是用于比較兩個(gè)蛋白質(zhì)結(jié)構(gòu)的相似度，其計(jì)算了兩個(gè)蛋白質(zhì)對(duì)應(yīng)CA原子之間的均方根偏差．RMSD越小，代表兩個(gè)蛋白質(zhì)結(jié)構(gòu)越相似．由表2看出，TRIOFOLD方法在GDT-TS和RMSD上均更優(yōu)于ROSETTA，但是單個(gè)能量函數(shù)的能量值略微偏高，還需要進(jìn)一步的優(yōu)化算法．

圖7 TRIOFOLD和ROSETTA的預(yù)測(cè)結(jié)構(gòu) 與真實(shí)結(jié)構(gòu)的比對(duì)Fig.7 Comparison of predicted structure and real structure of TRIOFOLD and ROSETTA

表2 蛋白質(zhì)T0958的預(yù)測(cè)結(jié)構(gòu)能量值

3 結(jié)論

文中提出了基于多目標(biāo)優(yōu)化的TRIOFOLD方法，并通過實(shí)驗(yàn)驗(yàn)證了該方法的可行性．TRIOFOLD采用3個(gè)能量函數(shù)評(píng)估蛋白質(zhì)結(jié)構(gòu)，能夠提供更魯棒的優(yōu)化方向．然而由于能量函數(shù)本身存在著一些誤差，3個(gè)能量函數(shù)雖然能夠?yàn)樯傻牡鞍踪|(zhì)結(jié)構(gòu)提供良好的優(yōu)化方向，但是由于碎片替換的機(jī)制，對(duì)于優(yōu)化的步長(zhǎng)卻很難定量的把握．這導(dǎo)致TRIOFOLD生成的結(jié)構(gòu)能獲得有效的收斂，但是與該蛋白的原始結(jié)構(gòu)還是存在著一定的差距．對(duì)于一個(gè)蛋白質(zhì)而言，不同的能量函數(shù)對(duì)它的靈敏度也不相同，實(shí)驗(yàn)中往往存在一個(gè)能量函數(shù)已經(jīng)收斂而其他能量函數(shù)依舊處于可優(yōu)化狀態(tài)的現(xiàn)象，如何有效的為不同能量函數(shù)分配不同的優(yōu)化方式也是待研究的問題．另外GDT-TS分?jǐn)?shù)也存在著一定的與能量函數(shù)值非正相關(guān)的情況，表2中的數(shù)據(jù)有一定的體現(xiàn)，足夠低的能量評(píng)估值僅是生成更好結(jié)構(gòu)的必要條件．TRIOFOLD將直接采用GDT-TS分?jǐn)?shù)作為優(yōu)化指標(biāo)的方式進(jìn)行構(gòu)象的迭代更新，將其與3個(gè)能量函數(shù)值并列加入的優(yōu)化過程中，而不是僅作為最后評(píng)估的指標(biāo)．

另一方面，生成整體水平更優(yōu)的備選解集是下一步重要的研究方向．受限于當(dāng)前算法基于碎片拼接的限制，對(duì)于一些具有很少同源結(jié)構(gòu)的蛋白，TRIOFOLD難以搜索到優(yōu)質(zhì)的碎片結(jié)構(gòu)，導(dǎo)致最終的解結(jié)構(gòu)與真實(shí)結(jié)構(gòu)相差甚遠(yuǎn)．對(duì)此，TRIOFOLD將探索加入更多的約束條件，例如序列的接觸圖信息等特征，為難以預(yù)測(cè)的蛋白質(zhì)結(jié)構(gòu)提供更多的計(jì)算信息，這樣在已有的先驗(yàn)條件下，將縮小優(yōu)化的搜索空間，并且已知的數(shù)據(jù)能夠提供更高的準(zhǔn)確性．同時(shí)TRIOFOLD也將探索采用深度學(xué)習(xí)的方式，通過統(tǒng)計(jì)蛋白質(zhì)數(shù)據(jù)庫(kù)中已有的17萬個(gè)蛋白的結(jié)構(gòu)特征，為序列找出更多的可行同源結(jié)構(gòu)和訓(xùn)練出更準(zhǔn)確的能量函數(shù)，從而優(yōu)化現(xiàn)在的算法框架．

綜上所述，文中提出的基于多目標(biāo)函數(shù)并行優(yōu)化的蛋白質(zhì)三維結(jié)構(gòu)預(yù)測(cè)方法，相對(duì)于傳統(tǒng)的基于單目標(biāo)能量函數(shù)的優(yōu)化方法，TRIOFOLD充分利用了3個(gè)能量函數(shù)的勢(shì)能面差異，結(jié)合進(jìn)化算法NSGA-II的思想進(jìn)行構(gòu)象集合的迭代更新與篩選，在CASP13數(shù)據(jù)集上取得了較好的效果．是領(lǐng)域內(nèi)已有蛋白質(zhì)三維結(jié)構(gòu)算法的有益補(bǔ)充．未來的工作將側(cè)重于解決能量函數(shù)的精度問題和有效的組合使用問題，以及提高TRIOFOLD對(duì)于長(zhǎng)序列蛋白和具有較低同源性蛋白的預(yù)測(cè)能力．