李章維，余寶昆，胡 俊，周曉根，2，張貴軍

1(浙江工業(yè)大學(xué) 信息工程學(xué)院，杭州 310023)

2(密西根大學(xué) 計(jì)算醫(yī)藥與生物信息學(xué)院，美國密西根州安娜堡 48109)

E-mail：zgj@zjut.edu.cn

1 引 言

作為人體內(nèi)分布最廣、功能最復(fù)雜的一類生物大分子，蛋白質(zhì)在人體的生命活動(dòng)中發(fā)揮著極其重要的作用.研究蛋白質(zhì)的結(jié)構(gòu)是研究其生物功能及活性機(jī)理的基礎(chǔ)，并且對新蛋白的發(fā)明和藥物標(biāo)靶蛋白的設(shè)計(jì)具有十分重要的指導(dǎo)意義[1].傳統(tǒng)實(shí)驗(yàn)測定蛋白質(zhì)結(jié)構(gòu)的方法(如：X射線晶體衍射法、核磁共振成像與冷凍電鏡等)所需時(shí)間長且花費(fèi)巨大，無法滿足大規(guī)模蛋白質(zhì)結(jié)構(gòu)測定的需求[2,3].因此，從蛋白質(zhì)的一維氨基酸序列出發(fā)，直接預(yù)測蛋白質(zhì)三維結(jié)構(gòu)是結(jié)構(gòu)生物學(xué)領(lǐng)域的一個(gè)重要問題[4].

自由模板(FM)是蛋白質(zhì)結(jié)構(gòu)預(yù)測領(lǐng)域的最具挑戰(zhàn)性的一類研究課題.對于任意一條蛋白質(zhì)的氨基酸序列，其可能形成的空間結(jié)構(gòu)數(shù)目是極其巨大的，現(xiàn)有的計(jì)算資源難以在巨大的構(gòu)象空間中搜索到能量最低的構(gòu)象[5].針對蛋白質(zhì)構(gòu)象空間優(yōu)化問題，研究者們提出了很多構(gòu)象空間優(yōu)化算法，如進(jìn)化算法[6-10]、蒙特卡洛算法[11-13]、副本交換算法[14-16]、構(gòu)象空間退火[17]等方法.為了進(jìn)一步減小構(gòu)象搜索空間，基于片段組裝策略[18]的構(gòu)象空間搜索在從頭預(yù)測中得到廣泛應(yīng)用，其主要代表有Rosetta[19-21]，Quark[22]，STRATCH[23]，PROFESY[24]，F(xiàn)RAGFOLD[25]等一系列方法，特別是Rosetta和Quark方法在國際蛋白質(zhì)結(jié)構(gòu)預(yù)測競賽(CASP)[26，27]中表現(xiàn)突出.

為了提高蛋白質(zhì)結(jié)構(gòu)預(yù)測精度，CONFOLD[28,29]使用預(yù)測的二級結(jié)構(gòu)和殘基接觸[30]轉(zhuǎn)化成空間約束，然后使用這些空間約束構(gòu)建蛋白三維結(jié)構(gòu)模型.Filb-Coevo[31]使用殘基接觸圖約束產(chǎn)生高質(zhì)量的片段庫[32，33]，進(jìn)而使用片段組裝方法搜索構(gòu)象.RMA[34]算法在遺傳算法的框架下使用預(yù)測的二級結(jié)構(gòu)增強(qiáng)對構(gòu)象采樣空間的探索.DPDE[35]算法使用距離譜[36]指導(dǎo)差分進(jìn)化進(jìn)行蛋白質(zhì)結(jié)構(gòu)預(yù)測.SCDE[37]算法使用基于二級結(jié)構(gòu)和殘基接觸的選擇策略指導(dǎo)構(gòu)象空間采樣.在進(jìn)化計(jì)算框架下[38-40]，RMA[34]、DPDE[35]、SCDE[37]算法使用蛋白質(zhì)結(jié)構(gòu)的先驗(yàn)知識輔助蛋白質(zhì)結(jié)構(gòu)預(yù)測，有效地提升了蛋白質(zhì)結(jié)構(gòu)預(yù)測精度.除了蛋白質(zhì)結(jié)構(gòu)的先驗(yàn)知識，種群更新中的過程信息同樣十分重要，由于能量函數(shù)的不精確，在種群進(jìn)化的過程中，一些結(jié)構(gòu)合理但能量較高的構(gòu)象可能會(huì)被淘汰.

為了保留結(jié)構(gòu)合理構(gòu)象的片段信息達(dá)到指導(dǎo)種群進(jìn)化目的，提高蛋白質(zhì)結(jié)構(gòu)預(yù)測精度，本文提出了一種接觸圖輔助的過程重采樣蛋白質(zhì)構(gòu)象空間優(yōu)化算法(CMPR).在種群更新的每一代，首先，根據(jù)殘基接觸打分模型選擇結(jié)構(gòu)合理的構(gòu)象，然后根據(jù)選出的構(gòu)象構(gòu)建過程片段庫，基于過程片段庫使用過程重采樣策略指導(dǎo)種群變異.12個(gè)測試蛋白的實(shí)驗(yàn)結(jié)果表明，CMPR方法縮小了構(gòu)象搜索空間，增強(qiáng)了對近天然態(tài)構(gòu)象區(qū)域的采樣，相對于Rosetta和Quark方法提升了蛋白質(zhì)結(jié)構(gòu)預(yù)測的精度.

2 算法設(shè)計(jì)

2.1 過程片段庫構(gòu)建

在采樣蛋白質(zhì)構(gòu)象的過程中，從采樣過程構(gòu)象中提取出來的片段庫稱之為過程片段庫.圖1所示是目標(biāo)序列第一個(gè)位置的3～9殘基過程片段庫的構(gòu)建.對于序列長度為L的N個(gè)蛋白質(zhì)構(gòu)象，記錄N個(gè)構(gòu)象第一個(gè)3～9殘基窗口內(nèi)的二面角(φ，φ)構(gòu)成3～9殘基片段庫，然后窗口向下滑動(dòng)記錄第二個(gè)位置的3～9殘基片段庫，直到構(gòu)建目標(biāo)序列所有位置上的過程片段庫.

圖1 構(gòu)建過程片段庫

2.2 殘基接觸打分模型

當(dāng)?shù)鞍踪|(zhì)序列上兩個(gè)殘基間Cβ-Cβ(甘氨酸Cα-Cα)的空間距離小于8?，就認(rèn)為這兩個(gè)殘基之間有殘基接觸[41].為了在種群更新的過程中提取構(gòu)象的片段信息，本文設(shè)計(jì)了殘基接觸的打分模型來選擇構(gòu)象構(gòu)建過程片段庫.根據(jù)目標(biāo)蛋白序列，使用RaptorX-Contact[42，43]預(yù)測目標(biāo)蛋白的殘基接觸圖，使用如下公式計(jì)算蛋白質(zhì)構(gòu)象的殘基接觸得分：

(1)

(2)

根據(jù)設(shè)計(jì)的殘基接觸打分模型，如果一個(gè)蛋白質(zhì)構(gòu)象的殘基接觸得分較低，換言之，這個(gè)蛋白質(zhì)構(gòu)象的殘基接觸和預(yù)測的殘基接觸很符合，那么就稱這個(gè)蛋白質(zhì)構(gòu)象的整體結(jié)構(gòu)是相對合理的.

2.3 過程重采樣策略

為了使用過程片段庫指導(dǎo)種群變異，本文設(shè)計(jì)了過程片段庫的片段組裝策略.對一個(gè)構(gòu)象進(jìn)行片段組裝之前，首先使用Rosetta的能量函數(shù)score3[19]計(jì)算構(gòu)象的能量E，公式(1)和公式(2)計(jì)算構(gòu)象的殘基接觸分?jǐn)?shù)Scon，然后在構(gòu)象上隨機(jī)選擇一個(gè)長度為[3，9]殘基長度的片段插入窗口，從窗口對應(yīng)位置和片段長度的過程片段庫中隨機(jī)選擇一個(gè)片段，使用該片段中的二面角φ和φ替換當(dāng)前構(gòu)象中的二面角.當(dāng)片段替換完成后，重新計(jì)算片段插入后構(gòu)象的能量E′和殘基接觸分?jǐn)?shù)S′con，根據(jù)如下條件判斷是否接收片段插入：

a)若E′-E≤0，S′con-Scon≤0，則接收這次片段插入；

b)若E′-E≤0，S′con-Scon>0，則根據(jù)概率Pcon判斷是否接收，Pcon公式如下：

(3)

kT為溫度常數(shù)；

c)若E′-E>0，S′con-Scon≤0，則根據(jù)概率Pe判斷是否接收，Pe公式如下：

(4)

d)若E′-E>0，S′con-Scon>0，則拒絕這次片段插入.

若連續(xù)拒絕插入次數(shù)達(dá)到最大連續(xù)拒絕次數(shù)Rmax，令kT=kT+s，s表示升溫步長.

2.4 交叉操作

圖2所示為個(gè)體xi交叉操作示意圖，對于個(gè)體xi，首先從種群中隨機(jī)選擇一個(gè)個(gè)體xj，xi≠xj，然后在xj上隨機(jī)選擇一個(gè)3殘基長度的片段，使用這個(gè)3片段上的二面角替換xi相同位置片段上的二面角，得到交叉后的個(gè)體x′i，交叉操作增加了種群的多樣性.

圖2 交叉操作

2.5 算法描述：

算法1.CMPR算法流程

輸入：目標(biāo)蛋白序列，片段庫，殘基接觸圖

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

其中NP表示種群規(guī)模，CR表示交叉選擇概率，p表示變異選擇概率，Q表示片段庫提取比例.第1行表示種群初始化過程，對NP個(gè)全伸展的構(gòu)象使用Rosetta的一階段和二階段進(jìn)行初始化.第24行種群選擇操作中，把交叉變異前的父代種群和交叉變異后的子代種群合并，使用殘基接觸打分模型對合并種群個(gè)體打分，保留殘基接觸得分最低的前一半個(gè)體進(jìn)入下一代種群.第25行根據(jù)殘基接觸打分模型從NP個(gè)個(gè)體中選擇殘基接觸分?jǐn)?shù)最低的Q*NP個(gè)體構(gòu)建過程片段庫.

3 實(shí)驗(yàn)結(jié)果與分析

3.1 測試蛋白和參數(shù)設(shè)置

實(shí)驗(yàn)選擇CASP12競賽中的12個(gè)測試蛋白，如表1所示，這些測試蛋白序列長度從39到146，折疊類型包括α、β和α/β類型，所有測試蛋白從Robetta[44]服務(wù)器下載去除同源的片段庫.本算法在Rosetta平臺(tái)上實(shí)現(xiàn)，使用Rosetta3.4版本的ClassicAbinitio協(xié)議.

設(shè)置種群規(guī)模NP=100，種群更新次數(shù)Gmax=100，交叉概率CR=0.5，變異概率p=0.5.最大連續(xù)拒絕次數(shù)Rmax=150，溫度參數(shù)初始值kT=2，升溫步長s=2，片段庫提取比例Q=0.1.通常研究者會(huì)在L，L/2，L/5，L/10之間選擇使用的殘基接觸數(shù)量[45]，L為目標(biāo)蛋白的序列長度，殘基接觸數(shù)據(jù)根據(jù)置信度排序，置信度越高排序越靠前，本實(shí)驗(yàn)選擇前L個(gè)殘基接觸數(shù)據(jù)計(jì)算蛋白質(zhì)構(gòu)象的殘基接觸分?jǐn)?shù).

表1 測試蛋白集

在CMPR算法運(yùn)行的過程中，記錄除種群初始化階段之外構(gòu)象搜索過程中接收的所有構(gòu)象，使用聚類工具SPICKER[46]對這些構(gòu)象聚類，記錄聚類的第一個(gè)類心點(diǎn)作為預(yù)測結(jié)果.Rosetta方法使用Rosetta默認(rèn)的參數(shù)，保持和CMPR相同的計(jì)算代價(jià)，記錄第三第四階段蒙特卡洛軌跡接收的所有構(gòu)象經(jīng)過聚類得到的第一個(gè)聚類中心作為預(yù)測結(jié)果.Quark的預(yù)測結(jié)果通過Quark服務(wù)器在線提交任務(wù)得到5個(gè)預(yù)測蛋白結(jié)構(gòu)模型，選擇第一個(gè)模型作為預(yù)測結(jié)果和其他兩種方法比較.

實(shí)驗(yàn)采用均方根偏差(RMSD)和TM-score[47]兩個(gè)評價(jià)指標(biāo)比較目標(biāo)蛋白的預(yù)測結(jié)構(gòu)和天然態(tài)蛋白的結(jié)構(gòu)相似度，RMSD值越小表明預(yù)測結(jié)構(gòu)與天然態(tài)蛋白質(zhì)結(jié)構(gòu)越相似，TM-score的值在[0，1]之間，TM-score值越大表明預(yù)測結(jié)構(gòu)與天然態(tài)蛋白質(zhì)結(jié)構(gòu)越相似.

3.2 結(jié)果分析

表2中是CMPR算法與Rosetta和Quark的預(yù)測結(jié)果的比較.在12個(gè)測試蛋白中，CMPR預(yù)測結(jié)果的TM-score大于0.5的有8個(gè)，其中有11個(gè)測試蛋白的預(yù)測精度在TM-score和RMSD上優(yōu)于Rosetta，有7個(gè)測試蛋白在TM-score和RMSD上都優(yōu)于Quark.其中測試蛋白1GYZ、1AIL、1GB1、1SAP、1TIG、1BQ9、1WAP在TM-score和RMSD上，CMPR的預(yù)測精度比其他兩種方法高.CMPR算法的平均TM-score比Rosetta高0.16，比Quark高0.04，平均RMSD比Rosetta低3.13?，比Quark低1.51?.總體而言，CMPR算法比Rosetta和Quark預(yù)測精度更高.

表2 算法對比結(jié)果

為了驗(yàn)證殘基接觸信息的使用和過程重采樣策略的有效性，本文設(shè)計(jì)了僅使用殘基接觸信息和僅使用過程重采樣策略的對比實(shí)驗(yàn).CM算法把CMPR算法中的變異操作更改為僅使用原片段庫進(jìn)行片段組裝.PR算法把CMPR算法中的種群篩選中殘基接觸打分模型替換成Rosetta的score3，使用score3選擇低能量構(gòu)象構(gòu)建過程片段庫，過程片段庫的片段組裝中，若E′-E>0，使用公式(4)判斷是否接收片段插入，否則，接收片段插入.對比實(shí)驗(yàn)結(jié)果如表3所示.

表3 CMPR組件比較

對比表2和表3的結(jié)果可以看出，僅使用殘基接觸的CM算法，所有12個(gè)測試蛋白的TM-sore比Rosetta高，平均TM-score值比Rosetta高0.11，平均RMSD比Rosetta低1.68?.僅使用過程重采樣的PR算法，12個(gè)測試蛋白中9個(gè)測試蛋白的TM-score比Rosetta高，其平均TM-score值比Rosetta高0.06，其平均RMSD比Rosetta低0.91?.CMPR算法平均TM-score比CM算法高0.05，比PR算法高0.1，平均RMSD比CM算法低1.45?，比PR算法低2.22?.實(shí)驗(yàn)結(jié)果表明，僅使用殘基接觸信息和過程重采樣方法相比較于Rosetta都提高了蛋白質(zhì)結(jié)構(gòu)預(yù)測精度，而CMPR結(jié)合殘基接觸和過程重采樣策略得到了更高的預(yù)測精度.

3.3 算法采樣能力分析

部分蛋白的構(gòu)象采樣結(jié)果如圖3所示，圖中橫坐標(biāo)表示構(gòu)象與天然態(tài)蛋白質(zhì)比較的RMSD值，縱坐標(biāo)表示構(gòu)象的能量值.從圖中測試蛋白4UEX、1SAP可以看出，雖然CMPR算法總體的構(gòu)象采樣范圍比Rosetta小，但是比Rosetta搜索到更多能量更低和RMSD更低的構(gòu)象區(qū)域，因此預(yù)測精度更高.1BQ9和1WAP是兩個(gè)結(jié)構(gòu)較復(fù)雜的β折疊型蛋白，Rosetta雖然搜索到了更低能量的構(gòu)象區(qū)域，但是由于能量函數(shù)的不精確，低能構(gòu)象的結(jié)構(gòu)與天然態(tài)蛋白質(zhì)結(jié)構(gòu)相差較大，CMPR算法雖然無法采樣到比Rosetta更低能量的構(gòu)象區(qū)域，但是得到的最低能量區(qū)域的構(gòu)象RMSD更小，所以得到結(jié)果的預(yù)測精度更高.

圖3 構(gòu)象空間采樣比較

圖4 構(gòu)象RMSD分布

圖4所示橫坐標(biāo)表示構(gòu)象與天然態(tài)蛋白質(zhì)結(jié)構(gòu)之間的RMSD，縱坐標(biāo)表示每個(gè)區(qū)域內(nèi)構(gòu)象數(shù)目所占構(gòu)象總數(shù)的百分比.從圖中4個(gè)測試蛋白的構(gòu)象RMSD分布圖可以看出，CMPR采樣到的低RMSD構(gòu)象占總構(gòu)象數(shù)目的比例比Rosetta采樣的更大，因此CMPR算法比Rosetta有更好的近天然態(tài)采樣能力.

4 結(jié) 論

本文提出了一種接觸圖輔助的過程重采樣蛋白質(zhì)構(gòu)象空間優(yōu)化算法，在進(jìn)化計(jì)算的框架下，基于殘基接觸圖設(shè)計(jì)了殘基接觸的打分模型，使用殘基接觸打分模型篩選種群，提取結(jié)構(gòu)合理構(gòu)象的片段信息構(gòu)建過程片段庫，結(jié)合殘基接觸圖的先驗(yàn)知識和種群進(jìn)化的過程信息指導(dǎo)構(gòu)象空間采樣.12個(gè)測試蛋白的實(shí)驗(yàn)結(jié)果表明，本文提出的CMPR算法有效地縮小了構(gòu)象搜索空間，增強(qiáng)了對近天然態(tài)構(gòu)象區(qū)域的搜索，有較高的預(yù)測精度，是一種有效的蛋白質(zhì)結(jié)構(gòu)從頭預(yù)測方法.下一步的研究中，將結(jié)合更多的先驗(yàn)知識構(gòu)建更精準(zhǔn)的過程片段庫進(jìn)行蛋白質(zhì)結(jié)構(gòu)從頭預(yù)測.