基于多視圖矩陣補全的蛋白受體功能預測

黃瑋翔 ,丁 季 ,劉夏栩 ,殷 勤 ,蘭闖闖 ,吳建盛*

（1.南京郵電大學地理與生物信息學院，南京，210023；2.南京郵電大學通信與信息工程學院，南京，210023）

蛋白受體是細胞信號轉導和基因調控的重要組成部分，也是人類主要的藥物靶點，其中G 蛋白偶聯受體（G Protein Coupled Receptors，GPCRs）占絕大多數.GPCRs 是一類具有七跨膜螺旋的膜蛋白受體［1］，是細胞信號轉導的重要組成部分，可以激活細胞內信號轉導通路，最終激活細胞反應［2］.目前，市場上大約34%的藥物都以GPCRs 作為靶點［3-4］，因此，準確注釋GPCR 蛋白的生物學功能對于理解GPCR 蛋白參與的生理過程及其相關的藥物開發具有重要價值.計算的方法是預測蛋白質生物學功能最常用的一種方式［5］，而且蛋白質的生物學功能有多種描述方法，其中基因本體學（Gene Ontology，GO）［6-7］的使用最廣泛.GO 指用來對基因及其產物的功能進行注釋的本體，它包含三個方面：分子功能（Molecular Function，MF）、生物過程（Biological Process，BP）和細胞成分（Cellular Component，CC）.

過去的研究已經開發了大量的計算方法來預測蛋白質的GO 生物學功能，主要分四類.第一類是基于序列同源性搜索的方法，即通過對已知生物學功能的蛋白進行序列相似性搜索來對目標蛋白質進行功能注釋.2017 年Porfiti et al［8］開發了BAR3.0，描述了一個新的具有嚴格度量的非層次聚類過程，該聚類技術捕獲同源和遠緣蛋白質序列之間共享的基本信息來對未鑒定的蛋白質進行功能注釋.2023 年Yuan et al［9］提出一個基于序列的蛋白質功能預測方法SPROF-GO，通過在序列相似性的同源網絡上使用標簽擴散來提高蛋白質功能預測的準確性.第二類是基于序列組合的方法.2018 年CAFA3［10］冠軍方法GoLabeler［11］被提出，該方法從序列輸入中提取同源信息并將它們集成到一個預測因子中來構建預測模型.2020 年Hong et al［12］指出基于序列的深度學習方法可同時提高蛋白質功能注釋的穩定性、準確性和錯誤發現率.2022 年Lai and Xu［13］提出GAT-GO 方法，利用預測的結構信息和蛋白質序列嵌入，大大改善了蛋白質功能預測.同年，Dhanuka et al［14］提出一種基于深度學習的蛋白質功能預測方法，利用一組自動編碼器以半監督的方式用蛋白質序列進行訓練，得出每個自動編碼器對應的蛋白質功能，該方法可擴展到預測任何數量的具有充分支持的蛋白質序列的功能.第三類是基于結構模板的方法，旨在從結構相似的功能模板推斷目標蛋白的功能.2020 年Swenson et al［15］提出一個端到端可訓練的深度學習模型Pers-GNN，結合圖形表示學習和拓撲數據分析來捕獲蛋白質的局部和全局結構特征.2021 年Smaili et al［16］提出QAUST 方法，利用基 于結構相似性、蛋白質相互作用和功能基序的三個信息源來預測蛋白質功能.2022 年Rojano et al［17］提 出DomFun，通過從各種注釋數據庫獲得的蛋白質結構域和功能之間的關聯來預測給定蛋白質的功能，關聯計算由結構域、蛋白質和功能注釋組成的三方網絡得出.第四類是基于生物網絡的方法，主要基于蛋白質-蛋白質相互作用（PPI）、基因組鄰域和共表達模式等信息.2020 年Gumerov and Zhulin［18］將蛋白質特征和基因鄰域信息與系統發育聯系起來，提出一個基于樹的鄰域和域探索新平臺TREND，使基于進化的蛋白質功能分析更加有效.2021 年Barot et al［19］提出一種基于多物種網絡的深度學習方法NetQuilt，有效集成了PPI 網絡信息和同源性用于蛋白質功能預測.2022 年Jagtap et al［20］提出一種用于蛋白功能分析的生物網絡集成方法BraneMF，有效整合了基因共表達網絡、PPI 網絡、遺傳互作網絡、代謝網絡信息來對蛋白質功能進行分析.2022 年Sengupta et al［21］提出PFP-GO 方法，結合序列相似性、PPI 網絡和結構域預測的信息，并利用功能富集得出GO 術語的共識預測，還可以識別功能活躍的蛋白質.2023 年Wu et al［22］提 出CFAGO 方法，通過多頭注意機制將PPI 網絡和蛋白質生物學屬性結合，用于蛋白質功能預測.

雖然基于計算學的蛋白質生物學功能預測方法已經取得很大進展，但還有些問題需要改進.

（1）蛋白質的GO 功能預測中往往只能得到正樣本，很難得到經過實驗驗證的負樣本，數據中大量的負樣本更適合理解為未標記樣本，因此蛋白質的GO 功能預測實質上更偏向是一個Positive-Unlabeled 學習問題（PU-Learning），而不是傳統的監督學習問題，適合用矩陣補全方法來進行解決.

（2）從蛋白質可以提取各種類型的屬性信息，這些屬性信息都會對蛋白質的功能預測有貢獻.多視圖學習［23-24］可以從不同視圖來融合這些屬性信息，從而更加全面地對蛋白質特征進行描述，提高模型預測性能.

（3）傳統的矩陣補全或機器學習方法往往更多地考慮蛋白質視圖中的樣本信息，容易忽略GO 標記視圖中的信息，在矩陣補全中融合GO 標記空間的視圖信息，有利于提升模型預測性能.

（4）目前存在大量的GPCR 蛋白質和GO 術語的文本描述信息，在矩陣補全中有效融合這些文本信息有利于提升模型的預測性能.

因此，本文提出一種基于多視圖歸納矩陣補全（Multi -View Inductive Matrix Completion，MVIMC）的方法對GPCR 蛋白質的GO 生物學功能進行預測.MVIMC 算法將多視圖特征表示方法與歸納矩陣補全技術進行結合，在模型中加入GPCR 蛋白質的多個視圖信息以及GO 術語的視圖信息.在收集到的所有GPCR 蛋白質的GO 生物學功能數據集上進行測試，結果顯示，MVIMC對于GPCR 蛋白質的分子功能和生物過程的GO生物學功能的預測概率分別達到68%和69%，明顯優于目前最好的矩陣補全方法以及CAFA 蛋白質功能預測比賽中的常用方法.

1 數據集和方法

1.1 數據集首先，從Uniprot 生物數據庫（https：//www.uniprot.org/）下載所有GPCRs 的Fasta 格式序列［25］，用NCBI 的blastclust 程序去冗余（相似度小于90%）［26］，得到最終的樣本數據集.然后，從UniProt 數據庫下載“gene_association.goa_ ref_uniprot”.該文件包含蛋白質具有的生物學功能，文件中P 表示生物學過程，F 表示分子功能.GPCR 蛋白均勻分布于細胞膜上，所以本文不考慮GO 的細胞組分預測.從該文件得到GPCR 蛋白質的分子功能和生物學過程的GO 條目（不考慮Evidence code 為IEA 的）.最后，從基因本體學網站（http：//geneontology.org/page/download-ontology）下載文件“go.obo”，運行網站提供的obo2csv.py 程序，得到文件“go.obo.F.is_a”和“go.obo.P.is_a”，得到GPCR 的GO 條目的父節點GO，即得到GPCR 對應的GO 標記空間.刪除樣本個數特別多和特別少的GO 條目.對于分子功能，最終得到GPCR 蛋白質樣本1167個，GO 標記192 個；對于生物學過程，最終得到GPCR 蛋白 質樣本1277 個，GO 標記1203 個.

1.2 特征本文的多視圖屬性包括GPCR 蛋白的視圖和GO 標記的視圖.GPCR 蛋白的視圖分為GPCR 文本信息和GPCR 結構域信息，其中結構域信息包括三聯氨基酸信息、氨基酸關聯信息、進化信息、二級結構關聯信息、物化信息、無序殘基信息、信號肽信息以及結構域文本信息.GO 標記的視圖從GO 術語的文本信息獲得.

1.2.1 GPCR 蛋白視圖

1.2.1.1 GPCR 文本信息從UniProt 數據庫中提取的GPCR 蛋白的文本信息包括蛋白質名稱（Protein Name）、物種信息（Organism）、分子功能關鍵字信息（Keywords for Molecular Function）和相關文獻的標題信息（Titles of Related Publications），在生物醫學文獻數據庫PubMed 中以“receptor（受體）”為關鍵詞進行搜索，得到100多萬篇文獻.以這些文獻的摘要作為比對數據庫，使用基于深度神經網絡的Word2Vec［27］工具，對GPCR 蛋白質的文本信息進行向量化表示［25］.對于一個GPCR 得到的多個向量，采用多示例學習方法miFV［28］將其轉變為一個示例向量，其維度為84.

1.2.1.2 GPCR 蛋白的結構域特征將去除多余信息的GPCRs 氨基酸序列的文件上傳至NCBI的Batch CD-Search 服務器［26，29］，得到其相關的結構域信息.對于蛋白質的每個結構域，提取以下特征信息.

（1）三聯氨基酸信息.按照其偶極矩和側鏈體積，可將20 種氨基酸分為A，B，C，D，E 和F［30-31］六類.對于每個結構域，計算其三聯體出現頻率（Conjoint_triad）［31］：

其中，a，b，c∈{A，B，…，F}，Mabc和l分別表示樣本三聯氨基酸的個數和長度.最終得到的三聯氨基酸信息特征維度為216.

（2）氨基酸關聯信息.氨基酸關聯信息描述結構域中氨基酸間的相關性.依據上面的六類氨基酸，可以得到氨基酸關聯信息（Amino Acid Correlation，AAC）［26，33］為：

其中，m，n∈{A，B，…，F}；Pm表示第m類氨基酸出現的概率，Pn表示第n類氨基酸出現的概率；Pmm(k) (k∈{2，4，8，16})是聯合概率，表示這兩個氨基酸殘基在序列上間隔的殘基數量.最終得到的氨基酸關聯信息特征維度為144.

（3）進化信息，用psiblast 軟件［33］獲得的結構域的位置特異性得分矩陣（Position-Specific Scoring Matrics，PSSMs）來表示：

對于每個結構域，Consortium 陣包含的元素為42×l，其中，l是氨基酸序列的長度.設定42個元素組成的向量為一個示例，則每個結構域的進化信息PSSM 矩陣為l個示例組成的示例包.采用多示例學習方法miFV［28］將其轉變為一個示例向量，最終得到一個84 維的特征向量.

（4）二級結構關聯信息.蛋白質的二級結構通常包含螺旋、折疊和轉角三種狀態，利用PSIPRED 在線分析工具［34］得到結構域的每個殘基的預測二級結構.計算結構域的二級結構關聯信 息（Secondary Structure Element Correlation，SSC）如下：

其中，m，n∈{H，E，C}，為氨基酸的二級結構；k∈{2，4，8，16}，為兩個氨基酸殘基在序列上間隔的殘基數量.

（5）物化屬性.采用SciDBMaker（SDK）［35］軟件得到結構域的各種物化屬性，并利用logistic 函數進行歸一化處理，其維度為59.

（6）無序殘基信息.采用DISOPRED［36］軟件對結構域的無序殘基信息進行預測，并將結構域中的每個氨基酸殘基的特征向量當作一個示例，采用多示例單示例化方法miFV［28］將其轉換為單個示例向量，其特征維度為84.

（7）信號肽信息.采用SignalP［37］軟件對結構域的信號肽信息進行預測，并將結構域中的每個氨基酸殘基的特征向量當作一個示例，采用多示例學習方法miFV［28］將其轉換為單個示例向量，得到的信號肽信息特征維度為84.

1.2.2 GO 標記視圖以“gene function”為關鍵字搜索PubMed 數據庫，得到約180 萬篇文獻，以這些文獻的摘要作為比對數據庫.從Gene Ontology 網站下載的go.obo 文件中得到每個GO 條目對應的文本描述信息，包括name，def，synonym 信息.最后，利用Word2Vec［27］將GO 條目的文本信息轉化成向量，即每個GO 條目表示為由多個示例向量組成的包，采用多示例學習方法miFV［28］將其轉換為單個示例向量，其特征維度為84.

1.3 MVIMC給定樣本集合X={x1，x2，…，xM}，其中，M表示樣本個數，樣本的視圖特征空間為Tv=(T1v，…，TMv)，其中v=1，…，m，Tiv∈R1×dv(i=1，…，M)表示 第i個樣本在第v個視圖上 的特征向量，dv表示第v個視圖的維度.Y=(y1，…，yN)表示標記集合，N表示標記個數，標記的特征空間為Q=(Q1，…，QN)，其中Qj∈R2c(j=1，…，N)表示第j個標記的特征向量.

假設樣本-標記關系矩陣為S∈RM×N，M表示樣本個數，N表示標記個數.Si，j=1 表示第i個樣本與第j個標記的關系已知，Si，j=0 表示關系未知.MVIMC 算法利用新型的歸納矩陣補全方法，根據觀測到的關系矩陣S在第v個視圖上補全潛在的真實關系矩陣Zv=Wv HvT，其中，Wv，Hv為分解后的兩個子矩陣，Zv∈Rdv×2c，Wv∈Rdv×k，Hv∈R2c×k.假設關系矩陣S是低秩的，MVIMC算法的目標函數如下：

其中，Ω表示具有已知關系的樣本-標記對的集合；l(·)是損失函數，用來衡量預測值與真實值之間的誤差.通常采用均方誤差作為損失函數，即l(a，b)=(a-b)2.第二項為正則化項，用來控制模型復雜度及避免過度擬合，其中參數λ用來平衡損失函數和正則化約束.

式（5）是一個非凸函數，為了求解目標函數，首先，隨機初始化W和H矩陣，然后使用交替最小化（即固定W求解H，再固定H求解W）方法求解W和H矩陣，直至達到收斂或局部最優.具體更新步驟如下.

固定Hv，更新Wv，式（5）等價于：

固定Wv，更新Hv，式（5）等價于：

通過固定求解的方法得到的式（6）和式（7）都是凸函數，可以采用共軛梯度法來求解.解得Wv和Hv后，對于任意的(a，b) ?Ω，都可以通過求解得到，即可預測關系矩陣S中的未知值.

MVIMC 首先得到模型在各個單視圖上的預測結果，對各個視圖上的預測性能進行排序，并對視圖進行不同組合作為模型的特征輸入，不斷進行優化求解，得到最優的視圖組合和模型.

1.4 評價指標采用矩陣補全中常用的預測概率（Probability，P）和相關錯誤率（Relative Error，rel.err）來對模型進行評價.

預測概率P指一個真實的樣本-標記關系對在得分前r位的預測中被發現的概率.P越大，說明預測越準確.

其中，X為預測的關系矩陣，M為真實的關系矩陣.rel.err越小，預測越準確.

采用三倍交叉驗證來評估模型的性能，即在構建模型的過程中，將數據集（即樣本-標記關系對）隨機分為三等份，每次使用二等份進行訓練，剩下的一等份進行測試.重復執行三次，保證每個關系對都被預測一次.

2 實驗與結果

2.1 不同視圖上的性能比較采用歸納矩陣補全算法（Inductive Matrix Completion，IMC）對GPCR 結構域的各個視圖進行單視圖建模，再進行視圖組合，并根據模型性能好壞對各視圖的性能進行排序.結果如圖1 所示.

圖1 不同的單視圖和組合視圖下GPCR 蛋白的GO 功能預測的比較：（a）分子功能；（b）生物學過程Fig.1 Performance of various views：(a) molecular function，(b) biological process

圖1 展示了不同的單視圖和組合視圖下的GPCR 蛋白的GO 功能預測的比較，圖中橫坐標表示預測得分最高的前r位樣本.由圖可見，當r一定時，基于三聯氨基酸信息（A）的GPCR-MF和GPCR-BP 關系矩陣預測模型的性能明顯優于其他視圖.對于分子功能（MF），基于三聯氨基酸信息（A）的預測概率為60%，基于其他視圖的模型性能從高到低的排序為：氨基酸關聯信息（B），進化信息（C），二級結構關聯信息（E），物化屬性（F），GPCR 文本信息（D），無序殘基信息（G），信號肽信息（H），結構域文本信息（I）.對于生物學過程（BP），基于三聯氨基酸信息（A）的預測概率為51%，其他特征模型的性能從高到低的排序分別為：GPCR 文本信息（D），氨基酸關聯信息（B），進化信息（C），二級結構關聯信息（E），物化屬性（F），結構域文本信息（I），無序殘基信息（G），信號肽信息（H）.

圖2 展示了不同的組合實現的GPCR 蛋白的GO 功能預測的比較，圖中橫坐標表示預測得分最高的前r位樣本.由圖可見，r=100 時，對于分子功能（MF），最優視圖組合為A+B+C+E，預測概率近67%.對于生物學過程（BP），最優視圖組合為A+D+B+C+E+F+G，預測概率近68%.以上兩個最優視圖組合的預測性能均優于所有單視圖的預測性能.

圖2 IMC 組合視圖方法的預測概率比較（生物過程）Fig.2 Prediction probabilities of IMC combined view method (biological process)

根據已知的關系矩陣與補全的關系矩陣，得到加入各個不同視圖樣本特征的相關錯誤率，如圖3 所示.由圖可見，當加入的特征信息為三聯氨基酸信息（A）時，模型預測的相關錯誤率最小，但整體上，各個視圖模型的相關錯誤率相差不大.

圖3 各視圖的相關錯誤率比較Fig.3 Relative error among different views

圖4 表明，對于分子功能（MF），當加入的視圖組合信息為A+B+C+E 時，模型預測的性能最好，相關錯誤率為1.0732.對于生物學過程（BP），各種視圖的組合，其模型的相關錯誤率差異不大.

圖4 各組合視圖模型的相關錯誤率比較Fig.4 Relative error among different combined view models

綜上，可以得到基于IMC 的最優視圖組合：對于分子功能（MF），最優視圖組合為A+B+C+E，預測概率近67%.對于生物學過程（BP），最優視圖組合為A+D+B+C+E+F+G，預測概率近68%.說明在蛋白質生物學功能的預測中，采用多視圖方法，組合多個單視圖，多角度加入樣本的特征信息，可以提高模型的預測性能.

2.2 不同矩陣補全算法的比較將本文提出的MVIMC 算法與目前最好的幾種傳統矩陣補全算法進行比較，包括Catapult［38］，katz［38］，ALM［39］，FPCA［40］，LmaFit［41］，SVT［42-43］和Maxide［44］.Catapult 和katz 利用樣本與樣本的關系矩陣以及標記與標記的關系矩陣.ALM 用增廣拉格朗日乘數法來精確恢復被破壞的低秩矩陣.FPCA 采用不動點和Bregman 迭代算法來解決線性約束矩陣秩最小化問題.LmaFit 通過非線性連續過松弛算法解決矩陣完成的低秩因子分解問題.SVT 通過奇異值閾值算法，用最小核范數來近似矩陣完成.所有對比算法均采用文獻中的默認參數.

圖5 展示了MVIMC 與對比算法的預測概率的比較，圖中橫坐標表示預測得分值最高的前r位樣本.由圖可見，對于GPCR 的GO 分子功能和生物學功能過程的關系預測，MVIMC 性能最優.對于GO 分子功能，在前100 位樣本中，MVIMC 的預測概率達到68%，比第二位的的LmaFit 高17%，比Catapult 高27%，比排名最末的FPCA 高52%.對于GO 生物學功能的預測，在前100 位樣本中，MVIMC 的預測概率達69%，比katz 和Catapult 高30% 左 右，比排名最末的ALM 算法高59%.

圖5 不同矩陣補全算法的預測概率的比較Fig.5 Prediction probabilities of different matrix completion algorithms

圖6 展示了不同的矩陣補全算法的相關錯誤率的比較.由圖可見，對GPCR 的GO 分子功能和生物學過程的預測，MVIMC 的相關錯誤率最低，約為1，第二位的Catapult 的相關錯誤率約為1.2，第三位的katz 的相關錯誤率約為1.3，最差的FPCA 的相關錯誤率為9 左右，約為MVIMC的9 倍.證明本文提出的MVIMC 算法對于GPCR 的GO 分子功能和生物學功能過程的關系預測明顯優于其他的矩陣補全方法，這是因為和傳統的矩陣補全方法相比，MVIMC 不僅加入了樣本和標記的特征進行模型訓練，同時還考慮了多視圖的特征，提高了預測性能.

2.3 視圖組合方法的比較對三種經典的視圖組合方法Concate，Max 和Ave_score 進行了比較，其中，Concate 將最優視圖組合連接成一個長向量，Max 使用最優的單個視圖來表示，Ave_score 使用所有視圖的平均值來進行衡量.

圖7 展示了不同的多視圖方法的預測概率的比較，圖中橫坐標表示預測得分值最高的前r位樣本.由圖可見，對于分子功能（MF）和生物過程（BP）的預測，Concate 性能均最優，預測概率分別為68%和69%.

圖7 不同多視圖方法的預測概率比較Fig.7 Prediction probabilities of different multi-view methods

圖8 展示了不同的多視圖方法的相關錯誤率的比較.由圖可見，這三種多視圖方法的相關錯誤率相差不大，均在1 左右.綜上，在相關錯誤率均較低的情況下，Concate 比Max 和Ave_score 具有更高的預測概率.本文提出的MVIMC 算法正是采用了Concate 多視圖方法，將最優視圖組合的特征拼接成一個長向量，然后利用神經網絡學習融合這些特征.在這個過程中信息不會損失，因而該方法的預測概率也就高于其他兩種多視圖方法.

圖8 不同多視圖方法的相關錯誤率比較Fig.8 Relative error of different multi-view methods

2.4 與CAFA 蛋白質功能預測方法的比較功能注釋關鍵評估（Critical Assessment of Functional Annotation，CAFA［10，45］）挑戰是國際上最權威的蛋白質功能注釋比賽，已經舉辦了四屆.CAFA中用于預測蛋白質功能的三種基本方法為Naive［45］，BLAST［45］和PSI-BLAST［45］.為了驗證本文提出的MVIMC 算法對于GPCR 蛋白質功能預測的有效性，將MVIMC 算法與這些方法進行了比較.GO 性能主要通過所有閾值的精度和召回之間的最大諧波平均值來評估，其值越大越好.計算如下：

其中，t是判定閾值，范圍在0～1；pr(t)和rc(t)分別表示閾值t處的精度和召回值.

圖9 展示了本文提出的MVIMC 算法與CAFA 預測平臺的性能比較.由圖可見，對于GPCR的GO 分子功能（MF）和生物學過程（BP）的預測，MVIMC 表現出更好的性能，Fmax分別達到31%和38%，均遠高于其他三種預測蛋白質功能的基本方法，證明了MVIMC 算法對于GPCR 蛋白的GO 功能預測的優越性.

圖9 MVIMC 算法與CAFA 預測平臺的性能比較Fig.9 Performance of MVIMC algorithm and CAFA prediction platform

3 結論

本文提出一種基于多視圖的歸納矩陣補全方法MVIMC，將多視圖表示與歸納矩陣補全技術相結合，并加入樣本多個視圖的特征以及標記的特征信息，實現了對分子功能和生物過程兩方面的GPCR 蛋白的GO 功能預測.在包含1167 個GPCR 的GO 生物學功能預測數據集上進行了測試，結果證明MVIMC 優于目前的矩陣補全算法.對于分子功能的預測，與排名第二的LmaFit相比，在前100 位樣本中，MVIMC 的預測概率平均提升17%；對于生物功能的預測，與排名第二的katz 相比，在前100 位樣本中，MVIMC 的預測概率平均提升29%.MVIMC 算法的預測結果還優于CAFA 挑戰賽中用于預測蛋白質功能的三種基本方法，對于分子功能和生物過程，MVIMC的預測概率分別提高24%和31%.但MVIMC 模型的訓練時間較長，下一步將改進歸納矩陣補全模型，在保證性能的同時提高訓練速度.