馮子健，黃偉鴻，姜博文

（浙江理工大學(xué) 信息科學(xué)與工程學(xué)院，浙江 杭州 310018）

0 引言

蛋白質(zhì)是一種由多個氨基酸分子組成的鏈，它是生物界的重要組成物質(zhì)，參與了大多數(shù)生命活動。預(yù)測、理解和分析蛋白質(zhì)對醫(yī)學(xué)、遺傳學(xué)、藥學(xué)等領(lǐng)域具有重要意義［1］。隨著新一代測序技術(shù)的發(fā)展，研究者能很容易獲得蛋白質(zhì)的氨基酸序列，然而蛋白質(zhì)的功能卻主要由其三維結(jié)構(gòu)決定。

傳統(tǒng)生化實驗方法獲取蛋白質(zhì)結(jié)構(gòu)既昂貴又耗時，為了能快速、高效獲得蛋白質(zhì)三維結(jié)構(gòu)，計算方法便成為了一種重要的補充手段，在近年來得到了快速發(fā)展。然而，目前仍未存在一種方法的預(yù)測結(jié)果能完全達到生化實驗的準確性，只能在特定環(huán)境下達到理想效果，因此評估計算方法便輸出的蛋白質(zhì)模型成為了實驗中不可缺少的一環(huán)。

質(zhì)量評估（Quality Assessment，QA）的目的就是將預(yù)測的蛋白質(zhì)模型與天然結(jié)構(gòu)進行比較，進而輸出預(yù)測模型在各種評價指標上的得分，便于后期篩選出更接近天然結(jié)構(gòu)的優(yōu)秀模型。

AlphaFold2［2］在CASP14 中的表現(xiàn)標志著深度學(xué)習(xí)正顛覆性地改變蛋白質(zhì)結(jié)構(gòu)預(yù)測領(lǐng)域，但其在質(zhì)量評估領(lǐng)域的應(yīng)用還有待開發(fā)。為此，本文提出一種基于動態(tài)圖卷積神經(jīng)網(wǎng)絡(luò)［3］（Dynamic Graph CNN，DGCNN）的蛋白質(zhì)模型質(zhì)量評估方法DGCQA，并結(jié)合遷移學(xué)習(xí)思想，選用預(yù)訓(xùn)練模型ESM-1b［4］的編碼特征作為輸入。在CASP13［5］數(shù)據(jù)集上的實驗表明，所提模型相較于同類方法效果更好。

1 相關(guān)工作

蛋白質(zhì)模型質(zhì)量評估方法大致可分為共識方法和單模型方法。共識方法通過一個候選池中的其他模型信息來評估蛋白質(zhì)模型，它是質(zhì)量評估發(fā)展初期的主流方法之一，盡管共識方法可在一些情況下取得較好的效果，但很大程度上受候選池中模型的影響，如果模型較少或缺少一致性和相似性時，將難以對蛋白質(zhì)模型作出最正確的評價。單模型方法是將單個蛋白質(zhì)結(jié)構(gòu)作為輸入，首先提取結(jié)構(gòu)中的特征信息，然后通過機器學(xué)習(xí)、深度學(xué)習(xí)等方法進行訓(xùn)練。

從CASP13 開始，單一模型方法在質(zhì)量評估領(lǐng)域開始有逐漸超越共識方法的趨勢。ProQ2［6］、ProQ3［7］是質(zhì)量評估中機器學(xué)習(xí)的典型方法，使用人工微調(diào)的特征來訓(xùn)練支持向量機（Support Vector Machine，SVM），但由于機器學(xué)習(xí)的局限性，該類方法難以學(xué)習(xí)蛋白質(zhì)結(jié)構(gòu)中更為復(fù)雜的特征關(guān)系。

近年來，伴隨著深度學(xué)習(xí)的快速發(fā)展，一些基于深度學(xué)習(xí)的方法開始逐漸涌現(xiàn)。Hou 等［8］采用多任務(wù)學(xué)習(xí)訓(xùn)練了一個1DCNN 網(wǎng)絡(luò)來評估蛋白質(zhì)的全局和局部質(zhì)量，但蛋白質(zhì)存在多個殘基，難以捕獲殘基間的長距離依賴關(guān)系。DeepAccNet［9］通過結(jié)合2D 和3D 卷積來預(yù)測每個殘基的準確度，其中3D 卷積用于評估當(dāng)前蛋白質(zhì)的局部原子環(huán)境，2D 卷積則提供全局上下文信息。Nie 等［10］提出一種基于多尺度卷積和雙向門控循環(huán)神經(jīng)網(wǎng)絡(luò)的質(zhì)量評估方法BMBQA，結(jié)合數(shù)據(jù)增強思想，在原數(shù)據(jù)集中加入大量同源蛋白質(zhì)來提升模型性能，但數(shù)據(jù)增強不僅增加了訓(xùn)練負擔(dān)，還可能引入額外噪聲。

圖數(shù)據(jù)結(jié)構(gòu)能方便、直觀地表示蛋白質(zhì)中的原子、氨基酸及他們之間的相互作用關(guān)系。ProteinGCN［11］、GraphQA［12］將蛋白質(zhì)模型表示為圖形式，通過圖卷積神經(jīng)網(wǎng)絡(luò)（Graph Convolutional Network，GCN）［13］提取蛋白中殘基的特征，再利用這些特征評估蛋白質(zhì)模型質(zhì)量。雖然，該方法效果較好，但圖神經(jīng)網(wǎng)絡(luò)在使用前需要構(gòu)建蛋白質(zhì)的圖結(jié)構(gòu)，而圖中的邊往往需要根據(jù)殘基間的距離進行構(gòu)建，如果構(gòu)建邊的閾值選取不合理或蛋白質(zhì)模型中的數(shù)據(jù)不夠準確，將造成圖中產(chǎn)生錯誤的邊或邊缺失等問題，會在一定程度上影響GCN 的效果。

為此，本文總結(jié)前人工作經(jīng)驗，針對圖網(wǎng)絡(luò)在蛋白質(zhì)質(zhì)量評估中存在的問題，提出基于動態(tài)圖卷積神經(jīng)網(wǎng)絡(luò)的單一模型質(zhì)量評估方法，并結(jié)合預(yù)訓(xùn)練模型ESM-1b 編碼特征進一步提升模型效果。

2 實驗方法

2.1 輸入特征

本文輸入為多種特征組合，可將其劃分為殘基特征和殘基之間的成對特征。為了方便表述，設(shè)輸入的蛋白質(zhì)中殘基個數(shù)為L，ri表示蛋白質(zhì)序列中索引為i的殘基，具體的殘基特征如下。

2.1.1 one-hot編碼

one-hot 編碼用來表示蛋白質(zhì)中殘基類型的二進制向量，其中一位為1，其余位置為0。本文工作中，每一個殘基被編碼為一個長度為21 的向量（包括20 種標準氨基酸和其他氨基酸），蛋白質(zhì)one-hot編碼大小為L× 21。

2.1.2 殘基相對位置編碼

該特征向量用來表示每個殘基在它所在蛋白質(zhì)序列中的相對位置，計算公式如式（1）所示。

2.1.3 三維坐標

本文從PDB 文件中提取所有原子的三維坐標(x，y，z)，并用Cβ原子的坐標代表殘基坐標。

2.1.4 統(tǒng)計量

本文參考Hurtado 等［14］的工作，從蛋白質(zhì)多序列比對（Multiple Sequence Alignment，MSA［15］）中提取自信息量（self-information）和部分熵（partial entropy），如式（2）、式（3）所示。

式中：pi表示殘基ri出現(xiàn)在當(dāng)前位置的頻率表示數(shù)據(jù)集的平均頻率，兩種統(tǒng)計量均為23維向量。

2.1.5 DSSP特征

DSSP 特征為DSSP［16］計算得到的二面角（dihedral angle）、相對可溶性（relative solvent accessibility）和二級結(jié)構(gòu)（secondary structure）類型，并他們拼接為一個6維向量。

2.1.6 ESM-1b編碼

本文受到遷移學(xué)習(xí)思想啟發(fā)，將Facebook AI Research提出的蛋白質(zhì)預(yù)訓(xùn)練模型ESM-1b 作為固定的特征編碼器，將蛋白質(zhì)序列編碼為一個L× 1 280 的向量，并在輸入前通過一維卷積降維到32 維，與其他特征拼接作為最終輸入。ESM-1b 是一個基于Transformer［17］的無監(jiān)督預(yù)訓(xùn)練模型，在具有上億個蛋白質(zhì)序列的數(shù)據(jù)庫UniProt［18］上進行訓(xùn)練，可挖掘蛋白質(zhì)序列中的結(jié)構(gòu)信息。

綜上，殘基對特征的殘基間序列間隔表示蛋白質(zhì)序列中的兩兩殘基間的索引距離（即兩殘基的索引差），用onehot 編碼表示距離間隔［1，2，3，4，5，6-10，11-15，16-20，＞20］，大小為L×L× 9。殘基間方位包括由trRosetta［19］定義的ω、θ二面角及φ平面角。

2.2 基于動態(tài)圖卷積的蛋白質(zhì)模型質(zhì)量評估方法

2.2.1 Inception模塊

本文選用Szegedy 等［20］提出的Inception 結(jié)構(gòu)提取殘基對特征，模型結(jié)果如圖1 所示。由此可見，該結(jié)構(gòu)為多尺度卷積核組成的網(wǎng)絡(luò)模塊，模塊中3×3 和5×5 卷積核可讓網(wǎng)絡(luò)捕獲不同大小鄰域內(nèi)的殘基對信息，為下游任務(wù)提供更豐富的特征表達。此外，為了減少網(wǎng)絡(luò)參數(shù)、提升訓(xùn)練速度，首先使用1×1 卷積核降維特征，后續(xù)進行大尺度卷積操作。

Fig.1 Incepiton module圖1 Inception模塊

2.2.2 EdgeConv層

GCN 等圖神經(jīng)網(wǎng)絡(luò)在輸入前就已構(gòu)建了鄰接矩陣，節(jié)點與節(jié)點間的鄰接關(guān)系固定不變。以蛋白質(zhì)為例，一個蛋白質(zhì)圖G一般通過以下方法構(gòu)建：

式中：e為圖的邊；v為圖的節(jié)點；i、j為殘基索引；C代表殘基中的Cβ原子坐標；dmax為構(gòu)建邊的閾值。

當(dāng)前，盡管GCN 在質(zhì)量評估領(lǐng)域取得了優(yōu)異成果，但隨著網(wǎng)絡(luò)層次加深，殘基間的鄰接關(guān)系逐漸從初始的歐幾里得距離變成高維的特征間距離，初始的鄰接矩陣不再適合表達這種深層聯(lián)系。因此，根據(jù)特征構(gòu)建動態(tài)的鄰接關(guān)系尤為重要。

EdgeConv 原本是一種面向點云學(xué)習(xí)的網(wǎng)絡(luò)模塊，本文將殘基類比為特征空間上的點云，應(yīng)用該模塊挖掘殘基的特征信息。設(shè)R={r1，r2…rn}?RF為蛋白質(zhì)的殘基云，任意一個殘基表示為ri，鄰域節(jié)點表示為rj，F(xiàn)為殘基特征的通道數(shù)。如圖2（a）所示，殘基ri在EdgeConv 層中計算與其他所有殘基的特征距離|ri-rj|，選取距離最小的k個殘基作為其在本層的鄰域（圖中k=4），整個過程動態(tài)變化，即每一個EdgeConv 層都會重新計算ri的鄰域來適應(yīng)特征的變化。圖2（b）展示了動態(tài)圖網(wǎng)絡(luò)的節(jié)點更新過程，計算公式如式（5）所示。

Fig.2 EdgeConv operation圖2 EdgeConv操作

式中：ε表示ri與鄰域節(jié)點rj構(gòu)成的邊集合。

hθ的定義為：

式中：為MLP，操作對象為ri的全局特征（ri）和局部邊特征（ri-rj）。

2.2.3 總體架構(gòu)

本文基于Inception、EdgeConv 模塊，搭建了網(wǎng)絡(luò)整體架構(gòu)DGCQA，如圖3 所示。網(wǎng)絡(luò)初始輸入的節(jié)點并不是殘基，而是原子，原子特征選用其三維坐標(x，y，z)。DGCQA 首先在原子尺度上進行兩次EdgeConv 操作（鄰域范圍k=40），以充分挖掘殘基的原子幾何特征，再通過Cβ原子特征代替殘基，降采樣到殘基尺度的點云后進行后續(xù)3 層的EdgeConv 操作（鄰域范圍k=10）。

Fig.3 DGCQA architecture圖3 DGCQA結(jié)構(gòu)

網(wǎng)絡(luò)中，同尺度下不同EdgeConv 層輸出都會進行拼接，以充分利用網(wǎng)絡(luò)在不同深度下的特征信息。最后，模型將特征壓縮至適應(yīng)大小，并與多尺度卷積模塊的輸出拼接，通過全連接層和sigmoid 函數(shù)（圖3 中的σ）得到預(yù)測結(jié)果。EdgeConv 層的網(wǎng)絡(luò)設(shè)計如圖4 所示，首先獲取殘基節(jié)點k鄰域范圍內(nèi)的節(jié)點信息，然后使用MLP 對所有節(jié)點及其鄰域特征進行處理，最后通過池化層聚合特征，Edge-Conv層前后節(jié)點數(shù)量保持一致。

Fig.4 EdgeConv layer圖4 EdgeConv層

DGCQA 結(jié)構(gòu)中的多尺度卷積模塊如圖5 所示，該模塊用于捕獲不同視角下的殘基對特征關(guān)系。以殘基對特征作為輸入，經(jīng)過3 層Inception 結(jié)構(gòu)的處理后通過pooling 層進行降采樣，然后將特征展平為輸出通道數(shù)Mout大小。

Fig.5 Multi-scale convolution module圖5 多尺度卷積模塊

3 實驗結(jié)果及分析

3.1 實驗環(huán)境

本文實驗基于Ubuntu 18.04 操作系統(tǒng)，Pytorch 1.13.0深度學(xué)習(xí)框架，處理器為AMD Ryzen 95900X，顯卡為NVIDIA RTX3090。

3.2 數(shù)據(jù)集

實驗數(shù)據(jù)集來自CASP 9～CASP 13 比賽中官方提供的目標蛋白質(zhì)和各參賽小組提供的預(yù)測蛋白質(zhì)。針對每一個目標蛋白質(zhì)（target），均有多個小組提交的候選模型（decoy）與之對應(yīng)。本文將CASP 9～CASP 12 整合成85 000 個蛋白質(zhì)的數(shù)據(jù)集，以目標蛋白質(zhì)為基準隨機劃分訓(xùn)練集（270 個目標蛋白質(zhì)）和驗證集（50 個目標蛋白質(zhì)）。同時，將CASP 13 階段兩個比賽中的14 000 個蛋白質(zhì)模型（對應(yīng)72 個目標蛋白質(zhì)）作為測試集，將本文所提方法與其他方法進行比較。

3.3 評價指標

本文選用皮爾森相關(guān)系數(shù)（Pearson Correlation Coefficient）、斯皮爾曼相關(guān)系數(shù)（Spearman Correlation Coefficient）、均方根誤差（Root Mean Square Error）、均方根誤差（Root Mean Square Error）、標準分數(shù)（standard score）分析DGCQA 的性能，如式（7）—式（10）所示。

皮爾森相關(guān)系數(shù)用于評價兩個向量之間的線性關(guān)系，取值范圍為［-1，1］，接近0 表示弱相關(guān)，接近-1 或1 表示強相關(guān)。

式中：Y代表真實值向量；代表預(yù)測值向量；M代表對應(yīng)向量的平均值。

斯皮爾曼相關(guān)系數(shù)用于統(tǒng)計兩個連續(xù)變量之間的單調(diào)關(guān)系，對異常值敏感度較低，其值與變量具體值無關(guān)，僅與變量間的大小關(guān)系有關(guān)。

式中：R為預(yù)測值的取值等級；S為真實值的取值等級；MR、MS分別為R、S的均值。

均方根誤差是回歸模型的典型評價指標，用于衡量預(yù)測值與真實值之間的偏差。

式中：n為預(yù)測的實例數(shù)量；yi為真實值為預(yù)測值。

標準分數(shù)也叫z 分數(shù)（z-score），本文在評估全局質(zhì)量分數(shù)GDT_TS［21］（Global Distance Test_Total Score）時，根據(jù)DGCQA 的預(yù)測結(jié)果，從每個目標蛋白質(zhì)的候選模型中選擇一個預(yù)測評分最高的蛋白質(zhì)模型，并計算z-score 的平均值，值越高代表DGCQA 對優(yōu)秀蛋白質(zhì)模型的選擇越準確。

式中：X為原始數(shù)據(jù)；為數(shù)據(jù)集的均值；σ為數(shù)據(jù)集的標準差。

3.4 結(jié)果分析

3.4.1 整體性能

由于蛋白質(zhì)的全局結(jié)構(gòu)精度分數(shù)GDT_TS 是CASP 評價蛋白質(zhì)模型質(zhì)量的重要指標，因此本文在全部蛋白質(zhì)目標上計算GDT_TS 的均方根誤差（RMSE）、皮爾森相關(guān)系數(shù)（R）、斯皮爾曼相關(guān)系數(shù)（ρ）和標準分數(shù)（z），并與12 種方法進行比較，具體數(shù)據(jù)如表1所示。

Table 1 GDT_ TS prediction evaluation indicators表1 GDT_TS預(yù)測評價指標

表1 中Rtarget為先在每個目標蛋白質(zhì)的候選模型中計算皮爾森相關(guān)系數(shù)，然后再求均值，值越高表明該方法能更好地根據(jù)具體目標蛋白質(zhì)候選模型的整體質(zhì)量進行排名。由此可見，本文所提方法雖然在標準分數(shù)z 上相較于GraphQA 較差，但也具有一定的競爭力，相較于其他方法在多項指標達到最優(yōu)。圖6 顯示了GDT_TS 的真實值與預(yù)測值的關(guān)系及其分布情況。

Fig.6 GDT_ TS marginal圖6 GDT_TS邊際

3.4.2 鄰域范圍影響

本文中動態(tài)圖卷積的鄰域范圍是通過模型中Edge-Conv 層的超參數(shù)k所設(shè)置，k的取值會在一定程度上影響模型性能。實驗中保持DGCQA 整體架構(gòu)不變，僅修改在殘基尺度上提取特征時的k鄰域大小，觀察模型在皮爾森系數(shù)下的得分情況。如圖7 所示，當(dāng)鄰域范圍k=10 時，模型的皮爾森相關(guān)系數(shù)達到最優(yōu)效果0.876，而當(dāng)鄰域選取變小或變大時性能呈下降趨勢。原因為當(dāng)鄰域選取過小時會導(dǎo)致動態(tài)圖卷積視野受限，無法充分捕獲節(jié)點間的特征關(guān)系；當(dāng)鄰域選取過大會使深層網(wǎng)絡(luò)下的節(jié)點特征過于相似，將影響模型對不同節(jié)點的區(qū)分能力。

Fig.7 Pearson correlation coefficient under different neighborhood ranges圖7 不同鄰域范圍下的皮爾森相關(guān)系數(shù)

3.4.3 遷移學(xué)習(xí)有效性分析

為了驗證遷移學(xué)習(xí)的有效性，保持模型框架和其他特征輸入不變，僅去除ESM-1b 編碼以評估該部分對模型的影響。如表2、圖8 所示，加入ESM-1b 特征能使DGCQA 在所有指標上存在不同程度的提升，尤其在Rtarget上ESM-1b特征幫助DGCQA 提升了0.032，超過了第二名GraphQA 的0.779，達到了目前最好效果。

Table 2 Comparison of predictions for the presence and absence of ESM-1b on the test set表2 測試集上有無ESM-1b的預(yù)測比較

Fig.8 ESM-1b feature impact圖8 ESM-1b特征影響

3.4.4 動態(tài)圖卷積效果分析

為了通過實驗公平比較傳統(tǒng)圖卷積和動態(tài)圖卷積在蛋白質(zhì)模型質(zhì)量評估上的性能。在動態(tài)圖卷積方法上選擇了僅使用EdgeConv 層搭建的簡單網(wǎng)絡(luò)DGCQA-RAW；在特征上使用蛋白質(zhì)的氨基酸one-hot 編碼作為GCN 和動態(tài)圖卷積神經(jīng)網(wǎng)絡(luò)的唯一輸入。

二者在5 個重要指標的得分如表3 所示。由此可知，DGCQA-RAW 的5 個指標中4 個超過GCN，證明在質(zhì)量評估任務(wù)中，圖節(jié)點的動態(tài)鄰域相較于固定鄰接關(guān)系對特征發(fā)掘更有效。

Table 3 Comparison of GCN and DGCNN performance in quality assessment表3 GCN和DGCNN在質(zhì)量評估中的性能比較

4 結(jié)語

本文提出基于動態(tài)圖卷積的質(zhì)量評估方法DGCQA 來解決傳統(tǒng)圖卷積的局限性，并結(jié)合遷移學(xué)習(xí)思想，引入蛋白質(zhì)預(yù)訓(xùn)練模型ESM-1b 編碼特征訓(xùn)練模型。實驗表明，DGCQA 的表現(xiàn)相較于傳統(tǒng)方法更優(yōu)秀，證明了所提方法的有效性。

然而，本文方法仍存在一定的缺陷。例如，動態(tài)圖雖然無需構(gòu)建固定鄰接關(guān)系，但構(gòu)建節(jié)點鄰域范圍的選擇依然是一個不確定因素，選取的大小會在一定程度上影響模型性能。下一步，將構(gòu)建一個對鄰域范圍具有自適應(yīng)能力的動態(tài)圖神經(jīng)網(wǎng)絡(luò)，以提升模型的魯棒性。