深度學(xué)習(xí)在蛋白質(zhì)結(jié)構(gòu)預(yù)測中的應(yīng)用及啟示

位置特異性打分矩陣(position-specific scoring matrix, PSSM)或位置權(quán)重矩陣(position weight matrix, PWM)^[17]是蛋白質(zhì)及生物信息學(xué)里非常重要的統(tǒng)計量. 它主要衡量了不同氨基酸(或核酸)在蛋白質(zhì)(或DNA)上某個特定序列位置上出現(xiàn)的概率. 在一些機器學(xué)習(xí)預(yù)測蛋白質(zhì)的二級結(jié)構(gòu)類型時^[18~25]，常會將PSSM作為網(wǎng)絡(luò)的輸入. 但注意PSSM只包含殘基絕對位置屬性的信息，不包含不同殘基配對關(guān)聯(lián)信息.

圖1(a)以DNA為例，給出了統(tǒng)計PSSM矩陣的示意流程. 首先給定一個序列庫(例如針對基因庫的所有DNA數(shù)據(jù)，或蛋白質(zhì)庫里所有可能的序列)，圖中給出了由10個假想DNA序列組成的DNA庫；然后統(tǒng)計不同的核酸在特定位置出現(xiàn)的頻次矩陣(position frequency matrix, PFM)；再根據(jù)PFM得到位置概率矩陣(position probability matrix, PPM)；最后根據(jù)圖中公式算出位置權(quán)重矩陣PWM.

  

Fig. 1  Illustration of (a) position-specific scoring matrix (PSSM) and (b) multiple sequence alignment (MSA). (a) In this illustrative example, PSSM is computed using the formula given in the top right corner based on a DNA database consisting of ten DNA sequences. (b) Sequence alignment (SA) is trying to match the fragment pairs from the two given sequences as much as possible. In the alignment, inserting gaps "-" is allowed. Multiple sequence alignment (MSA) is SA on multiple sequences.

目前大多蛋白質(zhì)結(jié)構(gòu)預(yù)測的深度學(xué)習(xí)算法的輸入中都有多重序列比對信息(multiple sequence alignment, MSA)^{[12,15,26~39]}.

序列比對(sequence alignment)主要任務(wù)是針對查詢序列(query sequence)從數(shù)據(jù)庫中，用基因信息學(xué)的方法找到進化樹上盡可能同源的序列，然后根據(jù)變異的氨基酸的相似程度，按照特定規(guī)則來給該序列與查詢序列的相似度打分.

某個序列的變異包括對序列中特定片段的插入、刪除和替換. 相對于查詢序列，當庫里的蛋白質(zhì)序列變異很少時，則兩者相似度高.

當變異多時，還需根據(jù)進化同源的特點分類對變異片段進行進一步分析. 變異的氨基酸片段可分為保守片段(功能及化學(xué)特性相同)、半保守片段(功能及化學(xué)特性相近)和非保守片段(化學(xué)特性相差甚遠). 顯然，若保守片段越多，表明與查詢序列越接近.

比對的目標是通過恰當?shù)夭迦肟掌?gap)，使得插入空片段后的2個序列盡量相似(如圖1(b)左圖所示). 比對的方法有許多^[27]，例如動態(tài)規(guī)劃(dynamic programming)和點陣法(dot-matrix method).

用上述比對方法對若干個給定的序列與查詢序列進行比對就稱為多重序列比對(multiple sequence alignment, MSA). 通?？捎密浖﨏lustalW, MAFFT, ClustalOmega以及MUSCLE等算法程序?qū)Χ鄠€序列進行MSA比對^[40~52].

而在蛋白質(zhì)預(yù)測中，通常會針對輸入的蛋白質(zhì)序列，從蛋白質(zhì)數(shù)據(jù)庫中找到與給定序列相近的若干個序列，然后再將這些MSA作為神經(jīng)網(wǎng)絡(luò)的輸入. 此信息相比于PSSM包含了更為豐富的信息. 可從MSA中看出目標序列大致從哪些序列變異而來. 在深度學(xué)習(xí)中，MSA數(shù)據(jù)維度為(N_seq, N_res, 21)，其中N_seq為MSA包含序列的數(shù)目，N_res為目標序列的長度，21用于分辨20種氨基酸和gap“-”的熱點表征(有時可能為22或23).

如圖2所示，圖2(b)與2(c)是一個HP蛋白質(zhì)模型結(jié)構(gòu)^[39]的接觸圖(Contact Map)與距離圖(Distogram). 其中接觸圖中只有2個殘基接觸時，才有值(黑)；而Distogram灰度值對應(yīng)于兩殘基的距離，當距離大于截斷閾值時，灰度為0(白色).

  

Fig. 2  Illustration of contact map and distogram. (a) A typical structure of a given HP protein. (b) The contact of the (c) structure where the black square indicates the matrix element corresponding to two contact residues. (c) The distogram of the (a) structure where the greyness indicates the distance between two residues.

同一序列中不同殘基間的接觸與否或距離是非常重要的信息，它基本蘊含了蛋白質(zhì)骨架的三維結(jié)構(gòu)所有的信息. 而且這個信息相比于純粹的結(jié)構(gòu)三維坐標信息有2個優(yōu)勢：(1) 具有旋轉(zhuǎn)平移不變性，而三維坐標會隨著蛋白質(zhì)的旋轉(zhuǎn)或平移而改變；(2) 表達更簡潔及更易標準化. 因為存在關(guān)聯(lián)變異(correlated mutation)現(xiàn)象，有些接觸的兩氨基酸會同時變異以保證變異后仍接觸，故接觸圖或距離圖信息就顯得相當重要^[35,36,39].

基于上述原因，在最近的深度學(xué)習(xí)預(yù)測蛋白質(zhì)結(jié)構(gòu)的實踐中^{[12,15,53~62]}，大多都會采用此信息去提高預(yù)測準確性或預(yù)測給定蛋白質(zhì)的Contact Map或Distogram.

目前最著名的蛋白質(zhì)數(shù)據(jù)庫為PDB^[63],即Protein Data Bank，收藏了約1×10⁵多條蛋白質(zhì)的三維結(jié)構(gòu)數(shù)據(jù). 這些結(jié)構(gòu)由X射線、NMR或電子顯微鏡等方法獲得.

Critical Assessment of Protein Structure Prediction (CASP)^[64,65]是蛋白質(zhì)結(jié)構(gòu)預(yù)測科學(xué)共同體舉辦的兩年一次的競賽，每次競賽優(yōu)勝者的水平基本代表了當前世界結(jié)構(gòu)預(yù)測的最高水準(benchmark progress). 在每次競賽中，舉辦方會給出若干個已知結(jié)構(gòu)但未曾公開的蛋白質(zhì)序列，參賽團隊在規(guī)定時間內(nèi)提交各自的結(jié)構(gòu)預(yù)測結(jié)果，同時不限制預(yù)測方法. CASP自1994年以來共舉辦了14屆，其中最近2屆的第一名皆來自deepmind的AlphaFold算法.

之前，人們通常用距離均方差root mean squared deviation (RMSD)衡量2個分子構(gòu)象的接近程度. 但現(xiàn)在模版建模得分template modelling score被認為是更準確的衡量方式^[66]. 其表達式如下：

其中

式中n為蛋白質(zhì)的殘基數(shù)，M為旋轉(zhuǎn)平移矩陣. 上式表達的含義是將預(yù)測得到的結(jié)構(gòu)與各種旋轉(zhuǎn)平移操作后的真實結(jié)構(gòu)進行比較，取最相近(極大)的那個作為最后的分值.

顯然TM score在0~1之間，分數(shù)越高表明越準確. 通常認為當TM>0.5時，預(yù)測與真實之間的折疊基本一致^[56]；而對同一蛋白質(zhì)，NMR與X射線測出結(jié)構(gòu)之間的TM分數(shù)為0.807 ± 0.107左右. 所以，可認為當TM分數(shù)>0.8時，預(yù)測的結(jié)果已經(jīng)完全正確.

而AlphaFold2(AF2)近2/3的預(yù)測結(jié)果達到中低分辨率的實驗精度^[12]. 也即AF2幾乎解決了單域蛋白質(zhì)折疊預(yù)測問題^[14].

由于多域蛋白質(zhì)各功能域之間可以相對獨立地移動旋轉(zhuǎn)，在評估多域蛋白質(zhì)結(jié)構(gòu)相似性上，局域距離差異性測試(local distance difference test)是一個比TM分數(shù)更佳的評分方式. lDDT不同于TM，不依賴于骨架α碳原子的重疊，能夠不受功能域間位移的影響，更加有效地評估結(jié)構(gòu)之間的局域相似性^[67].

本小節(jié)僅羅列結(jié)構(gòu)預(yù)測涉及到的深度學(xué)習(xí)技術(shù)及原理, 具體請參考相關(guān)文獻^[68].

神經(jīng)網(wǎng)絡(luò)(neural network)可抽象成一個函數(shù)y=f(x;w)，它關(guān)聯(lián)了2組信息數(shù)據(jù)x與y(比如蛋白質(zhì)的序列x及其結(jié)構(gòu)y)，分別稱為網(wǎng)絡(luò)輸入與輸出；w為網(wǎng)絡(luò)的參數(shù). 神經(jīng)網(wǎng)絡(luò)訓(xùn)練的目標是為了找到恰當?shù)?em data-v-838f3892="" class="italic" style="padding: 0px 0.14em 0px 0px; margin: 0px;">w使得網(wǎng)絡(luò)能夠根據(jù)x準確地預(yù)測y.

普適近似原理(universal approximation theorem)^[69]表明單隱藏層的神經(jīng)網(wǎng)絡(luò)，只要其激活函數(shù)為非線性且神經(jīng)元數(shù)目足夠多，便可無限精確近似任意非線性映射. 普適近似原理表明NN可用于擬合任意未知關(guān)聯(lián).

神經(jīng)網(wǎng)絡(luò)設(shè)計要點：考察待預(yù)測的量y與哪些量有關(guān)聯(lián)，即找出哪些信息可足夠推導(dǎo)出y，然后將這些信息與y之間架接合適的神經(jīng)網(wǎng)絡(luò)便可. 信息間的關(guān)聯(lián)如果能用現(xiàn)有知識進行關(guān)聯(lián)就用現(xiàn)有知識將其關(guān)聯(lián)；未知關(guān)聯(lián)用神經(jīng)網(wǎng)絡(luò)代替.

神經(jīng)網(wǎng)絡(luò)選擇需要考慮輸入輸出信息數(shù)據(jù)特點，目前結(jié)構(gòu)預(yù)測中常用的網(wǎng)絡(luò)結(jié)構(gòu)主要有下面幾種.

殘差網(wǎng)絡(luò)(resnet)^[70]的基本思想是不斷地將未處理過的信息直接復(fù)制并疊加到下面幾層由網(wǎng)絡(luò)抽取出的特征上去. 殘差網(wǎng)絡(luò)于2015年提出，后來被廣泛運用于圖像處理中.

基于自注意力機制的transformer^[28]近幾年備受人工智能領(lǐng)域喜愛，它幾乎完全取代循環(huán)神經(jīng)網(wǎng)絡(luò)^[68]，其基本思想是從不同位置對之間提取信息，適合處理文本類、時序性的信息，不過近年也常用于圖像處理. AlphaFold2^[12]中大量使用了自注意力機制.

TBM方法^[14,29~38]大致步驟如圖3所示，通常可分為以下幾步. 第一步，通過數(shù)據(jù)庫檢索，得到目標蛋白質(zhì)的一組同源性序列(MSA)，并根據(jù)MSA獲得1個或多個折疊結(jié)構(gòu)模板. 第二步，比對目標序列和模板對應(yīng)序列，兩序列一致的片段直接使用模板的對應(yīng)折疊結(jié)構(gòu). 第三步，對于目標序列與模板對應(yīng)序列不一致的區(qū)域，采用碎片組裝、優(yōu)化算法或是數(shù)據(jù)庫方法等單獨預(yù)測. 當然通常最后還會用諸如分子動力學(xué)的優(yōu)化方法進行模型的精細化(model refinement)，以優(yōu)化全局結(jié)構(gòu)^[13,14]. 歷史上，TBM方法^[14]可以細分成comparative modelling (CM)和threading 2種方法^[29~38]. 其中在CM中，模板與目標序列的同源性較近.

  

Fig. 3  Illustration of basic strategies of template-based method (TBM) and template-free method (FM).

無模板方法^[72~77]的流程見圖3右半側(cè)，F(xiàn)M從蛋白質(zhì)數(shù)據(jù)庫中依MSA比對結(jié)果找到一些片段的結(jié)構(gòu)并將其放入片段庫中，然后找到評分較高的片段結(jié)構(gòu)拼成初始結(jié)構(gòu)^[72~74]，接著采用FM里非常重要的片段組裝(fragment assembly)方法^[72]，大致凍結(jié)片段的結(jié)構(gòu)并以片段結(jié)構(gòu)為單元來演化全局結(jié)構(gòu)，比如可根據(jù)粗粒化的勢能函數(shù)用梯度下降化進行能量優(yōu)化.

人們發(fā)現(xiàn)在蛋白質(zhì)變異過程中經(jīng)常出現(xiàn)關(guān)聯(lián)變異(correlated mutation)的現(xiàn)象：一條蛋白質(zhì)鏈內(nèi)若發(fā)生變異，總是2個氨基酸成對地變異；因為演化壓力會迫使蛋白質(zhì)維持一致構(gòu)型，原本接觸的氨基酸對在變異過程中繼續(xù)保持接觸，可以避免其形狀發(fā)生劇烈變化. 因此，這就使得殘基接觸對(inter-residue contact map)的信息極為重要^[56~61].

早期有許多傳統(tǒng)方法致力于預(yù)測殘基接觸對. 處理該問題的早期算法，傾向于以一次一對的形式、孤立地預(yù)測每個接觸對是否可能. 由于忽視了蛋白質(zhì)包含的全局信息：一個殘基對是否接觸受到序列中其他殘基的影響，早期算法陷入了困境，預(yù)測效果糟糕. 而之后研究者提出了充分利用全局信息的預(yù)測方法，例如基于Markov隨機場模型MRF的direct coupling method (DCA)^[58~61]，在殘基接觸預(yù)測上獲得了突破性的成就.

深度神經(jīng)網(wǎng)絡(luò)在預(yù)測殘基接觸對問題上，也表現(xiàn)出了異常優(yōu)異的性能，有時甚至還直接被用于預(yù)測鍵角等信息. 這些預(yù)測特征均可作為約束，輔助指導(dǎo)無模板方法.

比如，Raptor X-Contact深度學(xué)習(xí)模型^[39]將Contact Map的預(yù)測當成圖片分割任務(wù)來對待，Raptor X-Contact所采用的方法也被其他方法，如ResPRE^[54]所采納. ResPRE^[54]采用了圖片識別領(lǐng)域非常著名的殘差網(wǎng)絡(luò)(Resnet)模塊^[70]，殘差網(wǎng)絡(luò)的重要思想是不斷地將網(wǎng)絡(luò)前面的信息直接復(fù)制到網(wǎng)絡(luò)后面.

而AlphaFold1^[15]又將Contact Map拓展成距離直方圖(distogram)預(yù)測，基于此，它在2018年CASP13的比賽中獲得了巨大成功.

2020年的CASP14的比賽中，AlphaFold2 (AF2)^[12]取得了驕人的成績. 對來自89個域(domain)實驗測得的蛋白質(zhì)結(jié)構(gòu)，AlphFold2在88個域TM分數(shù)>0.5，59個域分數(shù)>0.914. 前者意味著預(yù)測結(jié)果與答案之間折疊基本一致. NMR、X射線晶體學(xué)測出的一組112個單域蛋白質(zhì)，序列相同率大于95%. NMR與X射線測出的結(jié)構(gòu)之間的TM值為0.807±0.107. 這說明AlphFold2的近60%的預(yù)測達到中低分辨率的實驗精度. 也就是說AlphFold2幾乎解決了單域蛋白質(zhì)折疊預(yù)測問題^[14].

AlphaFold2深度學(xué)習(xí)模型的結(jié)構(gòu)簡圖如圖4所示，具體參考文獻^[12]. 它分別借助了基因同源信息和蛋白質(zhì)結(jié)構(gòu)數(shù)據(jù)庫模板信息. 如圖所示，根據(jù)同源信息，可得到序列比對信息MSA，通過同源搜索得到與輸入序列同源相近的(s-1)條序列和輸入序列一起放到MSA數(shù)組里，再通過線性神經(jīng)網(wǎng)絡(luò)變換得到MSA表征，此表征的維度為(s，r，c)，其中r為蛋白質(zhì)的序列長度，c為表征的特征數(shù)(通道數(shù)). MSA表征包含了輸入序列與其他同源序列間的關(guān)系.

  

Fig. 4  Sketch of the AlphaFold2 model. Detailed description is referred to Ref.[12].

而另一輸入通道中，主要輸入與MSA相對應(yīng)的序列的結(jié)構(gòu)殘基對距離信息以及扭轉(zhuǎn)角的信息. 在具體輸入時，AF2將距離對長度劃分成64個離散塊(64 bins)，并將其轉(zhuǎn)化為概率的形式，故對應(yīng)數(shù)組形狀為(s，r，r，64)，取值為0~1. 注意配對表征中，只包含了MSA除輸入序列之外的某個序列自己結(jié)構(gòu)信息，不同序列之間并沒有進行信息的關(guān)聯(lián).

然后再將MSA表征與配對表征輸入一個稱為Evoformer的模塊，此模塊主要將MSA的信息(同源性差異)與結(jié)構(gòu)信息整合起來，最后得到輸入序列的MSA表征與輸入序列的配對表征. 此時，輸入序列的配對表征同時將演化信息與其他模板結(jié)構(gòu)信息有機地融合在了一起. Evoformer主要利用了自注意力機制來實現(xiàn)上述信息整合.

而下一個結(jié)構(gòu)模塊structure module主要的功能是將Evoformer預(yù)測的配對表征展開成三維空間結(jié)構(gòu)，同時亦承擔(dān)一定的預(yù)測調(diào)整功能. 此模塊的結(jié)構(gòu)大致如圖5所示. 一條蛋白質(zhì)骨架結(jié)構(gòu)可想象成一系列三角形的疊加，三角形的中心相當各個殘基α碳的坐標，三角形平面本身代表N-α-C-C構(gòu)成的三角形. 這樣，此骨架可由2個數(shù)組表示，數(shù)組形狀分別為(r，3×3)和(r，3)，分別表示每個三角形取向與位置.

  

Fig. 5  Illustration of how the pairing information is transformed into the 3D structure using neural networks in AlphaFold2^[12].

初始時，假設(shè)所有氨基酸都在原點，然后將此初始骨架與配對表征輸入結(jié)構(gòu)模塊，由于配對表征存有距離對及取向信息，故可通過一個稱為不變點注意力神經(jīng)網(wǎng)絡(luò)模塊將其初步還原成展開的骨架結(jié)構(gòu)，緊接著再加入側(cè)鏈原子從而得到全原子的三維結(jié)構(gòu).

如圖4所示，最后再將中間輸出的MSA信息、配對信息和3D結(jié)構(gòu)信息重新疊加輸入到Evoformer，如此反復(fù)迭代3次，最終到預(yù)測結(jié)果.

因為PDB中只有大約1×10⁵多個的序列有對應(yīng)的三維結(jié)構(gòu)數(shù)據(jù). 而在big fantastic database (BFD)蛋白質(zhì)序列數(shù)據(jù)有多達2,204,359,010個序列，雖然這些序列并不一定有對應(yīng)的三維結(jié)構(gòu)信息(無標簽)，但self-distillation dataset的訓(xùn)練技巧可以將這些無答案的題目作為作業(yè)進行訓(xùn)練，自己提高預(yù)測準確度，AlphaFold2用此擴大訓(xùn)練集并進一步提高了預(yù)測準確度.

后來有諸多研究團隊對AlphaFold2進行了拓展與提升. 例如：Baker團隊^[78]的RoseTTAFold發(fā)展了三通路神經(jīng)網(wǎng)絡(luò)(three-track neural network)，對AlphaFold2只包括1D序列信息和2D距離圖信息的兩通路神經(jīng)網(wǎng)絡(luò)模型進行了拓展，引入了3D結(jié)構(gòu)通路道網(wǎng)絡(luò)模塊；高毅勤團隊的MindSpore算法^[79]對AlphaFold2的計算速度進行了較大的提升.

真實蛋白質(zhì)結(jié)構(gòu)預(yù)測無論從訓(xùn)練數(shù)據(jù)準備還是模型構(gòu)建及訓(xùn)練都極其復(fù)雜. 因此，人們希望找一個簡單的蛋白質(zhì)模型，以便能快速地試驗他們的想法. 就如手寫數(shù)字識別(對應(yīng)數(shù)據(jù)集為MNIST)^[80]對于圖像識別一樣，所有的方法都會用MNIST數(shù)據(jù)集先來檢驗其有效性.

而HP蛋白質(zhì)模型就是這樣的模型^[2,6,8,9]. 它僅有2類氨基酸H和P，其中H代表厭水型氨基酸，P代表親水型.

我們基于此HP模型，提出了一個強關(guān)聯(lián)神經(jīng)網(wǎng)絡(luò)^[16]，如圖6所示，此神經(jīng)網(wǎng)絡(luò)有2個核心要素，一是不同于傳統(tǒng)的向量表征，它采用一個小的神經(jīng)網(wǎng)絡(luò)來代表每個氨基酸，每種氨基酸都用一個神經(jīng)網(wǎng)絡(luò)來表征，不同氨基酸對應(yīng)的網(wǎng)絡(luò)的權(quán)重亦不同，而相同的氨基酸共享網(wǎng)絡(luò)權(quán)重；二是它有一個自洽循環(huán)通路，這樣可使得輸出的信息(環(huán)境)與氨基酸的屬性發(fā)生強關(guān)聯(lián).

  

Fig. 6  Architecture of the strongly-correlated neural network (SCN) where r indicates number of residues and c indicates number of features or channels.

該研究發(fā)現(xiàn)與傳統(tǒng)向量表征方法相比，強關(guān)聯(lián)網(wǎng)絡(luò)極大提升了預(yù)測準確性，提高了約20個百分點.