秦琪琪,丁學(xué)明,王金雷

(上海理工大學(xué) 光電信息與計算機工程學(xué)院,上海 200093)

1 引 言

蛋白質(zhì)是生命的重要組成部分,具有不同功能的蛋白質(zhì)在不同的生命活動中發(fā)揮作用,一旦缺失了關(guān)鍵蛋白質(zhì),生命過程將無法繼續(xù)進行,并且蛋白質(zhì)功能的準(zhǔn)確預(yù)測可以幫助人類針對性的開發(fā)有效藥物,因而了解蛋白質(zhì)的功能對推進醫(yī)學(xué)發(fā)展至關(guān)重要.蛋白質(zhì)序列和結(jié)構(gòu)隱含了其功能信息,高通量測序技術(shù)的快速發(fā)展加快了蛋白質(zhì)序列數(shù)據(jù)的增長速度,但是蛋白質(zhì)序列獲得功能注釋的速度遠(yuǎn)低于序列數(shù)據(jù)增長的速度,造成了已知蛋白質(zhì)序列數(shù)量與已具有功能注釋的蛋白質(zhì)數(shù)量之間的差距不斷增大,并且新序列缺乏實驗驗證的結(jié)構(gòu)等生物信息.UniProt[1]數(shù)據(jù)庫中已有2億多條序列,而其中只有約為1%的序列是經(jīng)過實驗驗證的[2],再次表明實驗驗證的速度遠(yuǎn)低于序列數(shù)據(jù)增長速度,無法通過實驗驗證快速注釋蛋白質(zhì)功能,因此通過計算方法來預(yù)測蛋白質(zhì)功能開始成為探索蛋白質(zhì)的重要步驟.傳統(tǒng)的蛋白質(zhì)功能預(yù)測方法不僅耗時耗力而且成本較高,無法達(dá)到較高的預(yù)測準(zhǔn)確率,蛋白質(zhì)序列數(shù)據(jù)的快速增長,深度學(xué)習(xí)的不斷發(fā)展,都為準(zhǔn)確預(yù)測蛋白質(zhì)功能創(chuàng)造了條件,大量數(shù)據(jù)的學(xué)習(xí)和有效模型的構(gòu)建可以大大提高預(yù)測準(zhǔn)確率.

早期利用序列相似的方法預(yù)測蛋白質(zhì)功能,比如BLAST[3]通過序列比對進行功能注釋,然而高度不相似的蛋白質(zhì)序列也可能具有相同的功能,所以這種方法預(yù)測準(zhǔn)確率較低.針對該問題人們開始對BLAST方法進行改進,GOtcha[4]利用BLAST搜索結(jié)果進行功能注釋,PFP[5]和ESG[6]使用E-value加權(quán)預(yù)測和機器學(xué)習(xí)的方法優(yōu)化預(yù)測性能,都在BLAST的基礎(chǔ)上提高了預(yù)測準(zhǔn)確率,但是由于這些方法依然是依賴于序列相似性的,所以可靠性較低.因此人們開始結(jié)合蛋白質(zhì)多種生物信息預(yù)測其功能,CombFunc[7]將ConFunc與蛋白質(zhì)序列、結(jié)構(gòu)域或蛋白質(zhì)-蛋白質(zhì)相互作用等信息結(jié)合進行功能預(yù)測,使用SVM組合這些信息,多種信息的結(jié)合使其預(yù)測性能優(yōu)于BLAST,但仍然不能實現(xiàn)準(zhǔn)確預(yù)測功能.隨著數(shù)據(jù)的增長,機器學(xué)習(xí)算法和神經(jīng)網(wǎng)絡(luò)等方法被大量用于蛋白質(zhì)功能預(yù)測.根據(jù)相互作用蛋白質(zhì)共享功能的原理,蛋白質(zhì)相互作用網(wǎng)絡(luò)(PPI)[8-10]得到了人們的關(guān)注.Hishigaki[11]等人根據(jù)蛋白質(zhì)相互作用圖,引用n鄰近蛋白質(zhì)的方法,通過統(tǒng)計最大卡方值來預(yù)測蛋白質(zhì)功能.2020年Cai Y等人提出了集成深度學(xué)習(xí)模型SDN2GO[12],該模型利用蛋白質(zhì)序列、結(jié)構(gòu)域和PPI 網(wǎng)絡(luò)構(gòu)建神經(jīng)網(wǎng)絡(luò)子模型,并通過加權(quán)分類器整合子模型進行功能預(yù)測.這些方法都使用了蛋白質(zhì)結(jié)構(gòu)和PPI網(wǎng)絡(luò)等蛋白質(zhì)序列之外的生物信息,而新測序的序列往往缺失這些豐富的生物信息,所以這些方法會受到相關(guān)信息缺失的限制.

由于蛋白質(zhì)可具有多個功能,所以蛋白質(zhì)功能預(yù)測可視為是一個多標(biāo)簽分類問題,適用于多標(biāo)簽分類的代表技術(shù)有支持向量機、決策樹、K鄰近方法和人工神經(jīng)網(wǎng)絡(luò)等[13],DeepGO[14]采用卷積神經(jīng)網(wǎng)絡(luò)(CNN)和PPI率先實現(xiàn)了較高的功能預(yù)測準(zhǔn)確率,表明CNN可以從蛋白質(zhì)序列和結(jié)構(gòu)中提取有用信息,而圖卷積神經(jīng)網(wǎng)絡(luò)(GCN)可以較好的解決CNN只能處理歐式數(shù)據(jù)的問題,結(jié)合新測序序列缺失豐富的生物信息問題,本文提出了基于蛋白質(zhì)序列和組合圖卷積網(wǎng)絡(luò)的蛋白質(zhì)功能預(yù)測模型(Protein Function Prediction using Sequences and Combined Graph Convolutional Networks,PFP-SCGCN),對GO ontology(GO)和Enzyme Commission(EC)兩種蛋白質(zhì)功能大類進行預(yù)測,并通過Grad-CAM(Gradient-weighted Class Activation Mapping)[15]方法分析PFP-SCGCN中的組合圖卷積網(wǎng)絡(luò)層識別功能位點.

2 方法概述

本文所提的蛋白質(zhì)功能預(yù)測模型PFP-SCGCN主要由序列特征提取模塊、進化信息提取模塊、組合圖卷積和全連接層塊(FC blocks)組成,具體結(jié)構(gòu)見圖1,為了識別功能位點,利用Grad-CAM方法對組合圖卷積的輸出進行分析計算.

圖1 PFP-SCGCN模型結(jié)構(gòu)Fig.1 PFP-SCGCN model structure

2.1 序列特征提取

對于蛋白質(zhì)序列,可以通過神經(jīng)網(wǎng)絡(luò)將序列氨基酸映射成空間向量,以此來獲得豐富的信息.在PFP-SCGCN中先對蛋白質(zhì)序列進行one-hot編碼,在one-hot時考慮了20種常見氨基酸和空位gap,所以經(jīng)過one-hot編碼后的序列為S(L×21).由于考慮到當(dāng)?shù)鞍踪|(zhì)序列較長時,其one-hot后的矩陣將會非常稀疏,不利于網(wǎng)絡(luò)提取信息,所以將one-hot后的序列通過全連接層(Dense layer)進行信息稠密化,以此獲得序列的空間嵌入信息.為獲取更多的蛋白質(zhì)序列信息,通過預(yù)訓(xùn)練好的蛋白質(zhì)語言模型ESM2[16]將蛋白質(zhì)序列中的每個氨基酸映射成空間表達(dá)式,L個氨基酸組成的序列經(jīng)過ESM2處理后變?yōu)長個空間向量組成的特征矩陣(Feature matrix),維度為L×1280,該矩陣表達(dá)了序列氨基酸之間關(guān)聯(lián)的隱層特征,將該隱層特征利用全連接層(Dense layer)進行再次映射,與獲得的序列空間嵌入信息一起作為序列特征信息,即組合圖卷積的初始輸入.

2.2 進化信息提取

長期以來人們一直認(rèn)為蛋白質(zhì)的同源序列可以為了解其功能和結(jié)構(gòu)提供重要的信息,比如在進化過程中保守的氨基酸通常與功能密切相關(guān),所以同源序列隱含的進化信息可以為了解和設(shè)計蛋白質(zhì)功能提供幫助,而且通過同源序列提取的進化信息有助于蛋白質(zhì)的結(jié)構(gòu)預(yù)測[17].除了通過進化信息探索蛋白質(zhì)功能和結(jié)構(gòu),Ji-Yong An等人[18]利用進化信息和循環(huán)神經(jīng)網(wǎng)絡(luò)(RNN)預(yù)測自相互作用蛋白質(zhì),Shanwen Sun[19]等人利用進化信息,通過建立支持向量機分類器來識別抗凍蛋白,因此進化信息蘊含的豐富信息可以幫助人們更好地了解蛋白質(zhì).

目前已有較多的方法來提取蛋白質(zhì)的進化信息,統(tǒng)計耦合分析(SCA)[20]通過對MSA的分析來識別共同進化的氨基酸對,為設(shè)計蛋白質(zhì)折疊和識別功能區(qū)域提供了有用的信息.互信息(Mutual Information,MI)[21]可以識別蛋白質(zhì)三維結(jié)構(gòu)中一些直接相互作用的殘基對,但是MI在計算耦合分?jǐn)?shù)時只考慮了兩個對齊列的位置.直接耦合分析(DCA)也可以在大量序列中尋找成對氨基酸之間的耦合關(guān)系[22].本文通過Leri[23]來提取進化信息,先為每條蛋白質(zhì)序列搜索對應(yīng)的多序列比對(MSA),Leri將MSA作為輸入,考慮了MSA中所有位置氨基酸間的相互作用,通過序列產(chǎn)生器建模獲得序列中任意氨基酸對之間的耦合項eij,序列中所有氨基酸對耦合項eij的組合形成進化耦合信息矩陣(Evolutionary coupling information,EVCI),對EVCI進行特征值計算,推斷出與蛋白質(zhì)功能相關(guān)的氨基酸殘基群落(Residue Communities,RC).

2.2.1 多序列比對(MSA)

為了獲取進化信息,利用多序列比對(MSA)來尋找同源序列,對于數(shù)據(jù)集中的序列,通過HHblits[24]工具在默認(rèn)參數(shù)E-value為0.001的情況下,從Uniclust30[25]數(shù)據(jù)庫中5次迭代搜索得到每一條蛋白質(zhì)序列對應(yīng)的MSA.

2.2.2 進化耦合信息EVCI提取

進化耦合信息(EVCI)包含了蛋白質(zhì)序列中的所有氨基酸,EVCI中任意氨基酸對之間的耦合項eij,表達(dá)了進化過程中任意氨基酸對之間的相關(guān)程度,可以幫助人們判斷每個氨基酸對于蛋白質(zhì)功能或結(jié)構(gòu)的貢獻大小.獲得每一條蛋白質(zhì)序列的MSA后,利用Leri[23]提取對應(yīng)的EVCI,具體內(nèi)容如下:

1)通過HHblits工具得到每一條序列的MSA后,需要對MSA進行過濾處理,以去除MSA中低質(zhì)量位點和序列.處理標(biāo)準(zhǔn)為兩個,分別是當(dāng)在MSA中的單個位點超過90%是空位gap時,移除該位點;當(dāng)在MSA中的某條序列超過80%是空位gap時,移除該條序列.在MSA中每條序列的權(quán)重ω(τ)是由序列一致性(sequence identity,I)計算得到的,I是通過序列τ與其它所有序列之間的漢明距離(Hamming distance)DH(τ,τj)衡量的,具體權(quán)重計算如式(1)所示:

ω(τ)=(∑jI[DH(τ,τj)<θ])-1

(1)

其中參數(shù)θ為默認(rèn)值0.2.

2)將處理后的MSA作為輸入,利用基于偽似然最大化方法的全局概率模型[26]來捕獲序列中氨基酸對之間的進化耦合關(guān)系.對蛋白質(zhì)家族進化過程進行序列產(chǎn)生器建模,根據(jù)其家族中所有序列的空間分布,以概率P(τ)產(chǎn)生序列τ,P(τ)的計算為式(2):

(2)

其中Z為配分函數(shù),通過將家族中所有序列的玻爾茲曼因數(shù)相加來進行歸一化.每條序列的玻爾茲曼因數(shù)的測量是由進化統(tǒng)計能量E(τ)(馬爾可夫隨機場或Potts模型)定義的,具體如式(3)所示:

E(τ)=∑i

(3)

其中eij為MSA中處于位置i的氨基酸與處于位置j的氨基酸之間的耦合項,hi為特定位點偏差項.eij和hi表達(dá)了兩種類型的約束,分別是每對氨基酸之間的成對約束,以及特定位點偏差的約束,考慮了所有可能位置之間的相互作用,這是互信息沒有考慮到的.對于eij和hi參數(shù),利用位置分解的最大偽似然方法進行優(yōu)化,如式(4)所示:

(4)

(5)

2.2.3 氨基酸殘基群落RC提取

氨基酸殘基群落(RC)包含的是蛋白質(zhì)序列中的部分氨基酸,這部分氨基酸通常與蛋白質(zhì)功能相關(guān).RC用于表達(dá)進化過程中與功能相關(guān)的氨基酸殘基對之間的耦合程度,因此可以利用RC信息著重關(guān)注序列中的部分氨基酸,去除部分噪聲并減輕網(wǎng)絡(luò)負(fù)擔(dān),使網(wǎng)絡(luò)抓取與蛋白質(zhì)功能相關(guān)的有用信息.Leri[23]利用譜分解法[27]從進化耦合信息(EVCI)中提取得到RC,譜分解法可以識別EVCI中的強耦合殘基對,將氨基酸殘基分類為不同的氨基酸殘基群落(RC),具體內(nèi)容如下:

對EVCI進行特征值計算和排序,通過雅可比迭代過程確定EVCI的正特征值和對應(yīng)的特征向量,取其前5個特征值中的兩個及其對應(yīng)的特征向量來定義與功能相關(guān)的氨基酸殘基群落(RC),將5個特征向量定義為Vk|k=1,2,…,5,本文使用的RC由兩部分組成,分別定義為:

其中ε為默認(rèn)參數(shù)0.05.圖1中進化信息提取模塊展示了使用Leri提取EVCI和RC的主要流程,其中MSA為N×L,代表搜索到的N條序列;EVCI為L×L,代表長度為L的序列中每對氨基酸之間的耦合關(guān)系;RC為m×m,是從EVCI中根據(jù)特征值計算得到的,代表了與功能相關(guān)的氨基酸殘基之間的耦合關(guān)系,是蛋白質(zhì)序列中的部分氨基酸且具有強耦合關(guān)系,所以m

2.3 組合圖卷積和全連接層塊(FC blocks)

PFP-SCGCN模型采用組合圖卷積網(wǎng)絡(luò)分析蛋白質(zhì)序列,組合圖卷積網(wǎng)絡(luò)定義為:采用3個通道的圖卷積網(wǎng)絡(luò),先將其中兩通道圖卷積網(wǎng)絡(luò)并行,并行的兩通道圖卷積分別由兩層網(wǎng)絡(luò)構(gòu)成,每層網(wǎng)絡(luò)輸出的節(jié)點特征維度為512,進行并行的圖卷積采用相同的初始輸入,并且分別采用兩種鄰接矩陣來進行信息融合,實現(xiàn)信息特征的多樣提取;接著使用另一通道圖卷積網(wǎng)絡(luò)與并行的兩通道圖卷積網(wǎng)絡(luò)進行串聯(lián),用于處理并行圖卷積的輸出,進行串聯(lián)的單通道圖卷積由一層網(wǎng)絡(luò)構(gòu)成,該層網(wǎng)絡(luò)輸出的節(jié)點特征維度為512,初始輸入為并行兩通道圖卷積的輸出,鄰接矩陣采用并行圖卷積中使用的兩種鄰接矩陣中的一個,以此對并行圖卷積輸出的信息特征進行噪聲過濾和蛋白質(zhì)功能相關(guān)信息再次提取,加大了蛋白質(zhì)序列中可能與功能相關(guān)的氨基酸節(jié)點的關(guān)注度,使進行串聯(lián)的單通道圖卷積輸出包含豐富的功能相關(guān)的信息特征.3個通道圖卷積的使用實現(xiàn)了圖卷積網(wǎng)絡(luò)的串并聯(lián)組合,從而創(chuàng)建組合圖卷積網(wǎng)絡(luò).

圖卷積網(wǎng)絡(luò)(GCN)在卷積神經(jīng)網(wǎng)絡(luò)(CNN)的基礎(chǔ)上考慮了圖結(jié)構(gòu)信息[28],隨著圖數(shù)據(jù)的不斷豐富,GCN已經(jīng)被廣泛應(yīng)用于各個領(lǐng)域.圖G(V,E,A)中,V代表了圖中的節(jié)點集合,E表達(dá)了節(jié)點之間的連接邊,A稱為鄰接矩陣,表達(dá)了節(jié)點之間是否存在連接.GCN通過層與層之間的傳遞來更新節(jié)點特征,如式(6)所示:

fl+1=σ(fl,A)

(6)

其中fl為GCN中第l層網(wǎng)絡(luò)的輸出,當(dāng)l為0時,f0為網(wǎng)絡(luò)的初始特征輸入,每一層中的具體處理方式如式(7)所示:

(7)

本文將蛋白質(zhì)序列氨基酸視為圖中節(jié)點,將提取的EVCI和RC在同等級別上進行歸一化,值大于等于0.2閾值的對應(yīng)位置視為氨基酸對之間存在聯(lián)系,以此生成兩種鄰接矩陣,即鄰接矩陣1和鄰接矩陣2.EVCI衡量了蛋白質(zhì)序列中所有氨基酸對之間的耦合關(guān)系,其包含了非常豐富的信息,RC衡量了蛋白質(zhì)序列中部分氨基酸殘基對之間的強耦合關(guān)系,其通常與蛋白質(zhì)功能相關(guān),所以RC去除了EVCI中的部分噪聲.組合圖卷積中并行的兩通道GCN將序列特征提取模塊提取的蛋白質(zhì)序列特征信息作為初始輸入f0,其中一個通道GCN采用了根據(jù)EVCI生成的鄰接矩陣1,考慮序列中所有氨基酸之間的相互作用,而另一通道GCN采用根據(jù)RC生成的鄰接矩陣2,只認(rèn)為與功能相關(guān)的氨基酸之間存在聯(lián)系,忽略其它氨基酸之間的作用.兩通道GCN的并行將EVCI和RC的信息結(jié)合起來,實現(xiàn)了既通過RC去除EVCI中包含的部分噪聲,又考慮EVCI中包含的可能有用的豐富信息,其中并行兩通道GCN分別由兩層512通道數(shù)的網(wǎng)絡(luò)構(gòu)成.并行GCN之后又采用一層通道數(shù)為512的GCN進行串聯(lián),串聯(lián)GCN將并行GCN的輸出作為輸入,使用根據(jù)RC生成的鄰接矩陣2,實現(xiàn)了噪聲的再次過濾.組合圖卷積的輸出經(jīng)過一個由3層全連接網(wǎng)絡(luò)構(gòu)成的分類器(FC blocks)進行分析,提取出蛋白質(zhì)功能類別的分類特征,每層的通道數(shù)為1024,最后經(jīng)過Softmax激活函數(shù)輸出每個功能類別的預(yù)測概率.

設(shè)計電壓比較電路，供5V直流電，用于比較的電路電阻R34.5kΩ、R41kΩ，可變電阻 R5初值 0.5kΩ，25℃時熱敏電阻R210kΩ，UA點電壓2.5V，UB點電壓約3.75V，電壓比較器LM311正極電壓低于負(fù)極，不通，無輸出。溫度升高到60℃時，熱敏電阻R2降到3kΩ，電壓約1.2V，UA點電壓升到3.8V，高于UB點電壓，電壓比較器正極電壓高于負(fù)極，接通，輸出信號?？勺冸娮栌糜谠O(shè)置比較電壓的閾值，可變電阻阻值越大，UB點電壓越低，報警溫度越低。

2.4 功能位點識別

氨基酸殘基群落(RC)包含了與功能相關(guān)的氨基酸殘基信息,為了探討RC信息的實用性,利用Grad-CAM[15]方法計算特定氨基酸殘基為功能位點的概率.Grad-CAM是一個類判別定位技術(shù),為基于CNN模型的預(yù)測提供了視覺解釋,可以利用該方法計算出網(wǎng)絡(luò)對每個氨基酸的關(guān)注度,從而得到特定氨基酸殘基為功能位點的概率.采用Grad-CAM根據(jù)組合圖卷積的輸出進行梯度值計算,得到序列中每個氨基酸為功能位點的置信分?jǐn)?shù),分?jǐn)?shù)越高則表示其為功能位點的概率越高,并將預(yù)測功能位點結(jié)果與BioLip[29]進行對比,BioLip數(shù)據(jù)庫包含了蛋白質(zhì)真實功能位點信息.利用組合圖卷積網(wǎng)絡(luò)最后一層GCN輸出的特征,首先計算輸出特征圖H中每個通道對于特定功能類別的權(quán)重值,如式(8)所示:

(8)

(9)

其中Hk,c為網(wǎng)絡(luò)針對功能類別c的輸出特征圖H的k通道值,通過ReLU激活函數(shù)保留與功能類別c相關(guān)的特征.

3 結(jié)果分析

3.1 構(gòu)建數(shù)據(jù)集

對于眾多的蛋白質(zhì)功能,可廣泛分為GO ontology(GO)和Enzyme Commission(EC),其中GO 可劃分為3大類:分子功能(Molecular Function,MF)、生物過程(Biological Process,BP)、細(xì)胞成分(Cellular Component,CC),EC功能注釋通過EC編號表達(dá),EC編號由4個圓點隔開的數(shù)字組成.本文對GO和EC兩種功能大類進行預(yù)測,從SIFTS[30]數(shù)據(jù)庫中下載蛋白質(zhì)功能注釋信息和蛋白質(zhì)序列,該數(shù)據(jù)庫提供了來自蛋白質(zhì)數(shù)據(jù)庫PDB[31]和UniProt知識庫(UniProtKB)[1]的信息,結(jié)合了蛋白質(zhì)結(jié)構(gòu)和序列的信息,有利于在蛋白質(zhì)結(jié)構(gòu)和序列之間傳遞蛋白質(zhì)功能注釋.在從SIFTS下載的數(shù)據(jù)中,對功能類別和序列進行篩選.對于GO功能注釋保留根據(jù)EXP(Inferred from experiment)、IDA(Inferred from direct assay)、IPI(Inferred from physical interaction)、IMP(Inferred from mutant phenotype)、IGI(Inferred from genetic interaction)、IEP(Inferred from expression pattern)、TAS(Traceable author statement)和IC(Inferred by curator)獲得的注釋,并且為了保證每個功能類別能有足夠的蛋白質(zhì)數(shù)據(jù),當(dāng)具有某個功能類別的蛋白質(zhì)數(shù)量介于50～5000之間時,就把該功能類別添加到總的功能類別中,EC功能注釋只保留了三級和四級的EC編碼.對于所下載數(shù)據(jù)集中的蛋白質(zhì)序列,只保留含有20種常見氨基酸的序列,利用MMseqs2聚類工具以95%的序列一致性去冗余,并且序列長度小于等于1024.

處理后的數(shù)據(jù)集中,GO總共有34501個蛋白質(zhì),EC總共有18336個蛋白質(zhì),GO中MF,BP和CC總的功能類別數(shù)分別為5013個,12248個,1688個,EC的總功能類別數(shù)為577個.一個蛋白質(zhì)可以具有多個功能,具有GO功能的蛋白質(zhì)可以同時具有屬于MF,BP,CC這3大類別的功能,有的則可能只具有屬于其中兩大類或者一類的功能,對數(shù)據(jù)集隨機打亂后按照大約8∶1∶1的比例劃分?jǐn)?shù)據(jù)集,具體的數(shù)據(jù)分布如表1所示.

表1 GO和EC的數(shù)據(jù)集分布Table 1 GO and EC data set distribution

3.2 評價指標(biāo)

精準(zhǔn)率(Precision),召回率(Recall)和F1-score是常用的分類評價指標(biāo),Precision表達(dá)了對蛋白質(zhì)預(yù)測出的功能類別中為真實功能類別的概率,Recall表達(dá)了蛋白質(zhì)的真實功能類別中被正確預(yù)測出的個數(shù)占比,F1-score則根據(jù)Precision和Recall的值計算得出,計算方式為式(10):

(10)

F1-score綜合評價了Precision和Recall.

(11)

其中TP為正確預(yù)測為正樣本的個數(shù),TN為正確預(yù)測為負(fù)樣本的個數(shù),FP為錯將負(fù)樣本預(yù)測為正樣本的個數(shù),FN為錯將正樣本預(yù)測為負(fù)樣本的個數(shù).

3.3 實驗環(huán)境及參數(shù)設(shè)置

本文基于TensorFlow2框架搭建PFP-SCGCN模型,通過NVIDIA GeForce RTX 2080Ti GPU進行模型訓(xùn)練,選擇Adam優(yōu)化器來優(yōu)化模型參數(shù),設(shè)置學(xué)習(xí)率為0.001,batch_size為64,訓(xùn)練周期為100 epoch,為了防止模型過擬合,采用early stopping機制,周期為10,損失函數(shù)為二進制交叉熵?fù)p失函數(shù),計算方式如式(12)所示:

(1-yij)log(1-P(yij))]

(12)

其中N為總共的樣本數(shù),C指總共的功能類別數(shù),yij為蛋白質(zhì)i對于功能類別j的真實標(biāo)簽,P(yij)則為其對應(yīng)的預(yù)測概率.

3.4 對比實驗

模型PFP-SCGCN主要與5個方法進行對比,包括BLAST[3]、DeepGOPlus[32]、FunFams[33]、DeepFRI[34]和GAT-GO[35],具體如下:

1)BLAST[3]:屬于無監(jiān)督方法,利用序列相似性進行功能注釋,本文采用了與CAFA1[36]中相同的BLAST使用方法.設(shè)E-value為默認(rèn)值1e-3,將訓(xùn)練集中序列的功能注釋給對應(yīng)的測試集序列,預(yù)測分?jǐn)?shù)為序列一致性分?jǐn)?shù).

2)DeepGOPlus[32]:使用蛋白質(zhì)序列信息進行功能預(yù)測,為基于一維卷積神經(jīng)網(wǎng)絡(luò)(1DCNN)的模型,將通過one-hot編碼后的序列經(jīng)過16層不同卷積核大小的1DCNN進行蛋白質(zhì)功能預(yù)測.

3)FunFams[33]:屬于無監(jiān)督方法,是一種基于域的蛋白質(zhì)功能注釋方法,利用CATH superfamilies的功能分類進行蛋白質(zhì)功能注釋,對每個蛋白質(zhì)序列通過HMMER3[37]工具在CATH FunFams中搜索,將HMM得分最高的功能類別注釋給測試序列,計算搜索結(jié)果中功能類別出現(xiàn)的頻率為預(yù)測分?jǐn)?shù),并且GO功能類別預(yù)測分?jǐn)?shù)按照Das[33]等人所述在GO層次結(jié)構(gòu)中向上傳播.

4)DeepFRI[34]:使用蛋白質(zhì)序列信息和結(jié)構(gòu)信息作為模型的輸入,通過預(yù)訓(xùn)練好的 LSTM(Long Short Term Memory)模型提取序列特征,根據(jù)蛋白質(zhì)結(jié)構(gòu)構(gòu)建氨基酸之間的接觸圖,將LSTM的輸出與構(gòu)建的接觸圖一起輸入給GCN來預(yù)測蛋白質(zhì)功能.

5)GAT-GO[35]:使用蛋白質(zhì)序列信息和預(yù)測的殘基接觸圖作為模型的輸入.首先生成每條蛋白質(zhì)序列對應(yīng)的PSSM(Position Specific Scoring Matrix),再根據(jù)PSSM通過RaptorX[38]預(yù)測對應(yīng)的殘基接觸圖,對蛋白質(zhì)序列進行one-hot編碼,并利用ESM-1b[16]獲得殘基級嵌入信息和蛋白質(zhì)嵌入信息,通過一維卷積神經(jīng)網(wǎng)絡(luò)(1DCNN)處理序列one-hot編碼、PSSM和殘基級嵌入信息,將1DCNN的輸出和預(yù)測的殘基接觸圖通過圖注意力網(wǎng)絡(luò)(GAT)進行信息融合,最后將GAT的輸出與蛋白質(zhì)嵌入信息輸入給全連接層進行蛋白質(zhì)功能預(yù)測.

3.5 實驗結(jié)果對比分析

根據(jù)Fmax,AUPR和MCC值的比較,對比PFP-SCGCN與BLAST、DeepGOPlus、FunFams、DeepFRI和GAT-GO這5種方法的預(yù)測性能差異,所對比方法均按照其默認(rèn)參數(shù)進行設(shè)置,不同方法在MF,BP,CC和EC中的預(yù)測性能對比如表2,結(jié)果表明模型PFP-SCGCN在所有評價指標(biāo)中幾乎都達(dá)到了最好的性能.對比結(jié)果中DeepFRI各項指標(biāo)值排名靠前,這與其使用了蛋白質(zhì)結(jié)構(gòu)信息有關(guān),PFP-SCGCN的各項指標(biāo)值遠(yuǎn)高于DeepFRI,表明了PFP-SCGCN只用序列信息而不用結(jié)構(gòu)等其它生物信息就可以較準(zhǔn)確的預(yù)測功能.這是因為DeepFRI中使用的殘基接觸圖是蛋白質(zhì)三級結(jié)構(gòu)到二維空間的映射,在殘基接觸圖中存在連接的氨基酸對只能代表其在空間中距離較近,但其不一定是功能結(jié)合位點,而PFP-SCGCN中使用了進化耦合信息(EVCI)和氨基酸殘基群落(RC),其中EVCI包含了非常豐富的信息,衡量了蛋白質(zhì)序列在進化過程中氨基酸對之間的關(guān)聯(lián)程度,進化過程中保守的氨基酸很有可能與蛋白質(zhì)功能相關(guān),EVCI中較大的值對應(yīng)位置的氨基酸對具有強耦合關(guān)系,說明該氨基酸對在進化過程中關(guān)聯(lián)密切,則其可能與蛋白質(zhì)功能有很大的相關(guān)性,RC包含了EVCI中具有強耦合關(guān)系的氨基酸,故RC擁有蛋白質(zhì)功能相關(guān)信息,過濾了EVCI中的部分噪聲.PFP-SCGCN在組合圖卷積網(wǎng)絡(luò)中同時使用了EVCI和RC,既考慮了EVCI中的多樣信息,又通過RC引導(dǎo)網(wǎng)絡(luò)關(guān)注與功能相關(guān)的氨基酸節(jié)點,以此使網(wǎng)絡(luò)極大程度的開發(fā)與蛋白質(zhì)功能相關(guān)的重要信息,從而使PFP-SCGCN更好的預(yù)測蛋白質(zhì)功能.雖然GAT-GO模型輸入中采用了預(yù)測的殘基接觸圖,但是在GO中GAT-GO的預(yù)測性能整體優(yōu)于DeepFRI,這與其使用了GAT網(wǎng)絡(luò)有關(guān),GAT中的自注意力機制給圖中每個節(jié)點分配不同程度的注意力,減少了對圖結(jié)構(gòu)的依賴程度,因而使GAT-GO在使用預(yù)測的殘基接觸圖時也可以較準(zhǔn)確的預(yù)測蛋白質(zhì)功能,但是預(yù)測的殘基接觸圖難免含有噪聲,并且預(yù)測出具有接觸關(guān)系的殘基對不一定是與功能相關(guān)的.FunFams在EC中的預(yù)測準(zhǔn)確率排名第二,這與其使用的CATH superfamilies數(shù)據(jù)有關(guān),使EC功能類別出現(xiàn)的頻率較高.BLAST方法的預(yù)測性能在這5個對比方法中較差,這也表明僅基于序列相似性預(yù)測功能是不夠準(zhǔn)確的.DeepGOPlus的性能高于BLAST,表明CNN可以從序列中提取有用的信息,但對于準(zhǔn)確預(yù)測功能還有距離.

表2 不同方法的功能預(yù)測性能比較Table 2 Performance comparison of function prediction between different methods

3.6 模型魯棒性檢驗

為了檢驗PFP-SCGCN的魯棒性,計算測試集中的序列與訓(xùn)練集中所有序列的序列一致性,按照序列一致性的最大值將GO和EC的測試集劃分為了30%,40%,50%,70%,95%這5組,選擇最大序列一致性≤30%這組測試集的預(yù)測性能進行比較,表3展示了具體結(jié)果.從表3對比結(jié)果可以看出,當(dāng)測試集序列與訓(xùn)練集序列極度不相似時,所提模型PFP-SCGCN依然可以達(dá)到較高的預(yù)測準(zhǔn)確率,并且整體的評價指標(biāo)值都優(yōu)于對比方法,證明了PFP-SCGCN具有較好的魯棒性.PFP-SCGCN在≤30%序列一致性測試集的預(yù)測性能優(yōu)于DeepFRI,再次證明了PFP-SCGCN不使用蛋白質(zhì)結(jié)構(gòu)這種序列信息之外的生物信息,依然可以較準(zhǔn)確的預(yù)測蛋白質(zhì)功能.BLAST的預(yù)測性能在魯棒性對比結(jié)果中依然是較差的,表明了僅基于序列相似性進行功能預(yù)測的不可靠性.DeepGOPlus在魯棒性檢驗中的性能較穩(wěn)定,表明多層CNN結(jié)構(gòu)是可以提取部分有價值信息的,但仍需改進.

表3 在≤30%序列一致性中的預(yù)測性能比較Table 3 Comparison of prediction performance in ≤30% sequence identity

3.7 功能位點識別

圖2 功能位點識別示例Fig.2 Examples of functional site identification

圖3 功能位點識別示例對應(yīng)ROC曲線Fig.3 ROC curves of the examples of functional site identification

圖4 功能位點識別評估Fig.4 Functional site identification evaluation

4 結(jié)束語

序列提取技術(shù)的快速發(fā)展使蛋白質(zhì)序列數(shù)據(jù)增長速度加快,導(dǎo)致未功能注釋的蛋白質(zhì)數(shù)量不斷增多,并且新測序序列缺失豐富的生物信息,針對該問題,本文提出了僅利用蛋白質(zhì)序列信息的蛋白質(zhì)功能預(yù)測模型PFP-SCGCN.該模型利用one-hot編碼,預(yù)訓(xùn)練好的蛋白質(zhì)語言模型ESM2和全連接網(wǎng)絡(luò)提取蛋白質(zhì)序列特征信息,通過MSA提取進化耦合信息(EVCI)和氨基酸殘基群落(RC),獲得序列中氨基酸對之間的耦合關(guān)系,構(gòu)建并行兩通道GCN和串聯(lián)單通道GCN結(jié)合的組合圖卷積網(wǎng)絡(luò),根據(jù)EVCI和RC生成表達(dá)氨基酸之間連接關(guān)系的兩種鄰接矩陣,將這兩種鄰接矩陣與序列特征信息一起作為組合圖卷積網(wǎng)絡(luò)的輸入,最后通過多層全連接網(wǎng)絡(luò)和Softmax激活函數(shù)預(yù)測蛋白質(zhì)功能類別,并通過Grad-CAM方法識別功能位點.PFP-SCGCN在MF,BP,CC和EC數(shù)據(jù)上的預(yù)測性能均優(yōu)于對比方法,并且具有較好的魯棒性,表明PFP-SCGCN可以在僅使用序列信息的情況下較準(zhǔn)確的預(yù)測蛋白質(zhì)功能.Grad-CAM計算值與BioLip的重合結(jié)果表明RC為PFP-SCGCN提供了有用的信息,使PFP-SCGCN具有一定的功能位點注釋能力.

雖然PFP-SCGCN已具有較好的功能預(yù)測能力和功能位點注釋能力,但EVCI和RC信息的使用仍有改進空間,并且PFP-SCGCN的功能位點識別能力可以進一步的提高,未來將嘗試不同的信息提取方式和新的網(wǎng)絡(luò)架構(gòu),以此進一步優(yōu)化模型性能.PFP-SCGCN相關(guān)代碼可在https://github.com/psxz1/PFP-SCGCN獲得.