接觸圖輔助的過程重采樣蛋白質構象空間優化算法

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

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

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

E-mail：zgj@zjut.edu.cn

1 引 言

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

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

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

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

2 算法設計

2.1 過程片段庫構建

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

圖1 構建過程片段庫

2.2 殘基接觸打分模型

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

(1)

(2)

根據設計的殘基接觸打分模型，如果一個蛋白質構象的殘基接觸得分較低，換言之，這個蛋白質構象的殘基接觸和預測的殘基接觸很符合，那么就稱這個蛋白質構象的整體結構是相對合理的.

2.3 過程重采樣策略

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

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

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

(3)

kT為溫度常數；

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

(4)

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

若連續拒絕插入次數達到最大連續拒絕次數Rmax，令kT=kT+s，s表示升溫步長.

2.4 交叉操作

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

圖2 交叉操作

2.5 算法描述：

算法1.CMPR算法流程

輸入：目標蛋白序列，片段庫，殘基接觸圖

輸出：蛋白質三維結構

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

3 實驗結果與分析

3.1 測試蛋白和參數設置

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

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

表1 測試蛋白集

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

實驗采用均方根偏差(RMSD)和TM-score[47]兩個評價指標比較目標蛋白的預測結構和天然態蛋白的結構相似度，RMSD值越小表明預測結構與天然態蛋白質結構越相似，TM-score的值在[0，1]之間，TM-score值越大表明預測結構與天然態蛋白質結構越相似.

3.2 結果分析

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

表2 算法對比結果

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

表3 CMPR組件比較

對比表2和表3的結果可以看出，僅使用殘基接觸的CM算法，所有12個測試蛋白的TM-sore比Rosetta高，平均TM-score值比Rosetta高0.11，平均RMSD比Rosetta低1.68?.僅使用過程重采樣的PR算法，12個測試蛋白中9個測試蛋白的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?.實驗結果表明，僅使用殘基接觸信息和過程重采樣方法相比較于Rosetta都提高了蛋白質結構預測精度，而CMPR結合殘基接觸和過程重采樣策略得到了更高的預測精度.

3.3 算法采樣能力分析

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

圖3 構象空間采樣比較

圖4 構象RMSD分布

圖4所示橫坐標表示構象與天然態蛋白質結構之間的RMSD，縱坐標表示每個區域內構象數目所占構象總數的百分比.從圖中4個測試蛋白的構象RMSD分布圖可以看出，CMPR采樣到的低RMSD構象占總構象數目的比例比Rosetta采樣的更大，因此CMPR算法比Rosetta有更好的近天然態采樣能力.

4 結 論

本文提出了一種接觸圖輔助的過程重采樣蛋白質構象空間優化算法，在進化計算的框架下，基于殘基接觸圖設計了殘基接觸的打分模型，使用殘基接觸打分模型篩選種群，提取結構合理構象的片段信息構建過程片段庫，結合殘基接觸圖的先驗知識和種群進化的過程信息指導構象空間采樣.12個測試蛋白的實驗結果表明，本文提出的CMPR算法有效地縮小了構象搜索空間，增強了對近天然態構象區域的搜索，有較高的預測精度，是一種有效的蛋白質結構從頭預測方法.下一步的研究中，將結合更多的先驗知識構建更精準的過程片段庫進行蛋白質結構從頭預測.