王金雷 丁學(xué)明 秦琪琪 彭博雅

摘 要：蛋白質(zhì)的功能對(duì)于理解細(xì)胞和生物的活動(dòng)機(jī)制、研究疾病機(jī)理等至關(guān)重要。面對(duì)序列數(shù)據(jù)庫(kù)的快速增長(zhǎng)，傳統(tǒng)的實(shí)驗(yàn)和序列對(duì)比方法不足以支撐大規(guī)模的蛋白質(zhì)功能標(biāo)注。為此，提出EGNet（evolutionary graph network）模型，采用蛋白質(zhì)預(yù)訓(xùn)練語(yǔ)言模型ESM2和onehot編碼得到蛋白質(zhì)序列編碼，通過(guò)序列自注意力和物理計(jì)算整合出殘基間的協(xié)同進(jìn)化信息PI（paired interaction）和SPI（strong paired interaction）；之后將兩種進(jìn)化信息和序列編碼作為多層串聯(lián)圖卷積網(wǎng)絡(luò)輸入，學(xué)習(xí)序列編碼節(jié)點(diǎn)特征，實(shí)現(xiàn)端到端的蛋白質(zhì)功能預(yù)測(cè)。與早期方法相比，在ENZYME數(shù)據(jù)庫(kù)中的EC（Enzyme Commission）類別標(biāo)簽上，EGNet獲得了更好的性能，其Fscore達(dá)到0.89，AUPR值達(dá)到0.91。結(jié)果表明，EGNet僅僅采用單條序列來(lái)預(yù)測(cè)蛋白質(zhì)功能就可以得到良好的結(jié)果，從而能夠提供快速且有效的蛋白質(zhì)功能注釋。

關(guān)鍵詞：蛋白質(zhì)功能；深度學(xué)習(xí)；協(xié)同進(jìn)化信息；語(yǔ)言模型；圖卷積神經(jīng)網(wǎng)絡(luò)

中圖分類號(hào)：TP181?? 文獻(xiàn)標(biāo)志碼：A?? 文章編號(hào)：1001-3695（2023）12-008-3572-06

doi：10.19734/j.issn.10013695.2023.04.0166

Protein function prediction based on coevolutionary information and deep learning

Abstract：The function of protein is crucial for understanding the mechanisms of cellular and biological activities，as well as for studying the mechanisms of diseases.Traditional experimental and sequence alignment methods are insufficient to support largescale protein functional annotation when in the face of the rapid growth of sequence databases.For this situation，this paper proposed EGNet model，which utilized the protein pretraining language model ESM2 and onehot encoding to obtain the protein sequence encoding.The model integrated the coevolutionary information between residues，including PI and SPI，through sequence selfattention and physical calculations.Subsequently，the two types of coevolutionary information and the sequence encoding used in inputs for a multilayered cascaded graph convolutional network to learn the node features of the sequence encoding and achieve endtoend protein function prediction.Compared with earlier methods，EGNet achieves better performance on the EC category labels in the ENZYME database，which reaches 0.89 in the Fscore and 0.91 in the AUPR.The results indicate that EGNet can achieve good performance by using only a single sequence to predict protein function，providing a rapid and effective method for protein function annotation.

Key words：protein function；deep learning；coevolutionary information；language model；graph convolutional neural network

0 引言

蛋白質(zhì)是一種功能多樣性的生物大分子，是生命的物質(zhì)基礎(chǔ)［1］。隨著高效、低成本測(cè)序技術(shù)和計(jì)算方法的進(jìn)步，UniProt （UniProtKB，2022/04）［2］等關(guān)鍵蛋白質(zhì)序列數(shù)據(jù)庫(kù)中的可用序列數(shù)量大幅增長(zhǎng)，而大多數(shù)蛋白質(zhì)并未進(jìn)行功能注釋。根據(jù)SWISSPROT［3］，在人工注釋的蛋白質(zhì)中，約一半的蛋白質(zhì)是酶，通常使用酶委員會(huì)（EC）系統(tǒng)對(duì)如此大量的酶進(jìn)行分類。

酶作為生物體內(nèi)必不可少的蛋白質(zhì)，能夠催化生物體內(nèi)的化學(xué)反應(yīng)，在生命活動(dòng)中發(fā)揮重要作用。對(duì)酶功能進(jìn)行準(zhǔn)確標(biāo)注，在宏基因組學(xué)、工業(yè)生物技術(shù)以及酶缺乏性疾病的診斷等領(lǐng)域都有著廣泛的應(yīng)用。某些酶的功能障礙會(huì)引起嚴(yán)重的代謝疾病。例如，DNA修復(fù)酶能夠識(shí)別和糾正DNA中的物理?yè)p傷，DNA修復(fù)酶的缺乏會(huì)導(dǎo)致突變的積累，從而可能進(jìn)一步導(dǎo)致各種癌癥［4］。通過(guò)對(duì)酶功能進(jìn)行標(biāo)注，可以更加深入地了解生物體內(nèi)各種生物過(guò)程的機(jī)制和調(diào)控途徑，同時(shí)還有助于發(fā)現(xiàn)新的藥物靶點(diǎn)和治療方法。受驗(yàn)證蛋白質(zhì)功能實(shí)驗(yàn)的規(guī)模、設(shè)計(jì)和成本因素的影響，通過(guò)實(shí)驗(yàn)得到相對(duì)準(zhǔn)確的蛋白質(zhì)功能標(biāo)注是非常困難的。所以，對(duì)于蛋白質(zhì)功能的研究是當(dāng)今的熱點(diǎn)話題。

傳統(tǒng)的酶功能預(yù)測(cè)方法，通常利用蛋白質(zhì)的序列或結(jié)構(gòu)域的同源相似性比對(duì)來(lái)實(shí)現(xiàn)蛋白質(zhì)功能預(yù)測(cè)。使用的工具有BLAST（basic local alignment search tool）［5］和FunFams（functional families）［6］ 。BLAST是一種廣泛使用的序列比對(duì)算法，常用于對(duì)生物學(xué)序列進(jìn)行比對(duì)以確定它們之間的相似性。BLAST可以在數(shù)據(jù)庫(kù)中快速查找相似序列，從而幫助識(shí)別和注釋未知的蛋白質(zhì)功能。FunFams是一種生物信息學(xué)工具，使用聚類算法將結(jié)構(gòu)域聚集為具有相似功能的族群，之后通過(guò)具有相似結(jié)構(gòu)域的蛋白質(zhì)之間的相似性來(lái)推斷其功能。但基于序列相似性或結(jié)構(gòu)域比對(duì)的方法忽略了在生物活動(dòng)中一種蛋白需要與其他蛋白或者配體相互作用來(lái)完成特定的生物學(xué)功能［7］，所以使用同源性方法進(jìn)行預(yù)測(cè)時(shí)會(huì)缺失部分功能。

考慮到序列比對(duì)方法的局限性，研究者們開(kāi)始采用機(jī)器學(xué)習(xí)和深度學(xué)習(xí)方法來(lái)實(shí)現(xiàn)蛋白質(zhì)功能預(yù)測(cè)。一些基于殘基耦合關(guān)系從蛋白質(zhì)相互作用網(wǎng)絡(luò)（proteinprotein interaction network，PPIN）提取功能信息［8～10］，被廣泛用于預(yù)測(cè)酶委員會(huì)（EC）標(biāo)注。DEEPre［11］采用了PSIBlast［12］、HMMER［13］等工具進(jìn)行序列相似性比對(duì)，同時(shí)以位置特定評(píng)分矩陣（PSSM）特征信息作為深度學(xué)習(xí)網(wǎng)絡(luò)的特征輸入，進(jìn)而實(shí)現(xiàn)EC功能預(yù)測(cè)。Dalkiran等人［14］提出了ECPred，該模型實(shí)現(xiàn)了酶功能分層預(yù)測(cè)，每個(gè)EC號(hào)都具有獨(dú)立的機(jī)器學(xué)習(xí)模型。Ryu等人［15］將序列和其對(duì)應(yīng)的催化位點(diǎn)注釋信息送入深度殘差卷積神經(jīng)網(wǎng)絡(luò)（ResidueCNN），實(shí)現(xiàn)酶功能預(yù)測(cè)。Sarker 等人［16］提出的GrAPFI采用一種基于結(jié)構(gòu)域相似圖的酶功能預(yù)測(cè)方法，通過(guò)將蛋白質(zhì)表示為圖形結(jié)構(gòu)來(lái)捕捉蛋白質(zhì)的拓?fù)湫畔⒑拖嗷プ饔藐P(guān)系，并將其輸入到一個(gè)多層神經(jīng)網(wǎng)絡(luò)中進(jìn)行訓(xùn)練，之后執(zhí)行基于鄰域的標(biāo)簽傳播來(lái)進(jìn)行功能注釋。文獻(xiàn)［17］進(jìn)一步采用帶有功能標(biāo)記的有向無(wú)環(huán)圖實(shí)現(xiàn)蛋白質(zhì)功能預(yù)測(cè)。文獻(xiàn)［18］通過(guò)對(duì)構(gòu)建的PPINs進(jìn)行聚類，實(shí)現(xiàn)功能預(yù)測(cè)。Sanderson等人［19］提出了ProteInfer蛋白質(zhì)功能預(yù)測(cè)方法。首先，ProteInfer對(duì)蛋白質(zhì)氨基酸序列進(jìn)行onehot編碼，并通過(guò)深度殘差卷積得到序列嵌入特征，最后經(jīng)過(guò)帶有激活函數(shù)的線性層來(lái)輸出功能類別概率分布。但其仍然采用了BLAST［5］序列對(duì)比工具來(lái)提高預(yù)測(cè)精度，通過(guò)將BLAST得到的頂部命中相關(guān)聯(lián)序列的分?jǐn)?shù)和ProteInfer預(yù)測(cè)的分?jǐn)?shù)相乘，得到最終預(yù)測(cè)結(jié)果。

近年來(lái)，國(guó)內(nèi)外學(xué)者在蛋白質(zhì)功能預(yù)測(cè)研究上不斷探索，基于不同的方法實(shí)現(xiàn)了蛋白質(zhì)功能預(yù)測(cè)。雖然以上方法在蛋白質(zhì)EC功能預(yù)測(cè)中取得了一定的成果，但這些方法在網(wǎng)絡(luò)特征的提取上多采用基于序列或結(jié)構(gòu)域?qū)Ρ确椒?，其在全局搜索過(guò)程中時(shí)間成本過(guò)高，也存在數(shù)據(jù)信息的局限性。當(dāng)然，以上結(jié)果也表明，采用深度學(xué)習(xí)實(shí)現(xiàn)蛋白質(zhì)功能預(yù)測(cè)是可行的。

本文提出了EGNet網(wǎng)絡(luò)模型，采用蛋白質(zhì)預(yù)訓(xùn)練語(yǔ)言模型ESM2［20］和onehot編碼得到蛋白質(zhì)序列編碼；利用蛋白質(zhì)序列殘基之間的協(xié)同進(jìn)化信息，包括序列成對(duì)殘基之間的耦合關(guān)系（paired interaction，PI）和強(qiáng)耦合關(guān)系（strong paired interaction，SPI）矩陣作為圖網(wǎng)絡(luò)的鄰接矩陣，構(gòu)建出圖網(wǎng)絡(luò)的拓?fù)浣Y(jié)構(gòu)，并通過(guò)學(xué)習(xí)序列編碼節(jié)點(diǎn)特征實(shí)現(xiàn)端到端的蛋白質(zhì)功能預(yù)測(cè)。EGNet在預(yù)測(cè)過(guò)程中不需要其他輸入信息，具有預(yù)測(cè)精度高、速度快的優(yōu)勢(shì)，為研究人員提供了便捷、準(zhǔn)確的蛋白質(zhì)功能預(yù)測(cè)方法。

1 研究方法

1.1 網(wǎng)絡(luò)架構(gòu)

如圖1所示， EGNet模型主要包含序列編碼模塊、協(xié)同進(jìn)化信息（coevolutionary information）模塊、兩組串聯(lián)的圖卷積神經(jīng)網(wǎng)絡(luò)［21］（GCN）模塊和MLP分類模塊。

對(duì)于EGNet模型的預(yù)測(cè)過(guò)程主要包含以下幾個(gè)步驟：

a）序列編碼嵌入。給定一條蛋白質(zhì)序列，分別采用ESM2蛋白質(zhì)預(yù)訓(xùn)練語(yǔ)言模型和onehot編碼得到兩組序列編碼特征，具體方法將在1.2節(jié)中詳細(xì)闡述。之后將這兩種編碼信息融合得到最終的序列編碼嵌入，其大小為L(zhǎng)×1306，L為序列長(zhǎng)度。

b）提取協(xié)同進(jìn)化信息。序列經(jīng)過(guò)協(xié)同進(jìn)化信息提取模塊，得到序列殘基間的耦合關(guān)系矩陣PI，其大小為L(zhǎng)×L。為了多尺度地提取協(xié)同進(jìn)化特征，對(duì)PI進(jìn)行了物理特征提?。?2］，得到強(qiáng)耦合關(guān)系矩陣SPI，大小為m×m（m≤L）。對(duì)于PI和SPI的具體提取方法，將在1.3.1節(jié)和1.3.2節(jié)中詳細(xì)闡述。

c）圖卷積神經(jīng)網(wǎng)絡(luò)特征學(xué)習(xí)。將序列水平的編碼特征經(jīng)過(guò)線性層降維，作為圖卷積網(wǎng)絡(luò)的節(jié)點(diǎn)特征，PI和SPI分別作為兩組圖卷積網(wǎng)絡(luò)的鄰接矩陣，進(jìn)行蛋白質(zhì)殘基間的耦合特征學(xué)習(xí)。具體方法在1.4節(jié)中詳細(xì)闡述。

d）節(jié)點(diǎn)池化和功能預(yù)測(cè)。將最后一個(gè)圖卷積網(wǎng)絡(luò)層的節(jié)點(diǎn)輸出進(jìn)行加和池化操作，并采用多層感知機(jī)（MLP）和softmax函數(shù)對(duì)預(yù)測(cè)輸出進(jìn)行打分，實(shí)現(xiàn)蛋白質(zhì)EC功能類標(biāo)的預(yù)測(cè)。式（1）為softmax函數(shù)表達(dá)式。

其中：xi表示輸出類別向量中的第i個(gè)元素；n表示輸出類別向量的維度。

1.2 基于蛋白質(zhì)語(yǔ)言模型ESM2和onehot的序列編碼

ESM2是一種基于BERT［23］模型的蛋白質(zhì)語(yǔ)言模型，包含33個(gè)注意力層［24］，通過(guò)多頭注意力機(jī)制和前向傳播神經(jīng)網(wǎng)絡(luò)，ESM2能夠表示出序列殘基間的耦合關(guān)系，并提供有關(guān)蛋白質(zhì)殘基間遠(yuǎn)程同源性的信息。ESM2在訓(xùn)練過(guò)程中采用了大規(guī)模的無(wú)標(biāo)簽蛋白質(zhì)序列數(shù)據(jù)，包括UniRef50、UniRef90［25］、Pfam［26］、TrEMBL［27］等數(shù)據(jù)庫(kù)，覆蓋了廣泛的物種和功能類別。通過(guò)如此龐大的信息學(xué)習(xí)，可以綜合考慮多個(gè)層面的生物信息，如進(jìn)化同源性等。因此，通過(guò)ESM2對(duì)蛋白質(zhì)序列進(jìn)行編碼，能夠更好地獲得序列間的同源信息特征，有利于對(duì)蛋白質(zhì)功能的研究。

蛋白質(zhì)的結(jié)構(gòu)和功能可能因自然序列的進(jìn)化多樣性而存在相似性。為了提高蛋白質(zhì)功能預(yù)測(cè)的準(zhǔn)確性，EGNet采用了超過(guò)2.5億條訓(xùn)練序列樣本的蛋白質(zhì)語(yǔ)言模型ESM2對(duì)輸入序列進(jìn)行編碼。

EGNet采用了兩種序列編碼方式，包括ESM2和onehot序列編碼。

a）ESM2序列編碼。該編碼是一種用于蛋白質(zhì)序列嵌入的技術(shù)，其輸出為定長(zhǎng)的向量表示，每個(gè)向量表示包含序列中每個(gè)氨基酸殘基的信息。這種序列編碼嵌入技術(shù)能夠捕獲蛋白質(zhì)序列的全局和局部特征。具體而言，對(duì)于長(zhǎng)度為L(zhǎng)的蛋白質(zhì)序列qL，ESM2序列編碼輸出的大小為RL×1280。

qL→ESM2→RL×1280（2）

b）onehot序列編碼。該編碼是蛋白質(zhì)序列分析中常用的一種編碼技術(shù)，能夠?qū)被嵝蛄械奈恢帽硎緸楣潭ㄩL(zhǎng)度的向量，并且僅有一個(gè)元素為1，其余元素為0，如丙氨酸（A）編碼為［1，0，0，…，0，0，0］。EGNet采用26維的onehot編碼方式，對(duì)輸入序列進(jìn)行編碼，其中包括20種常見(jiàn)氨基酸和6種非常見(jiàn)氨基酸。這些非常見(jiàn)氨基酸在蛋白質(zhì)中的出現(xiàn)頻率較低，但在某些生物體中具有特殊的生物學(xué)功能，因此被納入編碼過(guò)程中。對(duì)于長(zhǎng)度為L(zhǎng)的蛋白質(zhì)序列qL，經(jīng)onehot編碼，得到大小為RL×26的編碼特征。

qL→onehot→RL×26（3）

1.3 協(xié)同進(jìn)化信息模塊

1.3.1 耦合關(guān)系矩陣

EGNet的耦合關(guān)系矩陣（PI）主要提取流程如圖2所示。

首先，一條序列經(jīng)過(guò)序列自注意力模塊，得到序列殘基間的自注意力輸出。自注意力（selfattention）是注意力機(jī)制的一種類型，又稱內(nèi)部注意力機(jī)制，其核心思想是在序列數(shù)據(jù)中對(duì)每個(gè)位置進(jìn)行關(guān)注，計(jì)算某一位置與其他位置之間的關(guān)聯(lián)程度。對(duì)于蛋白質(zhì)序列，可以確定每個(gè)殘基與序列中其他殘基的相關(guān)性。因此，采用自注意力機(jī)制能夠?qū)W習(xí)序列殘基間的關(guān)系。

如圖3所示，序列殘基間的自注意力操作主要包含以下幾個(gè)步驟。其中，網(wǎng)絡(luò)參數(shù)源自于ESM2預(yù)訓(xùn)練參數(shù)。

a）序列嵌入。將序列編碼嵌入到向量空間Rd，d為序列編碼維度，并將每個(gè)殘基i的d維特征向量編碼看做一個(gè)查詢向量Qi，其他殘基編碼為一個(gè)鍵向量Kj和一個(gè)值向量Vj，計(jì)算公式為

Qi=WqRid，Kj=WkRjd，Vj=WvRdj（4）

其中：Wq∈Euclid Math TwoRApd×d′，Wk∈Euclid Math TwoRApd×d′，Wv∈Euclid Math TwoRApd×d′是學(xué)習(xí)得到的權(quán)重矩陣，d′是自注意力機(jī)制中的超參數(shù)，即Q、K和V的維度。

b）計(jì)算殘基間的相似度S，其計(jì)算公式為

其中：d為序列編碼的維度。

c）對(duì)殘基i的相似度Si采用softmax函數(shù)進(jìn)行打分，并和值向量V進(jìn)行點(diǎn)積計(jì)算后進(jìn)行加權(quán)求和，得到殘基i的注意力輸出向量Att。

之后，將自注意力模塊得到的注意力特征送入線性層進(jìn)行特征提取。對(duì)于第i個(gè)注意力層的輸出向量Atti，經(jīng)過(guò)線性層變換，得到第i層特征表示向量PIi。linear為全連接神經(jīng)網(wǎng)絡(luò)。

PIi=linear（Atti）（7）

最后，進(jìn)行多層信息融合。如式（8）所示，通過(guò)對(duì)每層的注意力輸出PIi進(jìn)行均值化處理，得到最終的耦合關(guān)系矩陣PI。

其中：m為所采用的注意力層的個(gè)數(shù)。

由于自注意力機(jī)制能學(xué)習(xí)序列成員間的相關(guān)程度，所以，得到的PI矩陣可表示出序列中每個(gè)殘基之間的耦合關(guān)系。

1.3.2 強(qiáng)耦合關(guān)系矩陣

在蛋白質(zhì)序列中，殘基之間的相互作用可以影響其折疊和功能。本文提出強(qiáng)耦合關(guān)系矩陣（SPI）來(lái)描述蛋白質(zhì)殘基之間較強(qiáng)的相互作用關(guān)系。

為了得到殘基間的強(qiáng)相互作用關(guān)系，本文實(shí)現(xiàn)了Cheung等人［22］提出的RC（residue communities）求解方法，將相關(guān)性較高的殘基分為不同的社區(qū)組。通過(guò)對(duì)PI進(jìn)行計(jì)算處理得到SPI，主要流程如圖4所示。

具體計(jì)算步驟如下：

a）為了使SPI在統(tǒng)計(jì)上盡可能獨(dú)立，提取PI矩陣（RPI）的特征值Eij和特征向量vk，表示為

RPIvk=Eijvk（9）

采用Eij的前五個(gè)特征值中的兩個(gè)及其對(duì)應(yīng)的特征向量vk定義前三個(gè)殘基群落。其中，閾值ε=0.05，用來(lái)提取相關(guān)性系數(shù)大的殘基：

（a）社區(qū)Ⅰ（紅），第i個(gè)位置上的殘基滿足

vik=2>max（vik=4，ε）（10）

（b）社區(qū)Ⅱ（藍(lán)），第i個(gè)位置上的殘基滿足

vik=2<-max（vik=4，ε）（11）

（c）社區(qū)Ⅲ（綠），第i個(gè)位置上的殘基滿足

vik=4>max（vik=2，ε）（12）

經(jīng)過(guò)以上數(shù)值計(jì)算，將每條序列上具有強(qiáng)相關(guān)性的殘基分為了三個(gè)社區(qū)，分別表示不同的相關(guān)性強(qiáng)度（紅>藍(lán)>綠）。之后，根據(jù)殘基編號(hào)，將得到的兩兩殘基的強(qiáng)耦合關(guān)系矩陣Rm×m插值到L×L大小的零矩陣，得到最終的強(qiáng)耦合關(guān)系矩陣SPI。

1.4 圖卷積神經(jīng)網(wǎng)絡(luò)模塊

圖卷積神經(jīng)網(wǎng)絡(luò)（GCN）模塊由兩個(gè)串聯(lián)堆疊的GCN塊組成，每個(gè)網(wǎng)絡(luò)塊包含三個(gè)圖卷積層，每個(gè)GCN模塊的內(nèi)部結(jié)構(gòu)如圖5所示。

對(duì)于EGNet模型的每個(gè)GCN層，采用無(wú)向圖G={V，E，A}來(lái)表示，其中節(jié)點(diǎn)V由L個(gè)殘基組成，每一個(gè)殘基對(duì)應(yīng)GCN圖中的一個(gè)節(jié)點(diǎn)；A為鄰接矩陣（包括PI和SPI），定義邊緣E。如果兩個(gè)殘基ai和aj之間存在關(guān)聯(lián)，則鄰接矩陣Aij=1；否則，Aij=0，鄰接矩陣直接體現(xiàn)在PI和SPI矩陣中，圖6（a）（b）展示了蛋白質(zhì)序列A0A073CEA3的PI和SPI作為鄰接矩陣的示例，深色為兩個(gè)殘基在GCN中存在聯(lián)系。GCN的逐層正向傳播定義如式（13）所示。

GCN模塊的初始輸入特征包括PI、SPI和序列編碼嵌入特征。每個(gè)網(wǎng)絡(luò)塊的三層圖卷積層的隱藏層大小均設(shè)為512，輸出為三個(gè)圖卷積層的輸出并集，大小為1 536。其中，一個(gè)GCN塊使用耦合矩陣PI作為鄰接矩陣，從序列編碼嵌入中收集信息。序列編碼特征由每個(gè)殘基的氨基酸類型的onehot編碼（L×26）和ESM2的序列編碼嵌入（L×1280）組成，后經(jīng)過(guò)線性層將序列編碼1 306維度嵌入到1 024維度，作為該網(wǎng)絡(luò)塊的節(jié)點(diǎn)特征。另一個(gè)GCN塊采用強(qiáng)耦合關(guān)系矩陣SPI作為鄰接矩陣，將功能上重要的殘基與殘基群體區(qū)分開(kāi)來(lái)，學(xué)習(xí)對(duì)蛋白質(zhì)功能有顯著貢獻(xiàn)的殘基信息，采用上一圖卷積塊的輸出作為節(jié)點(diǎn)信息。通過(guò)采用兩種鄰接矩陣，EGNet能夠從多個(gè)角度和層面綜合挖掘蛋白質(zhì)序列的特征，提高蛋白質(zhì)功能預(yù)測(cè)的準(zhǔn)確性。

2 結(jié)果與分析

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

在生物化學(xué)領(lǐng)域，Enzyme Commission（EC）編號(hào)是用于對(duì)酶進(jìn)行分類的本體系統(tǒng)，由國(guó)際生物化學(xué)聯(lián)合會(huì)命名委員會(huì)（Nomenclature Committee of the International Union of Biochemistry）制定。每個(gè)EC編號(hào)由四位數(shù)字構(gòu)成，中間用句點(diǎn)分隔，例如，EC：3.1.3.16代表蛋白質(zhì)絲氨酸/蘇氨酸磷酸酶。EC編號(hào)系統(tǒng)包含四個(gè)級(jí)別。在第一級(jí)別中，系統(tǒng)將酶分為氧化還原酶、轉(zhuǎn)移酶、水解酶、裂解酶、異構(gòu)酶和連接酶六大類。EC編號(hào)的第一位數(shù)字表示酶屬于六大類中的哪一類，第二位數(shù)字表示子類，第三位數(shù)字表示亞類，第四位數(shù)字表示酶的底物［28］。目前，EC系統(tǒng)是公認(rèn)的生物數(shù)據(jù)庫(kù)中酶的注釋方式。

采用ENZYME數(shù)據(jù)庫(kù)（2014年12月）［29］的EC類標(biāo)和UniProt （UniProtKB，2022/04）對(duì)應(yīng)的蛋白質(zhì)序列來(lái)構(gòu)建數(shù)據(jù)集。構(gòu)建數(shù)據(jù)集的主要流程如下：

a）剔除了殘基數(shù)小于50的序列以避免片段數(shù)據(jù)，并將最大殘基數(shù)限制在1 024以內(nèi)。

b）為了減弱樣本不均衡帶來(lái)的影響，只保留對(duì)應(yīng)序列數(shù)大于4的EC功能類標(biāo)，其數(shù)據(jù)分布如圖7所示，共剔除了2 503個(gè)功能類標(biāo)，但整體序列樣本僅減少了1.63％。

c）將篩選后的序列采用MMseqs2［30］以50％的序列同一性進(jìn)行聚類，并從每個(gè)聚類中選擇至少一個(gè)帶有注釋功能的代表性蛋白質(zhì)序列。

經(jīng)過(guò)以上步驟，確定了本文所采用的基準(zhǔn)數(shù)據(jù)集，包括29 748條蛋白質(zhì)序列和2 763個(gè)EC功能類別。其中，一條序列對(duì)應(yīng)一個(gè)或多個(gè)功能類別，標(biāo)記不同的催化反應(yīng)。后將數(shù)據(jù)集劃分為訓(xùn)練集（23 742）、驗(yàn)證集（2 617）和測(cè)試集（3 389）三組。在三組數(shù)據(jù)中，每條蛋白質(zhì)序列同樣對(duì)應(yīng)一個(gè)或多個(gè)EC功能類別，其構(gòu)建的數(shù)據(jù)集示例如表1所示。

2.2 超參設(shè)置

在網(wǎng)絡(luò)模型訓(xùn)練過(guò)程中，將訓(xùn)練批次batch的大小設(shè)為64，最多訓(xùn)練100個(gè)epoch。訓(xùn)練過(guò)程采用Adam［31］優(yōu)化器，學(xué)習(xí)率設(shè)置為1×10-4，并加入了權(quán)重衰減操作，衰減率為1×10-4。為了避免過(guò)擬合，訓(xùn)練采用了early stopping操作，并在MLP層引入了dropout操作，參數(shù)設(shè)置為0.3。采用PyTorch來(lái)搭建深度學(xué)習(xí)網(wǎng)絡(luò)框架，訓(xùn)練采用單張GTX2080Ti顯卡，訓(xùn)練成本低。

2.3 評(píng)價(jià)指標(biāo)

本文采用損失函數(shù)值（loss）、召回率（recall）、精確率（precision）、Fscore作為模型泛化能力評(píng)價(jià)指標(biāo)，具體計(jì)算公式如下：

a）對(duì)于損失函數(shù)，采用了交叉熵?fù)p失。

b）召回率。

c）精確度。

d）Fscore。采用β=1，即F1score來(lái)評(píng)價(jià)模型。

為了更準(zhǔn)確地描述評(píng)價(jià)模型的召回率、精確率以及Fscore等指標(biāo)，引入了混淆矩陣的概念?；煜仃囍?，TP表示正例預(yù)測(cè)為正例的個(gè)數(shù)；FP表示負(fù)例預(yù)測(cè)為正例的個(gè)數(shù)；FN表示正例預(yù)測(cè)為負(fù)例的個(gè)數(shù)；TN表示負(fù)例預(yù)測(cè)為負(fù)例的個(gè)數(shù)。

2.4 模型性能對(duì)比分析

對(duì)于EGNet模型，實(shí)現(xiàn)了端到端的預(yù)測(cè)方式，可以直接預(yù)測(cè)蛋白質(zhì)的四級(jí)EC功能標(biāo)注，同時(shí)也可以預(yù)測(cè)多功能蛋白。對(duì)于測(cè)試集的3 389條序列，將EGNet和其他幾種蛋白質(zhì)EC功能預(yù)測(cè)方法，包括BLAST、GrAPFI、FunFams、ECPred、DeepEC和ProteInfer的性能進(jìn)行比較。其中，根據(jù)ProteInfer的研究?jī)?nèi)容，最佳預(yù)測(cè)結(jié)果綜合運(yùn)用了BLAST序列比對(duì)工具，該方法將BLAST所得的與頂部命中相關(guān)聯(lián)序列的分?jǐn)?shù)與ProteInfer預(yù)測(cè)分?jǐn)?shù)相乘，以得到最終的預(yù)測(cè)分?jǐn)?shù)，下述結(jié)果采用了該策略。

對(duì)于召回率（recall）和精確率（precision）評(píng)價(jià)結(jié)果，體現(xiàn)在PR曲線上。由于PR曲線的兩個(gè)指標(biāo)都聚焦于正例，對(duì)于蛋白質(zhì)功能預(yù)測(cè)的研究，主要傾向于預(yù)測(cè)的準(zhǔn)確性，即主要關(guān)心正例的預(yù)測(cè)結(jié)果。因此，采用了PR曲線來(lái)評(píng)價(jià)模型性能，結(jié)果如圖8所示。同時(shí)，也采用了Fscore作為性能評(píng)價(jià)指標(biāo)，表2中列出了每種方法的Fscore和AUPR值，AUPR為PR曲線下的面積。圖9是結(jié)果數(shù)據(jù)的可視化。由結(jié)果可得，EGNet明顯優(yōu)于其他模型。

2.5 PI和SPI不同閾值下的網(wǎng)絡(luò)性能

鄰接矩陣作為GCN的重要輸入，對(duì)其網(wǎng)絡(luò)性能具有一定的影響。為了更好地優(yōu)化模型參數(shù)，提高模型的性能，采用雙循環(huán)閾值設(shè)置（［0.1，0.4］）來(lái)提高PI和SPI對(duì)預(yù)測(cè)模型的正向影響。EGNet在不同閾值下的預(yù)測(cè)結(jié)果如表3所示。其中，當(dāng)PI閾值設(shè)為0.2，SPI閾值設(shè)為0.3時(shí)，EGNet的性能最佳，F(xiàn)score達(dá)到了0.89，AUPR達(dá)到了0.91。

2.6 消融實(shí)驗(yàn)

由表3結(jié)果可知，GCN模塊鄰接矩陣的閾值選擇對(duì)于EGNet性能存在影響。為了驗(yàn)證強(qiáng)耦合關(guān)系矩陣SPI對(duì)于網(wǎng)絡(luò)性能提升的作用，下面對(duì)GCN只采用PI作為鄰接矩陣來(lái)進(jìn)行實(shí)驗(yàn)，進(jìn)而分析SPI對(duì)于網(wǎng)絡(luò)性能的影響。去除SPI前后的實(shí)驗(yàn)對(duì)比結(jié)果如表4所示，訓(xùn)練的超參設(shè)置保持一致。實(shí)驗(yàn)對(duì)PI采用不同閾值進(jìn)行測(cè)試。其中，EGNet采用了PI和SPI，所得的Fscore和AUPR為當(dāng)前PI閾值設(shè)置下的最好結(jié)果。由實(shí)驗(yàn)結(jié)果可得，相較于單獨(dú)采用PI，同時(shí)采用PI和SPI的網(wǎng)絡(luò)性能更好。

3 結(jié)束語(yǔ)

本文提出了一種基于協(xié)同進(jìn)化信息和深度學(xué)習(xí)的蛋白質(zhì)功能預(yù)測(cè)模型。相較于傳統(tǒng)的序列比對(duì)方法，EGNet實(shí)現(xiàn)了端到端的設(shè)計(jì)，通過(guò)處理ESM2的輸出獲取序列間的編碼特征和協(xié)同進(jìn)化信息，并將其作為兩個(gè)串聯(lián)的GCN模塊的節(jié)點(diǎn)特征和邊緣特征，用于訓(xùn)練并實(shí)現(xiàn)蛋白質(zhì)功能預(yù)測(cè)。通過(guò)在ENZYME數(shù)據(jù)集上進(jìn)行實(shí)驗(yàn)，證明了該方法在蛋白質(zhì)功能預(yù)測(cè)上的可行性。實(shí)驗(yàn)結(jié)果表明，與現(xiàn)有模型BLAST、GrAPFI、FunFams、ECPred、DeepEC和ProteInfer相比，EGNet模型能夠更準(zhǔn)確地預(yù)測(cè)蛋白質(zhì)功能。后續(xù)研究將進(jìn)一步發(fā)掘蛋白質(zhì)序列信息，補(bǔ)充對(duì)蛋白質(zhì)基因本體論GO（gene ontology）的預(yù)測(cè)，提供更全面的蛋白質(zhì)功能預(yù)測(cè)信息。

參考文獻(xiàn)：

［1］Ciechanover A.The ubiquitinproteasome pathway：on protein death and cell life［J］.The EMBO Journal，1998，17（24）：71517160.

［2］The UniProt Consortium.UniProt：the universal protein knowledgebase in 2023［J］.Nucleic Acids Research，2022，51：523531.

［3］Bairoch A，Apweiler R.The SWISSPROT protein sequence database and its supplement TrEMBL in 2000［J］.Nucleic Acids Research，2000，28（1）：45-48.

［4］Wood R D，Mitchell M，Sgouros J，et al.Human DNA repair genes［J］.Science，2001，291（5507）：12841289.

［5］Rahman S A，Cuesta S M，F(xiàn)urnham N，et al.ECBLAST：a tool to automatically search and compare enzyme reactions［J］.Nature Methods，2014，11（2）：171174.

［6］Das S，Lee D，Sillitoe I，et al.Functional classification of CATH superfamilies：a domainbased approach for protein function annotation［J］.Bioinformatics，2015，31（21）：34603467.

［7］Xie Zhongru，Hwang M J.Methods for predicting proteinligand binding sites［J］.Molecular Modeling of Proteins，2015，1215：383-398.

［8］Zhao Bihai，Sai Hu，Li Xueyong，et al.An efficient method for protein function annotation based on multilayer protein networks［J］.Human Genomics，2016，10（1）：article No.33.

［9］Chua H N，Sung W K，Wong L.Exploiting indirect neighbours and topological weight to predict protein function from proteinprotein interactions［J］.Bioinformatics，2006，22（13）：16231630.

［10］Nabieva E，Jim K，Agarwal A，et al.Wholeproteome prediction of protein function via graphtheoretic analysis of interaction maps［J］.Bioinformatics，2005，21（S1）：302310.

［11］Li Yu，Wang Sheng，Umarov R，et al.DEEPre：sequencebased enzyme EC number prediction by deep learning［J］.Bioinformatics，2018，34（5）：760769.

［12］Altschul S F，Madden T L，Schffer A A，et al.Gapped BLAST and PSIBLAST：a new generation of protein database search programs［J］.Nucleic Acids Research，1997，25（17）：3389-3402.

［13］Finn R D，Clements J，Eddy S R.HMMER Web server：interactive sequence similarity searching［J］.Nucleic Acids Research，2011，39（2）：29-37.

［14］Dalkiran A，Rifaioglu A S，Martin M J，et al.ECPred：a tool for the prediction of the enzymatic functions of protein sequences based on the EC nomenclature［J］.BioMed Central Bioinformatics，2018，19（1）：113.

［15］Ryu J Y，Kim H U，Lee S Y.Deep learning enables highquality and highthroughput prediction of enzyme commission numbers［J］.Proceedings of the National Academy of Sciences，2019，116（28）：1399614001.

［16］Sarker B，Ritchie D W，Aridhi S.GrAPFI：predicting enzymatic function of proteins from domain similarity graphs［J］.BioMed Central Bioinformatics，2020，21（1）：115.

［17］吳建盛，唐詩(shī)迪，梅德進(jìn)，等.面向蛋白質(zhì)功能預(yù)測(cè)中有向無(wú)環(huán)圖標(biāo)記結(jié)構(gòu)的多示例多標(biāo)記學(xué)習(xí)［J］.國(guó)防科技大學(xué)學(xué)報(bào)，2022，44（3）：2330.（Wu Jiansheng，Tang Shidi，Mei Dejin，et al.Multiinstance multilabel learning for labels with directed acyclic graph structures in protein function prediction［J］.Journal of National University of Defense Technology，2022，44（3）：2330.）

［18］沈婷婷，劉靜，管驍.基于多層功能結(jié)構(gòu)的谷物蛋白質(zhì)功能預(yù)測(cè)［J］.農(nóng)業(yè)工程學(xué)報(bào)，2023，39（1）：261-268.（Shen Tingting，Liu Jing，Guan Xiao.Prediction of cereal protein function based on multilayer functional structures［J］.Transactions of the Chinese Society of Agricultural Engineering，2023，39（1）：261-268.）

［19］Sanderson T，Bileschi M L，Belanger D，et al.ProteInfer，deep neural networks for protein functional inference［J］.Elife，2023，12：e80942.

［20］Lin Zeming，Akin H，Rao R，et al.Evolutionaryscale prediction of atomiclevel protein structure with a language model［J］.Science，2023，379（6637）：11231130.

［21］Kipf T N，Welling M.Semisupervised classification with graph convolutional networks［EB/OL］.（20170222）.https：//arxiv.org/abs/1609.02907.

［22］Cheung N J，Huang Siyu.Residue communities reveal evolutionary signatures of γδ TCell receptor［EB/OL］.（20221230）［20230313］.https：//doi.org/10.1101/2022.12.29.522230.

［23］Devlin J，Chang Mingwei，Lee K，et al.BERT：pretraining of deep bidirectional transformers for language understanding［C］//Proc of Conference of the North American Chapter of the Association for Computational Linguistics：Human Language Technologies.Stroudsburg，PA：Association for Computational Linguistics，2019：41714186.

［24］Vaswani A，Shazeer N，Parmar N，et al.Attention is all you need［C］//Proc of the 31st International Conference on Neural Information Processing Systems.New York：ACM Press，2017：60006010.

［25］Suzek B E，Huang Hongzhan，McGarvey P，et al.UniRef：comprehensive and nonredundant UniProt reference clusters［J］.Bioinformatics，2007，23（10）：12821288.

［26］Bateman A，Coin L，Durbin R，et al.The Pfam protein families database［J］.Nucleic Acids Research，2004，32（S1）：138141.

［27］Bairoch A，Apweiler R.The SWISSPROT protein sequence data bank and its new supplement TREMBL［J］.Nucleic Acids Research，1996，24（1）：21-25.

［28］CornishBowden A.Current IUBMB recommendations on enzyme nomenclature and kinetics［J］.Perspectives in Science，2014，1（16）：74-87.

［29］Bairoch A.The ENZYME database in 2000［J］.Nucleic Acids Research，2000，28（1）：304-305.

［30］Steinegger M，Sding J.MMseqs2 enables sensitive protein sequence searching for the analysis of massive datasets［J］.Nature Biotechnology，2017，35（11）：10261028.

［31］Kingma D P，Ba J.Adam：a method for stochastic optimization［EB/OL］.（20170130）.https：//arxiv.org/abs/1412.6980.