張海倉 高玉娟 鄧明華,4,5 鄭偉謀 卜東波

1(中國(guó)科學(xué)院計(jì)算技術(shù)研究所 北京 100190)(中國(guó)科學(xué)院大學(xué) 北京 100049)(北京大學(xué)定量生物學(xué)中心 北京 100871)(北京大學(xué)數(shù)學(xué)科學(xué)學(xué)院 北京 100871)(北京大學(xué)統(tǒng)計(jì)科學(xué)中心 北京 100871) (中國(guó)科學(xué)院理論物理研究所 北京 100190)(zhanghaicang@ict.ac.cn)

蛋白質(zhì)中殘基遠(yuǎn)程相互作用預(yù)測(cè)算法研究綜述

張海倉1,2高玉娟3鄧明華3,4,5鄭偉謀6卜東波1

1(中國(guó)科學(xué)院計(jì)算技術(shù)研究所 北京 100190)(中國(guó)科學(xué)院大學(xué) 北京 100049)(北京大學(xué)定量生物學(xué)中心 北京 100871)(北京大學(xué)數(shù)學(xué)科學(xué)學(xué)院 北京 100871)(北京大學(xué)統(tǒng)計(jì)科學(xué)中心 北京 100871) (中國(guó)科學(xué)院理論物理研究所 北京 100190)(zhanghaicang@ict.ac.cn)

蛋白質(zhì)是由多個(gè)氨基酸殘基順序連接而成的長(zhǎng)鏈.在天然狀態(tài)下，蛋白質(zhì)并不是無規(guī)則的自由狀態(tài)，而是自發(fā)形成特定的空間結(jié)構(gòu)，以執(zhí)行其特定的生物學(xué)功能.驅(qū)動(dòng)蛋白質(zhì)形成特定空間結(jié)構(gòu)的主要因素是殘基間的非共價(jià)相互作用，包括疏水作用、靜電相互作用、范德華力等.因此，對(duì)殘基之間遠(yuǎn)程相互作用的準(zhǔn)確預(yù)測(cè)將有助于對(duì)蛋白質(zhì)空間結(jié)構(gòu)的預(yù)測(cè)，進(jìn)而有助于對(duì)蛋白質(zhì)生物學(xué)功能的了解.在蛋白質(zhì)進(jìn)化過程，有相互作用殘基對(duì)之間存在一種“共進(jìn)化”模式，即當(dāng)一個(gè)殘基發(fā)生變異時(shí)，與其有相互作用的殘基也要發(fā)生相應(yīng)的變異，以維持相互作用，進(jìn)而維持整體空間結(jié)構(gòu)以及生物學(xué)功能.基于上述生物學(xué)觀察，研究者開發(fā)了多個(gè)統(tǒng)計(jì)模型和算法以預(yù)測(cè)殘基對(duì)之間的相互作用:1)概述殘基之間遠(yuǎn)程相互作用的兩大類基本預(yù)測(cè)算法，包括無監(jiān)督學(xué)習(xí)方法和監(jiān)督學(xué)習(xí)方法；2)使用蛋白質(zhì)結(jié)構(gòu)預(yù)測(cè)CASP比賽結(jié)果來客觀比較上述各類算法的性能，分析各個(gè)算法的特點(diǎn)和優(yōu)勢(shì)；3)從生物學(xué)觀察和統(tǒng)計(jì)模型2個(gè)角度分析總結(jié)了未來的發(fā)展趨勢(shì).

殘基遠(yuǎn)程相互作用預(yù)測(cè)；蛋白質(zhì)三級(jí)結(jié)構(gòu)預(yù)測(cè)；圖模型；共進(jìn)化；機(jī)器學(xué)習(xí)

蛋白質(zhì)是生物體的重要組成成分，行使催化、免疫、細(xì)胞信號(hào)傳導(dǎo)等重要的生物學(xué)功能[1-2].

蛋白質(zhì)的基本組成單元是氨基酸，常見的氨基酸有20種.蛋白質(zhì)是以氨基酸為單元，脫水后由肽鍵連接而成的長(zhǎng)鏈(氨基酸脫水之后的剩余部分被稱為殘基).因此從計(jì)算的觀點(diǎn)看，可以將蛋白質(zhì)抽象表示成一個(gè)字符串序列，其字符集規(guī)模為20，其中每一個(gè)字符表示一種氨基酸殘基,如圖1所示:

Fig. 1 Illustration of amino acids,peptide bond,and protein sequence圖1 氨基酸、肽鍵以及蛋白質(zhì)序列示意圖

在天然環(huán)境下，蛋白質(zhì)呈現(xiàn)的并不是松散的、無規(guī)則的形態(tài)，而是自發(fā)折疊成特定的空間結(jié)構(gòu)，其中每個(gè)殘基(確切地說是殘基中的每個(gè)原子)都有其特定的空間坐標(biāo).蛋白質(zhì)的空間結(jié)構(gòu)決定了其生化功能，因此，認(rèn)識(shí)蛋白質(zhì)的空間結(jié)構(gòu)對(duì)了解其功能至關(guān)重要[2].

目前測(cè)定蛋白質(zhì)結(jié)構(gòu)的主要實(shí)驗(yàn)技術(shù)包括核磁共振[3]、X-ray晶體衍射[4]和冷凍電鏡[5]等.然而上述實(shí)驗(yàn)測(cè)定技術(shù)常常受限于蛋白質(zhì)大小、蛋白質(zhì)能否結(jié)晶以及結(jié)構(gòu)測(cè)定的高成本等因素，使得蛋白質(zhì)結(jié)構(gòu)的測(cè)定速度遠(yuǎn)遠(yuǎn)達(dá)不到蛋白質(zhì)序列的測(cè)定速度，因而通過計(jì)算的方法預(yù)測(cè)蛋白質(zhì)結(jié)構(gòu)具有重要的研究意義.另一方面，從序列出發(fā)進(jìn)行蛋白質(zhì)空間結(jié)構(gòu)預(yù)測(cè)是可行的[6]：Anfinsen的經(jīng)典實(shí)驗(yàn)表明在一般情況下，蛋白質(zhì)折疊是一個(gè)自發(fā)過程；或者換句話說，蛋白質(zhì)的結(jié)構(gòu)信息完全蘊(yùn)含于其序列之中，從而意味著蛋白質(zhì)結(jié)構(gòu)預(yù)測(cè)的可行性[7].

驅(qū)動(dòng)蛋白質(zhì)序列形成特定空間結(jié)構(gòu)的主要因素是殘基之間的大量非共價(jià)相互作用，包括疏水作用、范德華力(van der Waals forces)、離子鍵以及氫鍵等.從具有相互作用的殘基間序列距離來看，上述相互作用可以分作近程相互作用和遠(yuǎn)程相互作用2類，其中近程相互作用主導(dǎo)蛋白質(zhì)形成局部結(jié)構(gòu)，而遠(yuǎn)程相互作用則引導(dǎo)局部結(jié)構(gòu)的合理擺放，最終形成穩(wěn)定的蛋白質(zhì)空間結(jié)構(gòu)[8].相對(duì)于近程相互作用而言，遠(yuǎn)程相互作用具有決定整體結(jié)構(gòu)框架的重要作用，從而獲得了更多的關(guān)注和研究.本文著重討論殘基間遠(yuǎn)程相互作用預(yù)測(cè)問題.

蛋白質(zhì)中殘基間是否有相互作用可使用殘基之間的歐氏距離作判據(jù)，即有相互作用的殘基之間距離一般較小.通常采用的標(biāo)準(zhǔn)是：當(dāng)2個(gè)殘基的Cβ原子之間歐氏距離小于8?時(shí)，則認(rèn)為這2個(gè)殘基具有相互作用，稱之為殘基間接觸(contact)，而所有的接觸則形象地表示成接觸圖譜(contact map，見圖2).形式上，一個(gè)包含L個(gè)殘基的蛋白質(zhì)接觸圖譜可以表示成一個(gè)L×L的矩陣A*：

(1)

Fig. 2 Illustration of residue-residue contact in a protein 1a3a_A圖2 蛋白質(zhì)1a3a_A中殘基接觸示意圖

在獲得了殘基間距離信息之后，采用分子動(dòng)力學(xué)模擬(molecular dynamics simulation, MDS)[9]等技術(shù)可以有效地反推出各個(gè)殘基的空間位置.因此，殘基間相互作用的準(zhǔn)確預(yù)測(cè)成為蛋白質(zhì)空間結(jié)構(gòu)預(yù)測(cè)的關(guān)鍵環(huán)節(jié)之一.

在蛋白質(zhì)進(jìn)化過程，有相互作用殘基對(duì)之間存在一種“共進(jìn)化”模式，即當(dāng)一個(gè)殘基發(fā)生變異時(shí)，與其有相互作用的殘基也要發(fā)生相應(yīng)的變異，以維持相互作用，進(jìn)而維持整體空間結(jié)構(gòu)以及生物學(xué)功能.基于上述生物學(xué)觀察，研究者提出了多種統(tǒng)計(jì)模型和算法以預(yù)測(cè)殘基間相互作用.從統(tǒng)計(jì)學(xué)的角度講，由蛋白質(zhì)序列信息出發(fā)預(yù)測(cè)殘基間遠(yuǎn)程相互作用是一個(gè)典型的機(jī)器學(xué)習(xí)問題，即預(yù)測(cè)組成單元間關(guān)聯(lián)關(guān)系的結(jié)構(gòu)學(xué)習(xí)(structured learning)問題[10].

目前的預(yù)測(cè)方法主要分為2類：無監(jiān)督學(xué)習(xí)方法和監(jiān)督學(xué)習(xí)方法.簡(jiǎn)要地說，無監(jiān)督學(xué)習(xí)方法僅從序列出發(fā)抽取出待測(cè)蛋白質(zhì)的進(jìn)化歷史信息，進(jìn)而分析2個(gè)殘基在進(jìn)化過程中的共變程度，以共變程度的強(qiáng)弱來推斷殘基間是否存在相互作用.由于這類方法不依賴于已知相互作用的蛋白質(zhì)集合，因此屬于無監(jiān)督學(xué)習(xí)方法的范疇.另一類方法是基于已有的結(jié)構(gòu)信息，依據(jù)每個(gè)殘基的序列特征和結(jié)構(gòu)特征，采用神經(jīng)網(wǎng)絡(luò)[11]、支持向量機(jī)[12-13]等分類模型預(yù)測(cè)殘基間是否存在相互作用，從而屬于有監(jiān)督學(xué)習(xí)方法的范疇.

值得強(qiáng)調(diào)指出的是在生物信息學(xué)這門交叉學(xué)科中，重要的是如何將計(jì)算模型和生物學(xué)現(xiàn)象相結(jié)合.具體到殘基相互作用這個(gè)問題而言:1)殘基間相互作用預(yù)測(cè)大量應(yīng)用了現(xiàn)有的統(tǒng)計(jì)、組合最優(yōu)化、機(jī)器學(xué)習(xí)等領(lǐng)域的研究成果，是對(duì)現(xiàn)有成果的應(yīng)用和檢驗(yàn)；2)殘基間相互作用預(yù)測(cè)問題有其特殊性，主要體現(xiàn)于蛋白質(zhì)是進(jìn)化作用的結(jié)果.因此，在建模過程中不能簡(jiǎn)單照搬現(xiàn)有算法和模型，而是需要考慮進(jìn)化等生物學(xué)觀察，對(duì)現(xiàn)有算法和模型作必要的擴(kuò)展和改進(jìn).

換句話說，每一種統(tǒng)計(jì)模型和算法的設(shè)計(jì)都是基于特定的生物學(xué)觀察基礎(chǔ)之上的，是對(duì)生物學(xué)觀察的數(shù)學(xué)刻畫和描述；另一方面，我們對(duì)殘基相互作用的生物學(xué)觀察越深刻，則越有助于我們?cè)O(shè)計(jì)更有效的統(tǒng)計(jì)模型.

依據(jù)上述觀點(diǎn)，本文在介紹每一類算法時(shí)，都首先介紹生物學(xué)觀察，進(jìn)而介紹如何基于上述生物學(xué)觀察設(shè)計(jì)統(tǒng)計(jì)模型.

本文綜述了目前的蛋白質(zhì)殘基遠(yuǎn)程相互作用預(yù)測(cè)算法:1)介紹無監(jiān)督學(xué)習(xí)方法(又進(jìn)一步細(xì)分為局部模型和全局模型2類);2)介紹監(jiān)督學(xué)習(xí)方法，并以國(guó)際蛋白質(zhì)結(jié)構(gòu)預(yù)測(cè)競(jìng)賽(critical assessment of protein structure prediction, CASP)的結(jié)果來分析比較現(xiàn)有方法；3)總結(jié)分析未來的發(fā)展趨勢(shì).

1 無監(jiān)督學(xué)習(xí)方法的局部模型

1.1 基本思想

1) 生物學(xué)觀察.從進(jìn)化的角度來看，同源蛋白質(zhì)是指由同一個(gè)祖先蛋白質(zhì)進(jìn)化而來的后代蛋白質(zhì)，其結(jié)構(gòu)和功能有一定的保守性，這種保守性由殘基的遠(yuǎn)程相互作用來維持.在進(jìn)化過程中，當(dāng)有相互作用的殘基對(duì)中其中一個(gè)殘基發(fā)生突變時(shí)，另一個(gè)常常也發(fā)生相應(yīng)的突變，以維持相互作用，否則不利于蛋白質(zhì)整體結(jié)構(gòu)的穩(wěn)定.這種現(xiàn)象被稱之為“有相互作用殘基對(duì)的共進(jìn)化”,如圖3所示.

2) 數(shù)學(xué)模型設(shè)計(jì).通常采用多序列聯(lián)配的數(shù)學(xué)形式來刻畫同源蛋白質(zhì)之間的同源關(guān)系，其中每一列表示由祖先蛋白質(zhì)中的某個(gè)殘基進(jìn)化生成的殘基(如圖4所示).考慮長(zhǎng)度為L(zhǎng)的多序列聯(lián)配，通常用離散隨機(jī)變量Xi(1≤i≤L)表示第i個(gè)殘基(亦稱為位點(diǎn))的氨基酸種類，可取21個(gè)離散值(包括20種氨基酸和一個(gè)聯(lián)配空位)，則多序列聯(lián)配中的每條序列可看成這些變量的一個(gè)觀測(cè)樣本.

Fig. 3 The phylogenetic tree and MSA of PF00111圖3 序列家族PF00111的進(jìn)化樹和多序列聯(lián)配

Fig. 4 Illustration of the principal for contact prediction using unsupervised methods圖4 無監(jiān)督學(xué)習(xí)方法預(yù)測(cè)殘基遠(yuǎn)程相互作用原理示意圖[14]

基于殘基共進(jìn)化現(xiàn)象，無監(jiān)督學(xué)習(xí)方法的基本思想是首先檢索出與待預(yù)測(cè)蛋白質(zhì)序列同源的所有序列，計(jì)算出多序列聯(lián)配(multiple sequence alignment, MSA)，以此來表示待測(cè)蛋白質(zhì)的進(jìn)化歷史信息；進(jìn)而分析2個(gè)殘基在進(jìn)化過程中的共進(jìn)化(或者共變異)程度，以共變程度的強(qiáng)弱來推斷殘基間是否存在相互作用.從統(tǒng)計(jì)的角度來講，即是通過多序列聯(lián)配中列向量之間的相關(guān)性，即2個(gè)隨機(jī)變量Xi,Xj之間的相關(guān)性，從而推斷殘基之間的相互作用.

無監(jiān)督學(xué)習(xí)方法可以分為兩大類，即局部模型和全局模型[6].其中,局部模型假設(shè)一個(gè)殘基對(duì)內(nèi)部的相關(guān)性與其他殘基對(duì)是獨(dú)立的,從而每對(duì)殘基單獨(dú)計(jì)算其相關(guān)性;而全局模型則考慮了殘基對(duì)之間的關(guān)聯(lián)關(guān)系，對(duì)所有的殘基建立統(tǒng)一的全局模型.我們?cè)诒竟?jié)介紹局部模型，在第2節(jié)介紹全局模型.

1.2 典型方法

局部模型在計(jì)算某對(duì)殘基之間的相關(guān)性時(shí)，不考慮其他殘基對(duì)的影響，直接分別計(jì)算各殘基對(duì)之間的相關(guān)性；各種局部模型的差異主要體現(xiàn)在變量相關(guān)性的衡量方法不同.下面我們介紹3種典型的局部模型，并分析其優(yōu)缺點(diǎn).

1.2.1 典型方法1：共變相關(guān)系數(shù)

(2)

其中，M是多序列聯(lián)配的序列條數(shù)，〈·〉為矩陣元素均值，σi為矩陣si所有元素的標(biāo)準(zhǔn)差，如圖5所示.

Fig. 5 Illustration of co-mutation extraction for residue pair (i,j)[15]圖5 位點(diǎn)對(duì)(i,j)共變信息計(jì)算示意圖[15]

實(shí)驗(yàn)表明:將相關(guān)系數(shù)作為共變度量，在一定閾值下推斷殘基間相互作用，準(zhǔn)確率比隨機(jī)預(yù)測(cè)有顯著提高，從而表明由殘基共變性推斷其相互作用的可行性；然而Pearson相關(guān)系數(shù)只能度量隨機(jī)變量間的線性相關(guān)關(guān)系，從而使得該方法存在一定的局限性.

1.2.2 典型方法2：互信息

Martin等人[16]用互信息(mutual information, MI)識(shí)別共進(jìn)化殘基對(duì).這種方法源于信息論，具體來說，對(duì)于某個(gè)多序列聯(lián)配中的位點(diǎn)i和j，其互信息定義為

(3)

其中,fij(a,b)為位點(diǎn)i出現(xiàn)殘基a且位點(diǎn)j出現(xiàn)殘基b的頻率，fi(a)表示位點(diǎn)i出現(xiàn)殘基a的頻率.與Pearson相關(guān)系數(shù)相比而言，互信息可以度量變量之間的非線性關(guān)系，其值越大表示殘基對(duì)間的共進(jìn)化程度越大，互信息為0則表示2位點(diǎn)獨(dú)立進(jìn)化，或存在保守位點(diǎn).

Martin等人的實(shí)驗(yàn)結(jié)果表明互信息較高的位點(diǎn)對(duì)傾向于具有相互作用，但是其效果受限于序列條數(shù)和進(jìn)化偏差造成的背景噪音.因而欲提高預(yù)測(cè)準(zhǔn)確率，需降低背景噪聲的影響，減少對(duì)序列條數(shù)的依賴，從而為后續(xù)研究指明了改進(jìn)方向.

1.2.3 典型方法3：OMES

Kass等人[17]提出另外一種共變性的度量方法OMES(observed minus expected squared)，這種方法基于統(tǒng)計(jì)學(xué)中的卡方檢驗(yàn)，通過比較殘基對(duì)在2個(gè)位點(diǎn)上實(shí)際出現(xiàn)次數(shù)與期望出現(xiàn)次數(shù)之間的差異來定量刻畫殘基對(duì)的共進(jìn)化程度.具體地，其定義為

(4)

其中，Oij(a,b)和Eij(a,b)分別表示殘基對(duì)(a,b)出現(xiàn)頻數(shù)的觀測(cè)值和期望值，M是序列條數(shù).Oij(a,b)可以直接從多序列聯(lián)配中統(tǒng)計(jì)得到；Eij(a,b)是在假設(shè)殘基對(duì)間不存在相關(guān)性的前提下計(jì)算得到的，即Eij(a,b)=Mfi(a)fj(b)，其中fi(a),fj(b)分別表示相應(yīng)位點(diǎn)某氨基酸的出現(xiàn)頻率.OMES的值越大表示2個(gè)位點(diǎn)之間的共進(jìn)化程度越高；對(duì)于2個(gè)完全獨(dú)立進(jìn)化的位點(diǎn)，OMES的值為0.

我們?cè)贕REMLIN數(shù)據(jù)集[18]上測(cè)試OMES方法，并與MI進(jìn)行了比較；實(shí)驗(yàn)結(jié)果表明：MI與OMES的預(yù)測(cè)性能相當(dāng),詳細(xì)實(shí)驗(yàn)結(jié)果分析參見第4節(jié).

1.3 局部模型的實(shí)驗(yàn)結(jié)果及分析

局部模型的優(yōu)點(diǎn)是簡(jiǎn)單、計(jì)算速度快；但是也有較大的不足，主要表現(xiàn)為：1)各殘基對(duì)之間并不是獨(dú)立的，而是存在關(guān)聯(lián)傳遞的現(xiàn)象，局部模型并沒有考慮這種關(guān)聯(lián)傳遞現(xiàn)象；2)未考慮序列空間采樣偏差及樣本數(shù)不足的影響；3)相關(guān)性計(jì)算存在大量由進(jìn)化偏差產(chǎn)生的背景噪聲.

雖然局部模型普遍存在預(yù)測(cè)準(zhǔn)確率偏低的缺陷，然而在一定程度上提取了進(jìn)化信息，是由序列信息推斷結(jié)構(gòu)約束的早期嘗試，對(duì)后續(xù)研究具有重要的啟發(fā)和借鑒意義.

2 無監(jiān)督學(xué)習(xí)方法的全局模型

2.1 基本思想

1) 生物學(xué)觀察.局部模型單獨(dú)計(jì)算各個(gè)殘基對(duì)的相關(guān)性，其暗含的假設(shè)是某個(gè)殘基對(duì)的相互作用是獨(dú)立于其他殘基對(duì)的.然而一個(gè)殘基可能與多個(gè)殘基有作用，從而導(dǎo)致關(guān)聯(lián)傳遞這一普遍存在的現(xiàn)象.如圖6所示，如果殘基A和殘基D共變，殘基D和殘基C共變，那么從序列信息來看，殘基A和殘基D也表現(xiàn)出共變性，然而殘基A和殘基D之間的共變性源于傳遞效應(yīng)，并非源自殘基A和殘基D的相互作用.這種通過殘基共進(jìn)化的傳遞效應(yīng)導(dǎo)致的相關(guān)性稱為間接關(guān)聯(lián).

Fig. 6 Illustration of direct and indirect couplings圖6 直接關(guān)聯(lián)和間接關(guān)聯(lián)示意圖[25]

2) 數(shù)學(xué)模型設(shè)計(jì).局部模型假設(shè)任意2個(gè)殘基的共進(jìn)化和其他殘基是相互獨(dú)立的，決定了它只能探測(cè)相關(guān)性，會(huì)受到間接關(guān)聯(lián)噪聲的顯著影響，從而不能準(zhǔn)確提取出真實(shí)共進(jìn)化(接觸)的殘基對(duì).鑒于此，全局模型對(duì)所有位點(diǎn)建立全概率模型，同時(shí)考慮所有殘基對(duì)之間的關(guān)聯(lián)關(guān)系，試圖去除間接關(guān)聯(lián)的影響，從而避免局部模型的缺陷.

2.2 典型方法

迄今為止已經(jīng)發(fā)展了多種全局模型，比如Markov隨機(jī)場(chǎng)模型(Markov random fields, MRF)[19-21]、貝葉斯網(wǎng)絡(luò)模型(Bayesian network)[22]、高斯圖模型(Gaussian graphical model)[23-24]和網(wǎng)絡(luò)反卷積(network deconvolution)等[25].這些方法的不同之處主要體現(xiàn)在如何對(duì)多序列聯(lián)配建模，其中貝葉斯網(wǎng)絡(luò)模型采用有向圖模型進(jìn)行建模；Markov隨機(jī)場(chǎng)模型和高斯圖模型采用無向圖模型進(jìn)行建模；高斯圖模型可以看成Markov隨機(jī)場(chǎng)模型的特殊形式.

2.2.1 典型方法1：貝葉斯網(wǎng)絡(luò)模型

Burger等人提出使用貝葉斯網(wǎng)絡(luò)模型預(yù)測(cè)殘基間相互作用[22].在這種方法中，使用貝葉斯網(wǎng)絡(luò)把殘基間的共進(jìn)化關(guān)聯(lián)關(guān)系表示成依賴關(guān)系：當(dāng)位點(diǎn)i和位點(diǎn)j存在相互作用時(shí)，則第i個(gè)位置出現(xiàn)殘基Xi的概率依賴于第j個(gè)位置出現(xiàn)殘基Xj的概率.這種依賴關(guān)系形象地表示成貝葉斯網(wǎng)絡(luò)中的一條有向邊,如圖7所示.

Fig. 7 Bayesian network model of a given MSA圖7 貝葉斯網(wǎng)絡(luò)模型示意圖

在已知位點(diǎn)之間依賴關(guān)系的情況下，可以計(jì)算觀察到某個(gè)多序列聯(lián)配的條件概率；反過來，在給定多序列聯(lián)配的情況下，結(jié)合依賴關(guān)系的先驗(yàn)分布，可以推斷位點(diǎn)i和位點(diǎn)j之間存在依賴關(guān)系(共進(jìn)化)的后驗(yàn)概率，最后認(rèn)為后驗(yàn)概率高的殘基對(duì)具有相互作用.

本節(jié)首先描述給定依賴關(guān)系的情況下觀察到某個(gè)多序列聯(lián)配的條件概率計(jì)算過程，然后介紹依賴關(guān)系的后驗(yàn)概率的計(jì)算過程.

1) 已知?dú)埢g依賴關(guān)系計(jì)算多序列聯(lián)配的條件概率

假設(shè)殘基間所有的依賴關(guān)系形成有向圖T=(π,V)，其中π表示所有有向邊的集合，V={1,2,…,L}表示殘基位點(diǎn)集合.為簡(jiǎn)化計(jì)算，進(jìn)一步假設(shè)依賴關(guān)系圖T是樹狀圖，即在T中存在唯一根節(jié)點(diǎn)r,除根節(jié)點(diǎn)r外，其余節(jié)點(diǎn)i都存在唯一父節(jié)點(diǎn)，記其父節(jié)點(diǎn)為π(i).

給定殘基間依賴關(guān)系T，則觀察到多序列聯(lián)配D的條件概率為

(5)

其中,Sij=P(Di,Dj)(P(Di)P(Dj))；，表示MSA第i列中各氨基酸出現(xiàn)頻率向量.假設(shè)第i列中氨基酸α出現(xiàn)的概率為，且wi服從Dirichlet分布，則:

P(Di)=∫P(Di|w)P(w)dw=

(6)

同理可得:

(7)

其中參數(shù)λ,λ′是偽計(jì)數(shù).

2) 殘基間依賴關(guān)系后驗(yàn)概率的計(jì)算

(8)

則式(8)可重寫為

(9)

由式(9)可推出MSA的概率模型為

，

(10)

其中,Bjπ(j)=Sjπ(j)Wjπ(j).

使用貝葉斯公式計(jì)算T的后驗(yàn)分布為:P(T|D)=P(D|T)P(T)P(D),進(jìn)而對(duì)于特定的殘基對(duì)(k,l)之間有邊Ekl的后驗(yàn)概率通過對(duì)包含Ekl的所有樹的后驗(yàn)概率求和得到,如圖8所示:

Fig. 8 Illustration of the calculation of posterior probability圖8 后驗(yàn)概率計(jì)算示意圖[22]

(11)

實(shí)驗(yàn)結(jié)果表明:以后驗(yàn)概率P(Eij|D)作為殘基對(duì)(i,j)的相關(guān)度量能夠去掉間接關(guān)聯(lián)，和局部模型相比，顯著提高了預(yù)測(cè)準(zhǔn)確率.

2.2.2 典型方法2：Markov隨機(jī)場(chǎng)

Markov隨機(jī)場(chǎng)是一種無向圖模型，其形式可由最大熵原理推導(dǎo)得到，所以也被稱為最大熵模型.Markov隨機(jī)場(chǎng)的優(yōu)勢(shì)在于可以直接刻畫殘基間的遠(yuǎn)程相互作用[19-21,27].

設(shè)多序列聯(lián)配的長(zhǎng)度為L(zhǎng)，變量Xi表示第i個(gè)位置出現(xiàn)的氨基酸，則多序列聯(lián)配中的一條序列的生成概率為

P(X1,X2,…,XL)=

(12)

其中:

為配分函數(shù)，eij(Xi,Xj)表示位置i處氨基酸Xi和位置j處氨基酸Xj的耦合強(qiáng)度，hi(Xi)為位置i處觀察到殘基Xi的可能性，均為待確定的參數(shù).最終的推斷規(guī)則為：耦合強(qiáng)度強(qiáng)的殘基對(duì)被預(yù)測(cè)為具有相互作用.

假設(shè)給定包含M條序列的多序列聯(lián)配，上述待定參數(shù)eij(Xi,Xj)與hi(Xi)可以通過極大似然策略進(jìn)行估計(jì).然而上述極大似然策略涉及到配分函數(shù)Z的計(jì)算，其計(jì)算是NP難問題，所以發(fā)展出多種近似求解方法，包括置信傳播算法、平均場(chǎng)近似算法和偽似然最大化算法，簡(jiǎn)要介紹如下：

1) 置信傳播算法bpDCA

Weigt等人[20]用Markov隨機(jī)場(chǎng)模型研究蛋白質(zhì)-蛋白質(zhì)相互作用，并用置信傳播算法(bpDCA)近似求解模型參數(shù)，后來這種方法也被直接借用于殘基間相互作用的推測(cè).

置信傳播算法的基本思想是通過局部信息的多次傳播以逼近全局信息，以此計(jì)算配分函數(shù)或邊際概率.確切地說，在最大化似然函數(shù)過程中，梯度函數(shù)的計(jì)算涉及邊際概率，而置信傳播算法的核心是解決梯度計(jì)算問題.在給定初始參數(shù)的情況下，bpDCA迭代執(zhí)行2個(gè)步驟直至滿足收斂條件：

① 用置信傳播算法估計(jì)邊緣分布Pi(Xi)和Pij(Xi,Xj)

首先對(duì)于每個(gè)位置i，迭代求解信息傳遞Pi→j(Xi)：

(13)

此處fi(A)為經(jīng)驗(yàn)頻率.

然后可獲得邊際分布Pi(Xi)的估計(jì)：

Pi(Xi)～exp{hi(Xi)}

(14)

用類似的方法可得Pij(Xi,Xj)的估計(jì).

② 用梯度下降策略更新參數(shù)估計(jì)

首先，似然函數(shù)的梯度可估計(jì)為

Δeij(Xi,Xj)=fij(Xi,Xj)-Pij(Xi,Xj)-

Δhi(Xi)=fi(Xi)-Pi(Xi)，

(15)

其中,fi(A)和fij(A,B)為經(jīng)驗(yàn)頻率.然后更新參數(shù)估計(jì)為

這里ρ為迭代步長(zhǎng).

bpDCA主要有2個(gè)缺陷：①速度慢.該算法迭代1次的計(jì)算復(fù)雜度為O(212L4)，即使對(duì)長(zhǎng)度為60的短蛋白，bpDCA在4核CPU上也需大約運(yùn)行4 d.②收斂性差.該算法解的漸進(jìn)性質(zhì)不能得到有效控制，理論上不能保證其收斂性.

2) 平均場(chǎng)近似算法mfDCA

Morcos等人[19]提出使用平均場(chǎng)近似策略(mfDCA)來近似求解Markov隨機(jī)場(chǎng)的參數(shù).平均場(chǎng)近似的基本思想是由簡(jiǎn)單可分解分布近似復(fù)雜分布，因此其核心在于2個(gè)問題：①如何確定簡(jiǎn)單分布的形式;②如何衡量簡(jiǎn)單分布和復(fù)雜分布之間的差異，并找到最接近原始復(fù)雜分布的簡(jiǎn)單分布.

(16)

其中C可由經(jīng)驗(yàn)協(xié)方差矩陣計(jì)算，即Cij(A,B)=fij(A,B)-fi(A)fj(B).此處fi(A)和fij(A,B)為經(jīng)驗(yàn)頻率.

mfDCA的優(yōu)勢(shì)是速度快，通過求逆計(jì)算耦合參數(shù)的時(shí)間復(fù)雜度是O(L3)，比bpDCA速度提高上千倍，從而使大量蛋白質(zhì)家族的計(jì)算成為可能.

3) 極大偽似然算法plmDCA

L?vkvist等人[21]提出偽似然最大化方法估計(jì)MRF的參數(shù)，其基本思想是用偽似然函數(shù)近似似然函數(shù).由于計(jì)算偽似然函數(shù)梯度的時(shí)間復(fù)雜度是多項(xiàng)式的，所以可以有效地估計(jì)參數(shù).

偽似然函數(shù)定義如下：

(17)

上述模型參數(shù)數(shù)目過多，對(duì)長(zhǎng)為L(zhǎng)的序列來說，模型參數(shù)規(guī)模為(20×20)×L(L-1)2+20L.當(dāng)L=100時(shí)，模型將有近200萬的參數(shù).為避免過擬合問題，Ekeberg等人在偽似然函數(shù)中引入了正則項(xiàng)R(h,e)，即通過解決以下優(yōu)化問題求解參數(shù)：

{hPLM,ePLM}=arg min{lpseudo(h,e)+R(h,e)}，

(18)

其中，λh和λe分別表示單體項(xiàng)h和雙體項(xiàng)e的正則化參數(shù).

該方法避免了極大似然求解復(fù)雜配分函數(shù)的問題且當(dāng)樣本量足夠大時(shí)，極大偽似然估計(jì)是極大似然估計(jì)的一致估計(jì)[30]，從而能夠保證獲得準(zhǔn)確的參數(shù)估計(jì)，并且和mfDCA相比顯著提高了準(zhǔn)確率.

Kamisetty等人[18]在極大偽似然的基礎(chǔ)上，進(jìn)一步將結(jié)構(gòu)先驗(yàn)信息引入正則項(xiàng)，開發(fā)軟件GREMLIN.實(shí)驗(yàn)結(jié)果表明，由于引入結(jié)構(gòu)先驗(yàn)信息，GREMLIN方法的性能優(yōu)于plmDCA.

2.2.3 典型方法3：高斯圖模型

高斯圖模型假設(shè)多序列聯(lián)配中每一條蛋白質(zhì)序列都服從高斯分布N(μ,Σ)，其中高斯分布的協(xié)方差矩陣的逆稱為精細(xì)矩陣(recision matrix)，記作Θ=Σ-1.精細(xì)矩陣表征了變量之間的直接關(guān)聯(lián)信息[31]，因而可以通過精細(xì)矩陣來預(yù)測(cè)殘基間的相互作用.在統(tǒng)計(jì)學(xué)中，通過分析精細(xì)矩陣來推斷直接關(guān)聯(lián)的策略也稱為偏相關(guān)分析.

為防止過擬合，通常采用引入正則項(xiàng)的策略來控制模型的復(fù)雜度.典型的方法包括PSICOV所使用的圖Lasso以及CoinDCA所采用的成組Lasso，簡(jiǎn)要介紹如下：

1) PSICOV使用的圖Lasso策略

Jones等人[24]利用圖Lasso策略推斷精細(xì)矩陣，開發(fā)了軟件PSICOV.該方法的核心思想是優(yōu)化含有正則項(xiàng)的似然函數(shù)：

(19)

其中,COV為經(jīng)驗(yàn)協(xié)方差陣.前2項(xiàng)為高斯分布的對(duì)數(shù)似然，第3項(xiàng)為正則項(xiàng).正則項(xiàng)的引入有2個(gè)作用：①控制模型復(fù)雜度，防止過擬合，以避免模型參數(shù)過多導(dǎo)致參數(shù)推斷的困難；②保證精細(xì)矩陣的稀疏性，以此刻畫接觸圖譜的稀疏性.

Jones等人[24]在150個(gè)目標(biāo)蛋白進(jìn)行測(cè)試，結(jié)果顯示PSICOV一致優(yōu)于局部模型和貝葉斯網(wǎng)絡(luò)模型.

2) CoinDCA采用的成組Lasso策略

Ma等人[23]對(duì)高斯圖模型做了擴(kuò)展，以融合多個(gè)相關(guān)家族的進(jìn)化信息，并開發(fā)了軟件CoinDCA.

CoinDCA的基本思想是:假設(shè)對(duì)于目標(biāo)蛋白質(zhì)序列，與其具有相同折疊類型的共有K個(gè)家族(可通過同源搜索獲得)，由于這K個(gè)蛋白家族屬于同一折疊類型，所以可認(rèn)為它們有類似的蛋白質(zhì)接觸圖譜；相應(yīng)地，當(dāng)用K個(gè)高斯分布分別對(duì)K個(gè)蛋白家族建模時(shí)，它們有類似的精細(xì)矩陣；成組Lasso的目的是約束這K個(gè)高斯分布具有類似的精細(xì)矩陣.

具體地，通過優(yōu)化下式求解這些精細(xì)矩陣.

(20)

CoinDCA充分利用了相近家族的進(jìn)化信息且融合了監(jiān)督學(xué)習(xí)(隨機(jī)森林)的方法.在PSICOV，CASP10和CASP11數(shù)據(jù)集上測(cè)試的實(shí)驗(yàn)結(jié)果表明，這種方法對(duì)同源序列少的蛋白預(yù)測(cè)準(zhǔn)確率有顯著提高，減少了對(duì)同源序列數(shù)目的依賴；而單一地整合多家族信息或者采用隨機(jī)森林的方法并不能對(duì)預(yù)測(cè)性能有所改進(jìn).

2.2.4 典型方法4：網(wǎng)絡(luò)反卷積

在殘基相互作用中消除間接效應(yīng)，本質(zhì)上是網(wǎng)絡(luò)推斷領(lǐng)域直接作用和間接作用的區(qū)分問題[32].Feizi等人[25]提出網(wǎng)絡(luò)反卷積(network deconvolution, ND)策略推斷直接相互作用網(wǎng)絡(luò).

網(wǎng)絡(luò)反卷積的基本思想是：假設(shè)觀測(cè)到的相關(guān)矩陣Gobs是直接相關(guān)矩陣Gdir和間接相關(guān)矩陣的疊加，而間接相關(guān)可視為直接相關(guān)通過多步傳遞得到的(如圖9所示)，即：

注意到當(dāng)矩陣Gdir特征值絕對(duì)值小于1時(shí)，等式右邊收斂，上式有閉合形式：

(21)

Fig. 9 Illustration of network deconvolution圖9 網(wǎng)絡(luò)反卷積意圖[25]

網(wǎng)絡(luò)反卷積方法廣泛應(yīng)用于社交網(wǎng)絡(luò)、基因調(diào)控網(wǎng)絡(luò)等領(lǐng)域的推斷中.Wright等人[32]在殘基相互作用推斷中的結(jié)果表明，網(wǎng)絡(luò)反卷積策略可有效過濾掉互信息的間接關(guān)聯(lián)噪聲，而對(duì)全局模型mfDCA輸出的直接信息(已去除間接關(guān)聯(lián))進(jìn)行反卷積并沒有顯著改進(jìn)；而對(duì)互信息矩陣反卷積的預(yù)測(cè)效果不如mfDCA.上述結(jié)果說明網(wǎng)絡(luò)反卷積的策略具有普適性，但對(duì)于特定的殘基相互作用預(yù)測(cè)問題則仍然需要進(jìn)行相應(yīng)的改進(jìn).

Sun等人[33]提出了平衡網(wǎng)絡(luò)反卷積方法，該方法不需要像原始網(wǎng)絡(luò)反卷積方法那樣對(duì)Gobs進(jìn)行線性縮放,其假設(shè)

Gobs=Gdir+Gindir=

(22)

則:

(23)

可以導(dǎo)出對(duì)任意的λobs∈(-∞,+∞)，都有|λdir|<1，所以該方法不需要對(duì)Gobs進(jìn)行線性縮放.結(jié)果表明，平衡網(wǎng)絡(luò)反卷積方法的預(yù)測(cè)性能優(yōu)于原始網(wǎng)絡(luò)反卷積方法；如何進(jìn)一步提升預(yù)測(cè)性能，仍然需要后續(xù)有針對(duì)性的研究.

2.3 全局模型的實(shí)驗(yàn)結(jié)果及分析

我們?cè)贕REMLIN數(shù)據(jù)集上詳細(xì)測(cè)試了無監(jiān)督學(xué)習(xí)方法的性能，如表1、表2所示，結(jié)果顯示全局模型優(yōu)于局部模型，且在全局模型中plmDCA預(yù)測(cè)性能最優(yōu).具體地，我們從實(shí)驗(yàn)結(jié)果分析中獲得3個(gè)結(jié)論：

1) 由于有效地去除了間接效應(yīng)的影響，全局模型比局部模型革命性地提高了預(yù)測(cè)準(zhǔn)確率.

2) 全局模型之間的預(yù)測(cè)性能差距相對(duì)較小，且不同方法的預(yù)測(cè)結(jié)果有一定程度的不同，將不同方法得到的殘基相互作用信息有效地整合，并用于蛋白質(zhì)三級(jí)結(jié)構(gòu)預(yù)測(cè)，具有重要意義和廣泛的應(yīng)用前景.

Table 1 Denoising Performance of LRS for Three LocalMethods on GREMLIN Benchmark

Table 2 Denoising Performance of LRS for Three GlobalMethods on GREMLIN Benchmark

3) 全局模型普遍存在參數(shù)多的問題，要得到精確的參數(shù)估計(jì)需要大量樣本信息，預(yù)測(cè)效果強(qiáng)烈依賴于同源序列的數(shù)目；且復(fù)雜的參數(shù)估計(jì)導(dǎo)致全局模型運(yùn)行速度較慢，需要開發(fā)更加有效的參數(shù)估計(jì)方法.

3 無監(jiān)督學(xué)習(xí)方法中的預(yù)處理和后處理

除了統(tǒng)計(jì)模型之外，影響遠(yuǎn)程相互作用預(yù)測(cè)性能的因素還包括樣本的不獨(dú)立性以及背景噪聲的影響等，因此需要采取預(yù)處理和后處理步驟以消除這些因素的影響.通常采用的預(yù)處理和后處理步驟簡(jiǎn)要陳述如下.

3.1 克服樣本不獨(dú)立性的預(yù)處理過程

多序列聯(lián)配中的樣本是與待測(cè)蛋白同源的序列，這些序列來源于同一祖先蛋白質(zhì)，從而造成了觀測(cè)樣本之間的不獨(dú)立性，影響模型預(yù)測(cè)的準(zhǔn)確性.為提高準(zhǔn)確性，一般對(duì)多序列聯(lián)配進(jìn)行2方面預(yù)處理：

1) 去冗余.比如去掉與目標(biāo)蛋白質(zhì)高度相似(通常采用序列等同度大于90%)的序列.

2) 加權(quán)重.對(duì)于每條序列，都依據(jù)在多序列聯(lián)配中與其相似的序列條數(shù)賦予權(quán)重，其基本思想是：如果一條蛋白質(zhì)序列具有較多的相似序列，則權(quán)重較低；反之則設(shè)置較高權(quán)重.具體地，對(duì)第k條序列來說，首先統(tǒng)計(jì)與其序列等同度高于80%的序列數(shù)目：

mk=|{n∈{1,2,…,M}|seqid(Xk,Xn)>80%}|，

(24)

進(jìn)而在統(tǒng)計(jì)殘基和殘基對(duì)頻率向量時(shí)，將每條序列的權(quán)重設(shè)置為1mk.

3.2 去背景噪聲的后處理過程

同源序列中通常包含由進(jìn)化造成的較強(qiáng)的背景噪聲.具體來說，如果一個(gè)位點(diǎn)突變發(fā)生在進(jìn)化早期，其后代都將延續(xù)這個(gè)突變，從而導(dǎo)致過高地估計(jì)了此位點(diǎn)和其他位點(diǎn)之間的共變性.

通常采用后處理的方法消除這些背景噪聲對(duì)相互作用預(yù)測(cè)的影響，常用的策略簡(jiǎn)要介紹如下：

1) 均值乘積校正方法(average product correction, APC)

Dunn等人[34]基于信息論提出對(duì)互信息矩陣進(jìn)行歸一化去噪的方法APC.其基本思想是假設(shè)背景噪聲具有如下的均值乘積的形式：

(25)

其中MIi.表示位點(diǎn)i與其他位點(diǎn)互信息的平均值；MI..表示所有位點(diǎn)對(duì)互信息的平均值.

經(jīng)APC去噪后的互信息MIp為

MIpij=MIij-APCij.

(26)

實(shí)驗(yàn)結(jié)果表明采用APC技術(shù)去除背景噪聲后，能夠有效提高基于互信息對(duì)殘基相互作用預(yù)測(cè)的精度.目前，該策略被推廣應(yīng)用到其他相關(guān)矩陣，已成為全局統(tǒng)計(jì)模型標(biāo)準(zhǔn)的后處理步驟.

2) 譜去除方法(spectrum cleaning, SC)

3) 低秩稀疏矩陣分解方法(low rank and sparse matrix decomposition, LRS)

譜去除方法假設(shè)背景噪聲來源于相關(guān)性矩陣的第一主成分，其秩為1；然而當(dāng)多序列聯(lián)配中的序列是來源于多個(gè)家族，則相關(guān)性矩陣的其他主成分也會(huì)含有背景噪聲[35].基于上述認(rèn)識(shí)和觀察，我們團(tuán)隊(duì)假設(shè)背景噪聲是低秩的，同時(shí)真實(shí)的相互作用是稀疏的.大量統(tǒng)計(jì)數(shù)據(jù)支持上述假設(shè)的合理性，即真實(shí)相互作用僅占所有可能相互作用的3%～5%[11,24].基于上述生物學(xué)觀察，我們?cè)O(shè)計(jì)了統(tǒng)計(jì)模型將背景噪聲和真實(shí)的相互作用信號(hào)區(qū)分開來[36].

具體地，對(duì)于給定的殘基相關(guān)性度量矩陣M，我們認(rèn)為它是低秩噪聲矩陣和稀疏信號(hào)矩陣的疊加，進(jìn)而用低秩稀疏矩陣分解技術(shù)還原真實(shí)的相互作用信號(hào)矩陣，即求解如下的優(yōu)化問題：

(27)

4) 去背景噪音方法的實(shí)驗(yàn)結(jié)果及分析

我們?cè)贕REMLIN測(cè)試集上測(cè)試去背景噪音方法的性能，詳細(xì)實(shí)驗(yàn)結(jié)果表1、表2所示.從表1和表2中可以看出無論序列間隔取值，無論是局部模型或全局模型的具體方法，LRS的去噪效果都一致地優(yōu)于APC技術(shù)，且對(duì)局部模型的改善顯著高于全局模型.LRS技術(shù)的價(jià)值集中體現(xiàn)在對(duì)局部模型的顯著改進(jìn)，使其達(dá)到和全局模型mfDCA相近的性能.這是自全局模型提出以來首次發(fā)現(xiàn)基于局部度量的方法能夠達(dá)到和全局模型可比的效果，也說明只有當(dāng)有效地去除背景噪音之后，相關(guān)性度量才能提取出更加準(zhǔn)確的共進(jìn)化信息.下面我們將從理論上進(jìn)一步深入分析LRS優(yōu)于其他去背景噪音方法的原因，主要基于2個(gè)事實(shí)：

① APC和SC的等價(jià)性

SC認(rèn)為第一主成分表示相關(guān)性的整體一致性，能夠刻畫由進(jìn)化偏差引起的背景噪聲.第一特征值和對(duì)應(yīng)特征向量元素分別近似為

則背景噪音矩陣元素近似為

(28)

從上述分析可以看出SC關(guān)于背景噪音的近似與APC的平均乘積校正是等價(jià)的，都是秩為1的矩陣.

② LRS是SC和APC的擴(kuò)展和加強(qiáng)

LRS用低秩矩陣近似背景噪聲，是上述2種技術(shù)的擴(kuò)展；另外，用稀疏矩陣表征信號(hào)矩陣符合真實(shí)殘基接觸的稀疏性事實(shí).這從理論上保證了LRS方法的優(yōu)越性[36].

我們預(yù)期LRS將取代APC成為有效的去除背景噪聲的手段.

4 監(jiān)督學(xué)習(xí)方法

4.1 基本思想

1) 生物學(xué)觀察.殘基相互作用本身有一定的規(guī)律，蛋白質(zhì)殘基的性質(zhì)，如二級(jí)結(jié)構(gòu)、溶液可及性、疏水性等，對(duì)殘基間形成接觸有重要的作用.舉例來說，不同的二級(jí)結(jié)構(gòu)對(duì)于接觸的分布有重大影響，如Beta正平行和反平行片段之間的殘基接觸呈現(xiàn)出完全不同的模式.

2) 數(shù)學(xué)模型設(shè)計(jì).監(jiān)督學(xué)習(xí)方法將殘基間相互作用預(yù)測(cè)視為機(jī)器學(xué)習(xí)中的分類問題，首先對(duì)每個(gè)殘基對(duì)都提取多種特征(比如保守性、預(yù)測(cè)的二級(jí)結(jié)構(gòu)、溶液可及表面積等)，然后在已知?dú)埢嗷プ饔们闆r的集合上訓(xùn)練上述特征的權(quán)重.

4.2 典型方法

各類監(jiān)督學(xué)習(xí)方法的不同主要體現(xiàn)在2方面：1)選取的特征不同；2)采用的機(jī)器學(xué)習(xí)的模型不同.典型的方法包括整數(shù)規(guī)劃PhyCMAP[37]、隨機(jī)森林PconsC系列[38-39]和神經(jīng)網(wǎng)絡(luò)方法MetaPSICOV[11]，簡(jiǎn)要介紹如下.

4.2.1 典型方法1：整數(shù)線性規(guī)劃

Xu等人[37]考慮殘基性質(zhì)對(duì)殘基遠(yuǎn)程相互作用的重大影響，利用整數(shù)線性規(guī)劃將殘基相互作用需滿足的物理約束和共進(jìn)化信息整合起來，開發(fā)了軟件PhyCMAP.其基本思想是：1)采用隨機(jī)森林技術(shù)預(yù)測(cè)殘基間存在相互作用的概率；2)利用整數(shù)線性規(guī)劃選擇出概率較大的殘基間相互作用，同時(shí)要求這些殘基對(duì)滿足一些物理約束.

整數(shù)線性規(guī)劃的目標(biāo)函數(shù)為

(29)

整數(shù)線性規(guī)劃考慮以下物理約束：殘基i最多參與形成多少殘基相互作用；2個(gè)二級(jí)結(jié)構(gòu)之間最多會(huì)形成多少殘基相互作用；2個(gè)strand(形成sheet的片段單元)之間的相互作用具有連續(xù)性等.比如當(dāng)2個(gè)strand之間形成正平行sheet時(shí)，接觸的相鄰殘基對(duì)需要滿足以下約束：Yi,j≥Yi-1,j+Yi+1,j+1-1，其中，i,i±1表示其中一個(gè)strand上的殘基，j,j±1表示另一個(gè)strand上的殘基.該約束保證2個(gè)strand之間形成的相互作用具有連續(xù)性.

PhyCMAP由于同時(shí)考慮了真實(shí)殘基相互作用的限制和共進(jìn)化信息，其在CASP10和Set600數(shù)據(jù)集[37]上測(cè)試的結(jié)果表明PhyCMAP超過了當(dāng)時(shí)比較流行的其他監(jiān)督學(xué)習(xí)方法的軟件，例如NNcon[12],SVMcon[40]等.

同時(shí)PhyCMAP也有其局限性.其只考慮了局部模型MI輸出的共進(jìn)化信息，而并沒有考慮更加有效的全局模型輸出的共進(jìn)化信息.下面介紹的Pconsc和MetaPSICOV方法克服了PhyCMAP的局限性.

4.2.2 典型方法2：隨機(jī)森林方法PconsC系列

Skwark等人[39]發(fā)現(xiàn)不同的全局模型預(yù)測(cè)得到的殘基接觸集合有一定差異；而且不同的構(gòu)建多序列聯(lián)配的軟件輸出的多序列聯(lián)配也不同,這些不同的多序列聯(lián)配也會(huì)導(dǎo)致遠(yuǎn)程相互作用的預(yù)測(cè)結(jié)果不同.

基于上述觀察，2013年Skwark等人[39]將2種全局模型PSICOV和plmDCA對(duì)8種不同多序列聯(lián)配的預(yù)測(cè)結(jié)果與其他特征整合，提出了預(yù)測(cè)殘基間相互作用的隨機(jī)森林方法，并開發(fā)了軟件PconsC.多序列聯(lián)配由HHblits和jackhamme兩種比對(duì)軟件取定4種不同的閾值得到；其考慮的殘基對(duì)特征包括二級(jí)結(jié)構(gòu)預(yù)測(cè)值、殘基溶液可及表面積、殘基替代向量.實(shí)驗(yàn)結(jié)果表明：PconsC具有較高的預(yù)測(cè)精度，超過PSICOV和plmDCA的預(yù)測(cè)結(jié)果.

基于“相互作用殘基對(duì)的成簇性”這一認(rèn)識(shí)，Skwark等人[38]進(jìn)一步改進(jìn)PconsC，用多層隨機(jī)森林逐步過濾掉孤立的相互作用對(duì)，開發(fā)了軟件PconsC2.值得指出的是，PconsC2的另一個(gè)優(yōu)勢(shì)在于顯著減少了對(duì)樣本數(shù)的要求，從而首次實(shí)現(xiàn)當(dāng)同源序列少于1 000條時(shí)的準(zhǔn)確預(yù)測(cè).實(shí)驗(yàn)結(jié)果表明：PconsC2比已有的phyCMAP具有更優(yōu)的預(yù)測(cè)準(zhǔn)確率.

4.2.3 典型方法3：神經(jīng)網(wǎng)絡(luò)方法MetaPSICOV

如4.2.1節(jié)所述，PhyCMAP結(jié)合了結(jié)構(gòu)特征和共進(jìn)化信息，但是只是引入了局部的共進(jìn)化信息，并沒有引入更加有效的全局共進(jìn)化信息.Jones等綜合考慮結(jié)構(gòu)特征和全局的共進(jìn)化信息，提出了預(yù)測(cè)殘基間相互作用的神經(jīng)網(wǎng)絡(luò)模型，并開發(fā)了軟件metaPSICOV[11].

具體地，metaPSICOV是一個(gè)2層前向神經(jīng)網(wǎng)絡(luò)模型：第1層基于二級(jí)結(jié)構(gòu)、溶液可及性、殘基替代向量等特征以及PSICOV，mfDCA，plmDCA的預(yù)測(cè)結(jié)果，利用含55個(gè)隱單元的神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)出殘基相互作用概率的粗略估計(jì)；第2層以第1層的粗略估計(jì)為輸入特征，再加上部分結(jié)構(gòu)特征，使用相同的神經(jīng)網(wǎng)絡(luò)對(duì)相互作用概率的估計(jì)進(jìn)行校正.

metaPSICOV根據(jù)MSA質(zhì)量比較準(zhǔn)確地權(quán)衡共進(jìn)化特征和傳統(tǒng)特征(如二級(jí)結(jié)構(gòu)等)的權(quán)重，從而更加有效地整合多種信息，提高預(yù)測(cè)準(zhǔn)確率.結(jié)果顯示，metaPSICOV超過PSICOV，mfDCA，plmDCA的預(yù)測(cè)效果，并在第11屆CASP競(jìng)賽中獲第1名(詳細(xì)分析見第6節(jié)).

4.3 監(jiān)督學(xué)習(xí)方法的實(shí)驗(yàn)結(jié)果及分析

針對(duì)監(jiān)督學(xué)習(xí)方法的實(shí)驗(yàn)結(jié)果分析表明2點(diǎn)結(jié)論：

1) 早期的監(jiān)督學(xué)習(xí)方法，例如采用支持向量機(jī)模型的SVMcon[12]和采用神經(jīng)網(wǎng)絡(luò)模型的NNcon[40]等，由于沒有加入有效的共進(jìn)化特征，其效果并不比無監(jiān)督學(xué)習(xí)方法的全局模型好.

2) 近年來提出的監(jiān)督學(xué)習(xí)方法metaPSICOV和Pconsc2等，綜合了多種結(jié)構(gòu)特征和無監(jiān)督方法輸出的共進(jìn)化特征，其效果超過了無監(jiān)督學(xué)習(xí)方法和早期的監(jiān)督學(xué)習(xí)方法，從而表明將無監(jiān)督方法得到的結(jié)果整合到監(jiān)督學(xué)習(xí)方法中是當(dāng)前預(yù)測(cè)殘基相互作用最有效的策略.

5 現(xiàn)有軟件匯總

第2.4節(jié)所述的方法都有相應(yīng)的服務(wù)器為用戶提供殘基相互作用預(yù)測(cè)服務(wù)，我們將這些服務(wù)器匯總?cè)绫?所示：

Table 3 Overview of Existing Softwares for Protein Contact Prediction

6 CASP比賽中殘基相互作用預(yù)測(cè)性能分析

CASP競(jìng)賽是全球范圍內(nèi)的蛋白質(zhì)結(jié)構(gòu)預(yù)測(cè)比賽，現(xiàn)已作為客觀評(píng)估蛋白質(zhì)結(jié)構(gòu)預(yù)測(cè)質(zhì)量的標(biāo)準(zhǔn).從1994年開始，每?jī)赡?屆，迄今已舉辦11屆.目前，CASP競(jìng)賽包括結(jié)構(gòu)預(yù)測(cè)、殘基遠(yuǎn)程相互作用預(yù)測(cè)、接觸位點(diǎn)輔助結(jié)構(gòu)預(yù)測(cè)、結(jié)構(gòu)優(yōu)化、結(jié)構(gòu)質(zhì)量評(píng)估5個(gè)部分.

在蛋白質(zhì)結(jié)構(gòu)預(yù)測(cè)領(lǐng)域，大部分軟件采用開源軟件或者免費(fèi)預(yù)測(cè)服務(wù)的方式，商業(yè)軟件較少(比如DNAStar公司開發(fā)的NovaFold和BSI公司開發(fā)的RAPTOR等).國(guó)內(nèi)研究團(tuán)隊(duì)也多次參加CASP比賽，包括中科院生物物理研究所的Jiang-Server團(tuán)隊(duì)、上海交通大學(xué)的Shen-group團(tuán)隊(duì)以及中科院計(jì)算所的FALCON團(tuán)隊(duì).其中本課題團(tuán)隊(duì)開發(fā)的FALCON系列軟件在CASP-8中結(jié)構(gòu)預(yù)測(cè)FR-Hard類上獲得第3名，在CASP-11中結(jié)構(gòu)預(yù)測(cè)TBM類上獲得第9名.Shen-group在CASP-11殘基接觸預(yù)測(cè)的FM類蛋白上取得了第2名(以precision評(píng)價(jià)).

殘基遠(yuǎn)程相互作用預(yù)測(cè)作為CASP競(jìng)賽的重要部分.在2014年的第11屆CASP比賽中，共有29個(gè)軟件參加了殘基相互作用預(yù)測(cè)[41-45].本文提到的一些經(jīng)典算法參加了這次比賽，例如采用PhyCMAP方法的RaptorX-Contact軟件、采用MetaPSICOV的CONSIP2軟件以及采用PconsC2方法的MetaPSICOV軟件等.

在CASP11中，參賽軟件大多是基于監(jiān)督學(xué)習(xí)的方法.在監(jiān)督學(xué)習(xí)方法中，無監(jiān)督學(xué)習(xí)全局模型輸出的信息是其重要特征.根據(jù)是否含有全局共進(jìn)化模型信息，可以將這些監(jiān)督學(xué)習(xí)方法分為2類：1)不包含全局共進(jìn)化模型信息的方法：例如PhyCMAP，采用SVMcon方法的MULTICOM-construct、采用DNcon方法的MULTICOM-cluster軟件和采用NNcon方法的MULTICOM-novel軟件等；2)包含全局共進(jìn)化模型信息的方法：例如采用神經(jīng)網(wǎng)絡(luò)模型的CONSIP2方法、采用隨機(jī)森林模型的Pcons-net(PconsC2)、采用SVM模型的Shen-group和RBO_ALEPH[46]等.我們的實(shí)驗(yàn)結(jié)果分析表明包含全局共進(jìn)化模型信息的方法要好于不包含全局共進(jìn)化模型信息的方法.

CASP競(jìng)賽對(duì)于方法的衡量主要有4個(gè)分項(xiàng)：

1) 預(yù)測(cè)準(zhǔn)確率(precision)

根據(jù)各方法輸出的分?jǐn)?shù)進(jìn)行排序，選取TopN作為預(yù)測(cè)的接觸集合，其他的作為非接觸集合.其中N一般取L10,L5,L2或L等(L為目標(biāo)蛋白的序列長(zhǎng)度)，計(jì)算正陽性(true positive,TP)和假陽性(false positive,FP)，precision的定義為

2)Xd值

將氨基酸對(duì)之間的距離分成15個(gè)區(qū)間：(0,4]?，…，(56,60]?.Xd值的定義為

其中，Ppi表示預(yù)測(cè)的殘基接觸的距離在第i個(gè)區(qū)間的比例，Pai表示結(jié)構(gòu)所有的殘基接觸對(duì)的距離在i個(gè)區(qū)間的比例.Xd用來衡量真實(shí)結(jié)構(gòu)和預(yù)測(cè)接觸中殘基對(duì)距離分布的差異.

3) Matthews相關(guān)系數(shù)(MCC)

CASP比賽要求各參賽軟件給出每個(gè)殘基對(duì)接觸的概率.選取0.5作為閾值計(jì)算TP,TN，F(xiàn)P(false positive)和FN(false negative).

MCC=(TP×TN-FP×FN)

4) precision-recall曲線下的面積(AUC_PR)

由于上述共進(jìn)化分析方法給出相關(guān)性度量，而非殘基對(duì)的接觸概率，Matthews相關(guān)系數(shù)并不適用于這些方法，所以我們采用precision,Xd值和AUC_PR這3種分項(xiàng)對(duì)這些方法進(jìn)行評(píng)估.

我們?cè)贑ASP11比賽數(shù)據(jù)集合測(cè)試了本文綜述的典型方法以及參加CASP比賽的主要方法，對(duì)實(shí)驗(yàn)結(jié)果采用自行開發(fā)的程序進(jìn)行了詳細(xì)分析(程序和數(shù)據(jù)下載地址http:bioinfo.ict.ac.cnCOLORS)，詳細(xì)分析結(jié)果按照無模板建模(free modelling, FM)類和有模板建模(template based modelling， TBM)類分別描述如下：

6.1 在FM類目標(biāo)蛋白上的測(cè)試結(jié)果

在CASP11比賽中，共有17參賽隊(duì)伍提交了超過20個(gè)FM類蛋白質(zhì)域的預(yù)測(cè)結(jié)果.圖10列出了各方法的預(yù)測(cè)性能，包括precision,Xd和AUC_PR.從圖10我們可以得到如下結(jié)論.

1) 綜合全局共進(jìn)化模型信息的監(jiān)督方法領(lǐng)先于其他方法.比如CONSIP(metaPSICOV)在3種分項(xiàng)都取得第1名;Shen-group也處于領(lǐng)先地位.以precision和AUC_PR評(píng)價(jià)，RBO_ALEPH的排名也較靠前.

2) 監(jiān)督學(xué)習(xí)方法整體上優(yōu)于無監(jiān)督學(xué)習(xí)方法的效果.例如以precision作評(píng)價(jià)，無論是包含了全局共進(jìn)化信息的監(jiān)督方法(如CONSIP,Shen-group)還是不包含全局共進(jìn)化信息的監(jiān)督方法(如MULTICOM-novel,RaptorX,MULTICOM-cluster)，都超過了無監(jiān)督學(xué)習(xí)的方法(如plmDCA,PSICOV).

3) 在無監(jiān)督學(xué)習(xí)方法中，全局模型普遍優(yōu)于無監(jiān)督學(xué)習(xí)方法.

4) 以precision作為評(píng)價(jià)，LRS技術(shù)優(yōu)于APC技術(shù).

值得指出的是：不同預(yù)測(cè)方法在不同的評(píng)價(jià)指標(biāo)下表現(xiàn)不同，例如RaptorX以precision為評(píng)價(jià)排名為16，而以Xd為評(píng)價(jià)排名為2；對(duì)于無監(jiān)督學(xué)習(xí)方法的去噪音方法，以precision為評(píng)價(jià)，LRS好于APC；但以AUC_PR為評(píng)價(jià)，APC好于LRS.其原因在于：LRS只提取顯著的共進(jìn)化信號(hào)，將非顯著的共進(jìn)化信號(hào)的分?jǐn)?shù)設(shè)置為0；而APC卻可以對(duì)非顯著的共進(jìn)化信號(hào)進(jìn)行排名，從而造成如果以AUC_PR為評(píng)價(jià)，APC技術(shù)優(yōu)于LRS技術(shù).

Fig. 10 Prediction performance of different methods on CASP-11 FM targets圖10 典型方法對(duì)CASP-11 FM類蛋白的預(yù)測(cè)性能

6.2 在TBM類目標(biāo)蛋白上的測(cè)試結(jié)果

在CASP11比賽中，共有14個(gè)參賽組提交了多于60個(gè)TBM類別蛋白域的預(yù)測(cè)結(jié)果,如圖11所示:

Fig. 11 Prediction performance of different methods on CASP-11 TBM targets圖11 典型方法對(duì)CASP-11 TBM類蛋白的預(yù)測(cè)性能

我們可得到如下結(jié)論：

1) 與FM類目標(biāo)蛋白上的觀察相同，綜合了全局共進(jìn)化模型信息的監(jiān)督方法領(lǐng)先于其他方法.

2) 以precision為評(píng)價(jià)，整體上來講，無監(jiān)督學(xué)習(xí)方法中的全局模型優(yōu)于沒有結(jié)合全局模型信息的監(jiān)督學(xué)習(xí)方法.例如，plmDCA_LRS和PSICOV_LRS的預(yù)測(cè)性能優(yōu)于RaptorX和MULTICOM-cluster.這在很大程度上源于TBM目標(biāo)蛋白的多序列聯(lián)配比FM包含更多的同源序列，從而提供了更準(zhǔn)確的共進(jìn)化信息.

3) 與FM類目標(biāo)蛋白上的觀察相同，無監(jiān)督學(xué)習(xí)方法中的全局模型整體優(yōu)于局部模型.

6.3 預(yù)測(cè)準(zhǔn)確率與有效同源序列數(shù)目的關(guān)系

多序列聯(lián)配的有效同源序列的數(shù)目對(duì)預(yù)測(cè)準(zhǔn)確率的影響很大.無監(jiān)督學(xué)習(xí)方法中性能較好的全局模型往往具有較多的參數(shù)，從而需要大量的同源序列進(jìn)行參數(shù)估計(jì)；在監(jiān)督學(xué)習(xí)中使用的重要特征也受同源序列數(shù)目的影響，例如序列譜、預(yù)測(cè)的二級(jí)結(jié)構(gòu)和預(yù)測(cè)的可及水表面積等特征.

我們?cè)?9個(gè)FM類和TBM類目標(biāo)蛋白上測(cè)試同源序列數(shù)目對(duì)各方法預(yù)測(cè)準(zhǔn)確度的影響，選取經(jīng)典的無監(jiān)督學(xué)習(xí)方法和在CASP11中提交超過90個(gè)蛋白域的軟件進(jìn)行評(píng)價(jià).

我們根據(jù)Meff將這些目標(biāo)蛋白分成了3組：1)(0,100],共36個(gè)蛋白；2)(100,1 000]，共32個(gè)蛋白；3)Meff>1 000，共31個(gè)蛋白.

如圖12所示，我們可以得到以下結(jié)論：

1) 所有方法的預(yù)測(cè)準(zhǔn)確度都隨Meff的提高而提高，但不同方法的提高程度不同.

2) 總體來講，融合全局度量的監(jiān)督學(xué)習(xí)方法在3種Meff類別下都領(lǐng)先于其他方法，例如CONSIP2(MetaPSICOV),Pcons-net(Pcons2),RBO_Aleph等.

3) 當(dāng)Meff較低時(shí)，不融合全局度量的監(jiān)督學(xué)習(xí)方法優(yōu)于無監(jiān)督方法的全局模型.例如MULTICOM-CONSTRUCT，RaptorX-Conact，MULTICOM_CLUSTER優(yōu)于CCMpred，PSICOV等方法.

4 當(dāng)Meff較高時(shí)，全局模型逐漸超越不含全局模型信息的無監(jiān)督方法.

Fig. 12 Relationship of mean precision with Meff for different methods on CASP-11 targets(FM and TBM)圖12 典型方法對(duì)CASP-11 (FM和TBM)蛋白預(yù)測(cè)的平均準(zhǔn)確率與Meff的關(guān)系

7 殘基遠(yuǎn)程相互作用預(yù)測(cè)的應(yīng)用

殘基共進(jìn)化分析可在未知蛋白質(zhì)結(jié)構(gòu)時(shí)，僅依據(jù)序列推斷殘基間的相互作用，因此在蛋白質(zhì)結(jié)構(gòu)和功能預(yù)測(cè)上具有重要的研究意義.

1) 殘基間遠(yuǎn)程相互作用的信息能夠有效地促進(jìn)對(duì)蛋白質(zhì)結(jié)構(gòu)的預(yù)測(cè)，其典型工作是基于mfDCA的相互作用預(yù)測(cè)信息開發(fā)的蛋白質(zhì)結(jié)構(gòu)預(yù)測(cè)軟件DCA-fold[47]和EVfold[14].Hopf等人[48]考慮膜蛋白的結(jié)構(gòu)特性，使用EVfold對(duì)DrugBank數(shù)據(jù)庫中23個(gè)已知結(jié)構(gòu)的膜蛋白家族進(jìn)行結(jié)構(gòu)預(yù)測(cè)，發(fā)現(xiàn)其中20個(gè)可以得到高精度預(yù)測(cè)，進(jìn)而對(duì)Pfam中11個(gè)未知結(jié)構(gòu)蛋白質(zhì)進(jìn)行結(jié)構(gòu)預(yù)測(cè).

2) 由于蛋白質(zhì)結(jié)構(gòu)中的二硫鍵可以看成一種特殊的殘基相互作用，所以相互作用預(yù)測(cè)的信息也有助于二硫鍵的預(yù)測(cè).Yang等人提取遠(yuǎn)程相互作用預(yù)測(cè)軟件GREMLIN的輸出信息，作為支持向量回歸(support vector regression, SVR)模型的特征之一[49].結(jié)果表明其軟件Cyscon優(yōu)于其他同類軟件.

3) 殘基間相互作用能夠有助于對(duì)殘基功能的推斷.一般地，功能位點(diǎn)傾向于是蛋白質(zhì)中的保守位點(diǎn)；類似地，殘基間的關(guān)聯(lián)強(qiáng)度也有助于推斷功能位點(diǎn)，比如Marks等人[6,48]計(jì)算某特定殘基與其他殘基的累積耦合強(qiáng)度，作為該殘基功能性(功能選擇壓力)的度量，并發(fā)現(xiàn)了一些重要的功能元件.

8 相互作用預(yù)測(cè)方法發(fā)展趨勢(shì)分析

經(jīng)過20多年的發(fā)展，研究者已經(jīng)提出了多種殘基相互作用的預(yù)測(cè)模型，使得預(yù)測(cè)精度有了顯著提升.然而目前已有的算法依然存在假陽性率較高、對(duì)樣本數(shù)目顯著依賴等缺陷.在前期工作的基礎(chǔ)上，我們認(rèn)為殘基間相互作用有3個(gè)發(fā)展趨勢(shì)：

趨勢(shì)1. 改進(jìn)參數(shù)估計(jì)方法.經(jīng)典的統(tǒng)計(jì)模型如Markov隨機(jī)場(chǎng)能夠比較準(zhǔn)確地描述蛋白質(zhì)所有位點(diǎn)的全概率分布，但其參數(shù)估計(jì)的各種近似策略還有待改善，以進(jìn)一步縮小與極大似然估計(jì)之間的差異，并提高計(jì)算效率.

趨勢(shì)2. 預(yù)測(cè)二級(jí)結(jié)構(gòu)單元之間的相互作用.傳統(tǒng)遠(yuǎn)程相互作用預(yù)測(cè)方法都是預(yù)測(cè)殘基間相互作用，這在計(jì)算上具有很大便利，但是同時(shí)也造成顯著的假陽性.事實(shí)上，從預(yù)測(cè)蛋白質(zhì)整體結(jié)構(gòu)這個(gè)目標(biāo)來說，判斷二級(jí)結(jié)構(gòu)單元之間是否存在相互作用就能夠提供足夠有價(jià)值的信息，而不用細(xì)化到殘基間是否存在相互作用.

趨勢(shì)3. 改進(jìn)相互作用預(yù)測(cè)的評(píng)價(jià)方法.目前的評(píng)價(jià)方法中，所有殘基對(duì)是同等考慮的.然而在蛋白質(zhì)中，各個(gè)殘基的重要性是有顯著差異的，比如二級(jí)結(jié)構(gòu)單元之間的相互作用、疏水集團(tuán)與其他位點(diǎn)的相互作用等具有相對(duì)更高的重要性，而突變較多較隨機(jī)的不重要的位點(diǎn)則對(duì)結(jié)構(gòu)影響不大.一種可能的方案是首先基于多序列聯(lián)配給出位點(diǎn)重要性的先驗(yàn)概率，進(jìn)而在相互作用預(yù)測(cè)中有側(cè)重地考慮那些重要的殘基，這種有針對(duì)性地設(shè)置約束將能夠提高結(jié)構(gòu)預(yù)測(cè)效率和精度.

9 總 結(jié)

本文對(duì)殘基間遠(yuǎn)程相互作用預(yù)測(cè)進(jìn)行了綜述，主要介紹了殘基間相互作用預(yù)測(cè)的機(jī)器學(xué)習(xí)方法，分析了各方法的預(yù)測(cè)性能，并總結(jié)了未來的發(fā)展趨勢(shì).

值得指出的是，殘基間相互作用預(yù)測(cè)是機(jī)器學(xué)習(xí)中結(jié)構(gòu)學(xué)習(xí)(structured learning)的一個(gè)典型問題，因此這方面的研究不僅具有重要的生物學(xué)意義，同時(shí)能夠推動(dòng)機(jī)器學(xué)習(xí)領(lǐng)域的研究.

[1]Lodish H F, Berk A, Zipursky S L, et al. Molecular Cell Biology[M]. New York: WH Freeman, 2000

[2]Petsko G A, Ringe D. Protein Structure and Function[M]. London: New Science Press, 2004

[3]Wüthrich K. The way to NMR structures of proteins[J]. Nature Structural & Molecular Biology, 2001, 8(11): 923-925

[4]Kendrew J C, Bodo G, Dintzis H M, et al. A three-dimensional model of the myoglobin molecule obtained by X-ray analysis[J]. Nature, 1958, 181(4610): 662-666

[5]Taylor K A, Glaeser R M. Electron diffraction of frozen, hydrated protein crystals[J]. Science, 1974, 186(4168): 1036-1037

[6]Marks D S, Hopf T A, Sander C. Protein structure prediction from sequence variation[J]. Nature Biotechnology, 2012, 30(11): 1072-1080

[7]Anfinsen C B. Principles that govern the folding of protein chains[J]. Science, 1973, 181(4096): 223-230

[8]Kim De, Dimaio F, Wang R Y, et al. One contact for every twelve residues allows robust and accurate topology-level protein structure modeling[J]. Proteins: Structure, Function, and Bioinformatics, 2014, 82(S2): 208-218

[9]Haile J M. Molecular Dynamics Simulation[M]. New York: Wiley, 1992

[10]Anzai Y. Pattern Recognition and Machine Learning[M]. New York: Academic Press, 2012

[11]Jones D T, Singh T, Kosciolek T, et al. MetaPSICOV: Combining coevolution methods for accurate prediction of contacts and long range hydrogen bonding in proteins[J]. Bioinformatics, 2015, 31(7): 999-1006

[12]Cheng Jianlin, Baldi P. Improved residue contact prediction using support vector machines and a large feature set[J]. BMC Bioinformatics, 2007, 8(1): 11-13

[13]Chen Peng. Analysis and prediction of interactions between residues in proteins[D]. Hefei: University of Science and Technology of China, 2007 (in Chinese)(陳鵬. 蛋白質(zhì)殘基間的相互作用分析與預(yù)測(cè)[D]. 合肥: 中國(guó)科學(xué)技術(shù)大學(xué), 2007)

[14]Marks D S, Colwell L J, Sheridan R, et al. Protein 3D structure computed from evolutionary sequence variation[J]. PLoS ONE, 2011, 6(12): 1287-1296

[15]Gobel U, Sander C, Schneider R, et al. Correlated mutations and residue contacts in proteins[J]. Proteins: Structure, Function and Bioinfomatics, 1994, 18(4): 309-317

[16]Martin L C, Gloor G B, Dunn S D, et al. Using information theory to search for co-evolving residues in proteins[J]. Bioinformatics, 2005, 21(22): 4116-4124

[17]Kass I, Horovitz A. Mapping pathways of allosteric communication in GroEL by analysis of correlated mutations[J]. Proteins: Structure, Function, and Bioinformatics, 2002, 48(4): 611-617

[18]Kamisetty H, Ovchinnikov S, Baker D. Assessing the utility of coevolution-based residue-residue contact predictions in a sequence-and structure-rich era[J]. Proceedings of the National Academy of Sciences, 2013, 110(39): 15674-15679

[19]Morcos F, Pagnani A, Lunt B, et al. Direct-coupling analysis of residue coevolution captures native contacts across many protein families[J]. Proceedings of the National Academy of Sciences, 2011, 108(49): 1293-1301

[20]Weigt M, White R A, Szurmant H, et al. Identification of direct residue contacts in protein-protein interaction by message passing[J]. Proceedings of the National Academy of Sciences, 2009, 106(1): 67-72

[21]L?vkvist C, Lan Y, Weigt M, et al. Improved contact prediction in proteins: Using pseudolikelihoods to infer Potts models[J]. Physical Review E, 2013, 87(1): 12707-12929

[22]Burger L, van Nimwegen E. Disentangling direct from indirect co-evolution of residues in protein alignments[J]. PLoS Computational Biology, 2010, 6(1): 10006-10033

[23]Ma Jianzhu, Wang Sheng, Wang Zhiyong, et al. Protein contact prediction by integrating joint evolutionary coupling analysis and supervised learning[J]. Bioinformatics, 2015, 31(21): 3506-3513

[24]Jones D T, Buchan D W, Cozzetto D, et al. PSICOV: Precise structural contact prediction using sparse inverse covariance estimation on large multiple sequence alignments[J]. Bioinformatics, 2012, 28(2): 184-190

[25]Feizi S, Marbach D, Médard M, et al. Network deconvolution as a general method to distinguish direct dependencies in networks[J]. Nature biotechnology, 2013, 31(8): 726-733

[26]Meil? M, Jaakkola T. Tractable Bayesian learning of tree belief networks[C]Proc of the 6th Conf on Uncertainty in Artificial Intelligence. San Francisco, CA: Morgan Kaufmann, 2000: 380-388

[27]Lapedes A S, Bertrand G G, Liu L, et al. Correlated mutations in models of protein sequences: Phylogenetic and structural effects[J]. Lecture Notes-Monograph Series, 1999, 33(1), 236-256

[28]Plefka T. Convergence condition of the TAP equation for the infinite-ranged Ising spin glass model[J]. Journal of Physics A: Mathematical and general, 1982, 15(6): 1971-1985

[29]Georges A, Yedidia J S. How to expand around mean-field thery using high-temperature expansions[J]. Journal of Physics A: Mathematical and General, 1991, 24(9): 2173-2179

[30]Csisz X, R I, Talata Z. Consistent estimation of the basic neighborhood of Markov random fields[J]. The Annals of Statistics, 2006, 34(1): 123-145

[31]Lauritzen S L. Graphical Models[M]. Oxford, UK: Oxford University Press, 1996

[32]Wright S. Correlation and causation[J]. Journal of Agricultural Research, 1921, 20(7): 557-585

[33]Sun Haiping, Huang Yan, Wang Xiaofan, et al. Improving accuracy of protein contact prediction using balanced network deconvolution[J]. Proteins: Structure, Function, and Bioinformatics, 2015, 83(3): 485-496

[34]Dunn S D, Wahl L M, Gloor G B. Mutual information without the influence of phylogeny or entropy dramatically improves residue contact prediction[J]. Bioinformatics, 2008, 24(3): 333-340

[35]Halabi N, Rivoire O, Leibler S, et al. Protein sectors: Evolutionary units of three-dimensional structure[J]. Cell, 2009, 138(4): 774-786

[36]Zhang Haicang, Gao Yujuan, Deng Minghua, et al. Improving residue-residue contact prediction via low rank and sparse decomposition of residue correlation matrix[J]. Biochemical and Biophysical Research Communications, 2016, 472(1): 217-222

[37]Wang Zhiyong, Xu Jinbo. Predicting protein contact map using evolutionary and physical constraints by integer programming[J]. Bioinformatics, 2013, 29(13): 266-273

[38]Skwark M J, Raimondi D, Michel M, et al. Improved contact predictions using the recognition of protein like contact patterns[J]. PLoS Computatioal Biology, 2014, 10(11): 1003-1019

[39]Skwark M J, Abdel-Rehim A, Elofsson A. PconsC: Combination of direct information methods and alignments improves contact prediction[J]. Bioinformatics, 2013, 29(14): 1815-1816

[40]Tegge A N, Wang Z, Eickholt J, et al. NNcon: Improved protein contact map prediction using 2D-recursive neural networks[J]. Nucleic Acids Research, 2009, 37(Suppl 2): 515-518

[41]Kajan L, Hopf T A, Kalas M, et al. FreeContact: Fast and free software for protein contact prediction from residue co-evolution[J]. BMC Bioinformatics, 2014, 15(1): 158-164

[42]Seemayer S, Gruber M, S?ding J. CCMpred-fast and precise prediction of protein residue-residue contacts from correlated mutations[J]. Bioinformatics, 2014, 30(21): 3128-3130

[43]Baldassi C, Zamparo M, Feinauer C, et al. Fast and accurate multivariate Gaussian modeling of protein families: Predicting residue contacts and protein-interaction partners[J]. PLoS ONE, 2014, 9(3): 927-940

[44]Di Lena P, Nagata K, Baldi P. Deep architectures for protein contact map prediction[J]. Bioinformatics, 2012, 28(19): 2449-2457

[45]Monastyrskyy B, D’Andrea D, Fidelis K, et al. New encouraging developments in contact prediction: Assessment of the CASP11 results[J]. Proteins: Structure, Function, and Bioinformatics, 2015, 6(4): 126-140

[46]Schneider M, Brock O. Combining physicochemical and evolutionary information for protein contact prediction[J]. PloS ONE, 2014, 9(10): 1108-1120

[48]Hopf T A, Colwell L J, Sheridan R, et al. Three-dimensional structures of membrane proteins from genomic sequencing[J]. Cell, 2012, 149(7): 1607-1621

[49]Yang Jing, He Baoji, Jang R, et al. Accurate disulfide-bonding network predictions improve ab initio structure prediction of cysteine-rich proteins[J]. Bioinformatics, 2015, 31(23): 3773-3781

Zhang Haicang, born in 1987. PhD candidate from the Institute of Computing Technology, Chinese Academy of Sciences. His main research interests include bioinformatics, algorithm design and machine learning.

Gao Yujuan, born in 1992. PhD candidate from Peking University. Her main research interests include network inference and convex optimization algorithm.

Deng Minghua, born in 1969. Received his BS, MS, and PhD degrees in applied mathematics from Peking University. Professor in the School of Mathematical Sciences, Centre for Quantitative Biology and Center for Statistical Sciences, Peking University. His main research interests include bioinformatics and system biology.

Zheng Weimou, born in 1946. Received his BS degree in the Department of Physics from Peking University and PhD degree from Universite Libre de Bruxelles. Professor of the Institute of Theoretical Physics, Chinese Academy of Sciences. His main research interests include surface physics, stochastic process, nonlinear dynamics, biophysics and bioinformatics.

Bu Dongbo, born in 1973. Received his BS in computer science, MS and PhD degrees from the Institute of Computing Technology, Chinese Academy of Sciences. Professor of the Institute of Computing Technology, Chinese Academy of Sciences. Member of CCF. His main research interests include algorithm design and analysis, SAT problem, and bioinformatics (especially on genome sequencingassembly, protein sequencing via mass spectra, protein structure prediction).

A Survey on Algorithms for Protein Contact Prediction

Zhang Haicang1,2, Gao Yujuan3, Deng Minghua3,4,5, Zheng Weimou6, and Bu Dongbo11

(InstituteofComputingTechnology,ChineseAcademyofSciences,Beijing100190)2(UniversityofChineseAcademyofSciences,Beijing100049)3(CentreforQuantitativeBiology,PekingUniversity,Beijing100871)4(SchoolofMathematicalSciences,PekingUniversity,Beijing100871)5(CenterforStatisticalSciences,PekingUniversity,Beijing100871)6(InstituteofTheoreticalPhysics,ChineseAcademyofSciences,Beijing100190)

Proteins are large molecules consisting of a linear sequence of amino acids. In the natural environment, a protein spontaneously folds into specific tertiary structure to perform its biological functionality. The main factors that drive proteins to fold are interactions between residues, including hydrophobic interaction, Van der Waals’ force and electrostatic interactions. The interactions between residues usually lead to residue-residue contacts, and the prediction of residue-residue contacts should greatly facilitate understanding of protein structures and functionalities. A great variety of techniques have been proposed for residue-residue contacts prediction, including machine learning, statistical models, and linear programing. It should be pointed out that most of these techniques are based on the biological insight of co-evolution, i.e., during the evolutionary history of proteins, a residue’s mutation usually leads its contacting partner to mutate accordingly. In this review, we summarize the state-of-art algorithms in this field with emphasis on the construction of statistical models based on biological insights. We also present the evaluation of these algorithms using CASP (critical assessment of techniques for protein structure prediction) targets as well as popular benchmark datasets, and describe the trends in the field of protein contact prediction.

protein contact prediction; protein tertiary structure prediction; graphical model; co-evolution; machine learning

2015-12-10；

2016-04-14

國(guó)家“九七三”重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃基金項(xiàng)目(2012CB316502,2015CB910303)；國(guó)家自然科學(xué)基金項(xiàng)目(11175224,11121403,31270834,61272318,31171262,31428012,31471246)；中國(guó)科學(xué)院理論物理研究所理論物理國(guó)家重點(diǎn)實(shí)驗(yàn)室開放工程項(xiàng)目(Y4KF171CJ1) This work was supported by the National Basic Research Program of China (973 Program) (2012CB316502, 2015CB910303), the National Natural Science Foundation of China (11175224, 11121403, 31270834, 61272318, 31171262, 31428012, 31471246), and the Open Project Program of State Key Laboratory of the Institute of Theoretical Physics, Chinese Academy of Sciences (Y4KF171CJ1).

高玉娟(lacus2009@163.com，其對(duì)本文的貢獻(xiàn)同第一作者)；卜東波(dbu@ict.ac.cn)

TP399