基于多视图矩阵补全的蛋白受体功能预测

黄玮翔 ,丁 季 ,刘夏栩 ,殷 勤 ,兰闯闯 ,吴建盛*

（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 模型的训练时间较长，下一步将改进归纳矩阵补全模型，在保证性能的同时提高训练速度.