化合物-蛋白质对的交互预测任务旨在给出一个化合物分子和一个蛋白质的氨基酸序列, 判断二者是否交互, 实际上是一个二分类问题.形式化定义如下:给定一个化合物分子V=[v₁, v₂, …, v_N]; 一个氨基酸序列S=[x₁, x₂, …, x_M]; 交互的标签集合L={l₁, l₂}, 其中, l₁表示交互的类别标签, l₂表示不交互的类别标签.化合物-蛋白质对的交互预测任务是学习一个分类模型G:(V, S)→ L.模型输入化合物分子V和氨基酸序列S, 输出类别标签.

本文方法的整体流程如图1所示.主要过程如下.1)基于GAT-GRU的化合物分子的表示学习, 实现对化合物分子的特征向量化的表示.首先对化合物序列进行嵌入, 再利用图注意力层^[20]实现对化合物分子特征的提取与表示, 并利用门控循环单元^[21]实现对高阶邻居节点的特征进行抽取与过滤.2)基于ML-SRU的氨基酸序列的特征提取.首先对氨基酸序列进行嵌入保持与化合物的维度一致, 再利用多层简单循环单元提取氨基酸序列的特征并进行向量化表示.3)基于注意力机制的蛋白质向量表示, 拼接前两部分生成的向量表示, 作为网络输入, 通过注意力机制模拟相互作用强度, 形成蛋白质的向量表示, 送入多层前馈网络进行分类预测.

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 本文方法流程图Fig.1 Flow chart of the proposed method

基于GAT-GRU的化合物分子的表示学习模型主要由3个相同的基本模块堆叠而成, 每个模块由图注意力层-标准化-拼接-门控循环单元构成, 将各模块之间的计算结果逐层拼接以整合各模块的特征向量表示, 实现对高阶邻居节点的特征抽取与更新.另外, 考虑到数据经过图注意力卷积操作后其概率分布可能会发生偏移, 导致模型表达能力下降, 因此在进行图注意力卷积操作后引入实例标准化操作(Instance Normalization, IN)^[23], 使新的分布更切合数据的真实分布, 保证模型的非线性表达能力.模型框架如图2所示.

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 基于GAT-GRU的化合物分子的表示学习模型框架Fig.2 Framework of representation learning model for compound molecules based on GAT-GRU

由于网络输入是基于SMILES格式的化合物, 因此使用RdKit软件包将基于SMILES格式的数据转换为图结构的数据, 将化合物中的原子视为图中的节点, 将原子之间的化学键视为图中的边并对其进行建模.利用Rdkit提取相应原子特征信息, 存放于节点特征矩阵X∈ R^{N× C}中, 原子间的结构信息存放于邻接矩阵A∈ R^{N× N}中, X和A能为模型提供每个原子及其相邻原子的特征信息及结构信息.

节点特征矩阵X=[x_ij]_{N× C}, 其中, N为节点数, C为特征向量的维度, X的第i行元素x_i· 表示第i个原子v_i对应的特征向量, X的第j列元素x_{· j}表示各个原子的第j维特征.邻接矩阵A=[a_ij]_{N× N}描述原子之间的结构关系, a_ij表示原子v_i与v_j的连接关系, 当原子v_i、v_j之间有键连接时, a_ij=1, 无键连接时, a_ij=0.

具体实现过程如下.

1)对X进行嵌入操作, 转换为低维向量表示形式X'∈ R^{N× C'}, C'为目标嵌入维度.为了增强节点特征向量的表征能力, 对X'进行线性变换, 即通过线性层学习权重矩阵W₁∈ R^{C'× C'}, 获得新的节点特征矩阵X″=X'W₁.

2)将X″和A一起送入各模块的图注意力层进行图注意力卷积操作, 学习化合物分子图中潜在的特征表示, 经标准化操作后得到特征表示H₀∈ R^{N× C'}.

3)将H₀与未经过图注意力层处理的特征表示X″进行拼接, 送入门控循环单元进行进一步的特征提取与过滤, 得到特征表示H₁∈ R^{N× C'}.

4)将H₁与模块2中以同样方式计算得到的H₀进行拼接, 送入模块2的门控循环单元进行处理.重复此过程, 可依次得到H₂∈ R^{N× C'}, H₃∈ R^{N× C'}, 并对H₃进行线性变换, 使H=H₃W₂, 其中, W₂∈ R^{C'× C'}, H∈ R^{N× C'}.此时H已融合各模块的特征信息, 将H按照特征维度对所有原子在该维度的特征值求平均, 获得该化合物分子的特征向量表示H∈ R^{1× C'}.

具体公式如下:

H₀=norm(GATConv(Linear(x), A)),

H₁₌GRU(concat(H_0, Linear(x))),

H₂=GRU(concat(H₀, H₁)₎,

H₃=GRU(concat(H₀, H₂))_,

H=mean(Linear(H₃)),

其中, H₀表示由图注意力层和标准化操作得到的隐藏表示, H_i(i=1, 2, 3)表示各个模块计算的隐藏特征表示, norm表示标准化操作, GATConv表示图注意力层, GRU表示门控循环单元, concat表示拼接操作, Linear表示线性变换, mean表示对各维度的特征表示求均值, H表示最后得出的化合物分子的特征向量表示.

具体地, GAT使用自注意力机制解决基于图卷积模型的不足.通过注意力机制对化合物分子图中的相邻节点进行图卷积操作, 将相应的权重分配给不同的邻居节点, 并对邻域内节点特征进行加权求和, 充分考虑相邻节点间的结构关系, 将图中各节点特征之间的相关性融入模型中, 有助于模型学习结构信息.GATConv操作利用图的结构信息A和节点特征信息X进行节点向量表示, 将待更新的节点记为v_i, 对应的特征向量记为h_i, 经注意力机制的聚合操作输出新的特征向量表示h'_i, 计算每对节点v_i和v_j间的权重系数:

e_ij=LeakyReLU(a^T(Wh_i‖ Wh_j)),

其中, a为节点特征的参数, W为特征变换的权重.为了对不同邻居节点分配相应权重, 使用softmax对权重进行归一化操作, 归一化后权重系数为

α _ij=softmax(e_ij)= exp(eij)∑vk∈N(vi)exp(eik),

其中N(v_i)表示节点i的邻居节点, 则节点v_i更新后的特征向量

给定氨基酸序列S=[x₁, x₂, …, x_M], 首先利用n-gram方法划分给定的氨基酸序列, 本文使用n=3的划分, 形成连续且重叠的子序列:

[x₁x₂x₃]=s₁,

[x₂x₃x₄]=s_2,

︙

[x_M-2x_M-1x_M]=s_M-2.

在数据处理时将氨基酸的字符序列逐一映射为按其长度排序的索引形成字典, 按照字典对每个氨基酸子序列进行向量化表示.将划分后的字符序列S'=[s₁, s₂, …, s_M-2]按照字典映射为数值形式的向量Q=[q₁, q₂, …, q_M-2], 其中q_i表示第i个氨基酸序列s_i的向量形式.

多数现有方法使用CNN处理氨基酸序列, 但是CNN在处理长距离信息提取时存在困难, 因此考虑使用循环神经网络(Recurrent Neural Network, RNN)代替CNN, 而RNN在处理长序列时存在长距离依赖问题, 原因在于其串行的计算方式容易产生梯度消失或梯度爆炸.长短期记忆网络(Long Short-Term Memory, LSTM)和门控循环单元(Gated Recurrent Units, GRU)以门控方式控制信息传递, 各个门及当前隐状态的计算依赖于上一步隐藏状态的输出.为了降低对门控单元和隐藏状态计算的复杂程度, 并且能在长序列中并行计算, 对简单循环单元(Simple Recurrent Units, SRU)^[22]进行改进, 使各个门和隐状态的计算仅与当前输入有关, 不依赖于上一步的隐藏状态.另外, 在SRU的高速网络中引入跳跃连接, 使梯度直接向前传递, 使模型即便叠加多层, 仍能避免梯度消失并保持模型性能, 因此本文使用堆叠的多层SRU进行特征提取.

基于ML-SRU的氨基酸序列的特征提取模块的具体过程如下.将氨基酸序列嵌入实值向量空间中, 输入至堆叠的ML-SRU中进行特征提取.SRU处理特征信息的过程为t时刻, 读取输入向量q_t, 计算细胞状态序列:

c_t=f_t☉c_t-1+(1-f_t)☉ q~t,

并捕获序列的信息, 其中,

f_t=σ (W_f q~t+b_f),

表示遗忘门用于控制输入信息 q~t=Wq_t的大小.重置门

r_t=σ (W_r q~t+b_r),

用于将当前的输入信息q_t及当前时刻的细胞状态c_t进行组合, 输入前一时刻细胞的平均状态c_t-1及当前输入信息 q~t共同决定当前细胞状态c_t.隐藏层输出状态为

h_t=r_t☉g(c_t)+(1-r_t)☉q_t,

其中, W、W_f、W_b、b_f、b_r为可训练的参数矩阵.由于 q~t、 f_t不依赖于前一时刻的状态c_t-1, 因此具备独立性, 能执行并行运算.基于ML-SRU的氨基酸序列的特征提取模型框图如图3所示.

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 基于ML-SRU的氨基酸序列的特征提取框架Fig.3 Feature extraction framework of amino acid sequence based on ML-SRU

在经过前2个模型的特征提取后, 得到化合物分子的向量表示及氨基酸子序列的向量表示.由于不同的氨基酸子序列与化合物分子的作用强度存在差异, 在计算时需要考虑二者间的相互作用, 因此考虑使用软性注意机制^[24]模拟这种相互作用.注意力机制通过对数据进行加权处理, 突出显示更重要的特征信息.本文使用注意力机制学习各氨基酸片段的重要性, 即某个氨基酸子序列对化合物越重要, 则赋予该子序列越大的权重.使用注意力的分数表示相互作用的强度, 通过对氨基酸序列的加权求和得到蛋白质的向量表示.

给定经ML-SRU特征提取后蛋白质的氨基酸子序列的一组隐藏向量Q=[q₁, q₂, …, q_M-2], 蛋白质向量的表示过程如下所示:

其中:v_mol表示输入GAT-GRU的化合物分子的向量表示; q_i表示输入ML-SRU中的第i个氨基酸子序列; hidden_mol表示经过GAT-GRU处理后的化合物分子的隐藏向量表示, hidden_i表示经过ML-SRU处理后得到的第i个氨基酸子序列的隐藏表示, 将基于hidden_mol和hidden_i的点积计算的α _i作为注意力分数, 通过对hidden_i的加权求和计算得出最终的蛋白质序列的隐藏表示; W表示权重矩阵, b表示偏置; GAT-GRU(· )表示GAT-GRU模型, ML-SRU(· )表示ML-SRU模型.注意力分数是通过在网络结构中加入一个前馈神经网络进行学习的.使用神经注意机制, 可反映出在相互作用中蛋白质的哪些区域是重要区域.

将得到的蛋白质的向量表示与化合物分子的向量表示进行拼接, 输入多层前馈网络, 最后一层网络得到的特征向量可作为分类器所需的特征, 经过softmax激活函数进行激活, 输出化合物-蛋白质对交互的概率:

p_k= softmax(W_outconcat(hidden_mol, hidden_i)+b_out).

其中, W_out表示全连接层的权重矩阵, b_out表示全连接层的偏置.模型使用有监督学习的方式进行学习, 使用反向传播算法更新模型参数, 并利用最小化交叉熵训练模型.为了防止过拟合在交叉熵公式后加上L₂正则项, 交叉熵损失函数表示如下:

Loss=- ∑i=1np^kiln( pki)+ λ2‖ θ ‖ .

其中:k=0表示无交互, k=1表示交互; n为训练集的数据量; p^ki表示待分类数据的预测类别; pki为数据实际的类别; λ 表示正则项系数, 用于权衡正则项与交叉熵函数的比重; θ 表示所有权重矩阵和偏置向量的集合.

实验使用Human和C.elegans^[25] 2个生物领域常用的基准CPIs数据集进行评估, 具体数据集信息如表1所示.

表1

Table 1

表1(Table 1)

表1 实验数据集 Table 1 Experimental datasets 数据集 蛋白质 种类数 化合物 种类数 样本总计 正负样本 比例 Human 852 1052 6728 1:1 C.elegans 2504 1434 7786 1:1

表1 实验数据集 Table 1 Experimental datasets

每个数据集中记录化合物与蛋白质之间有无交互的关系, 每个化合物以SMILES格式表示, 每个蛋白质以氨基酸序列的形式表示, 数据样本中的标签分为有交互与无交互, 分别用1和0表示.

实验中每个数据集均按照80%、10%、10%的比例划分为训练集、验证集、测试集.

超参数设置如下:使用Adam优化器进行训练, 学习率设置为1e-4, 权重衰减设置为1e-6, 学习率衰减为0.9, 嵌入维度为30, 输出层层数及GAT-GRU中的模块数均为3, SRU层数选择为5, 迭代次数为100轮.

实验环境设置如下.操作系统为Windows10, 内存8.0 GB, 处理器Intel(R)Core(TM)i5-5350U CPU@1.8 GHz, 使用的框架是Pytorch, 编程语言为Python.

为了衡量方法性能, 采用如下评价指标, 包括曲线下面积(Area under Curve, AUC), 精确率(Preci-sion), 召回率(Recall), F₁值:

Precision= TPTP+FP,

Recall= TPTP+FN,

F₁= 2×Precision×RecallPrecision+Recall.

其中:TP表示真正例, 实际为正且预测为正; FP表示假正例, 实际为负但预测为正; FN表示假反例, 实际为正但预测为负; TN表示真反例, 实际为负且预测为负.F₁值用于衡量Precision和Recall, 是这两个值的调和均值, 综合反映整体情况.

在泛化能力的实验中使用正确率(Accuracy, ACC)作为评估标准:

Accuracy= TP+TNTP+TN+FP+FN.

为了验证本文方法预测的有效性, 选取目前最具代表性的2种基于深度学习^{[16, 17]}方法和6种基于机器学习^[25]的方法进行实验对比.为了准确复现基准方法, 部分对比实验使用原论文设置的相关参数并在相同数据集上运行该论文中开源代码, 将得到的结果作为实验对比依据.

用于对比的2种基于深度学习方法如下.

1)GraphDTA^[16].基于GNN和CNN的组合, 结合池化操作, 用于计算药物-靶标结合亲和力, 是一种回归任务, 实验时对其进行微调, 改为二分类问题, 用于CPIs预测.

2)GCN-CNN^[17].基于GCN和CNN的组合, 利用GCN提取化合物分子特征, CNN提取蛋白质序列特征, 并利用注意力机制模拟二者相互作用, 用于CPIs的预测.

用于对比的6个机器学习方法如下.

1)朴素贝叶斯(Naive Bayesian, NB).

2)k近邻算法(k-Nearest Neighbor, KNN).

3)随机森林(RandomForest, RF).

4)L1(L1-logistic).

5)L2(L2-logistic).

6)SVM.

实验设置如下:NB、KNN、RF由Weka3.7运行, L1、L2由liblinear1.94运行, SVM由libsvm3.17运行.各分类方法分别在各软件的默认环境下进行.

本文方法和基于深度学习方法在2个数据集上的性能对比如表2所示.由表可看出, 本文方法的预测结果更优.在Human数据集上, 相比GraphDTA、GCN-CNN, 本文方法的AUC值分别提升3.5%和2.7%, F₁值分别提升7.9%和6.5%.在C.elegans数据集上, 相比GraphDTA、GCN-CNN, 本文方法的AUC值分别提升2.2%和1.6%, F₁值分别提升5.2%和3.9%, 从而验证本文方法的有效性.

表2

Table 2

表2(Table 2)

表2 3种方法在2个数据集上的性能对比 Table 2 Performance comparison of 3 methods on 2 datasets 方法 Human C.elegans AUC F1 AUC F1 GraphDTA 0.962 0.906 0.972 0.920 GCN-CNN 0.970 0.920 0.978 0.933 本文方法 0.997 0.986 0.994 0.973

表2 3种方法在2个数据集上的性能对比 Table 2 Performance comparison of 3 methods on 2 datasets

本文方法与机器学习方法的性能对比如表3所示.

表3

Table3

表3(Table3)

表3 7种方法在2个数据集上的性能对比 Table3 Performance comparison of 7 methods on 2 datasets 方法 Human C.elegans AUC Precision Recall AUC Precision Recall NB 0.672 0.624 0.575 0.645 0.613 0.465 KNN 0.860 0.798 0.927 0.858 0.801 0.827 RF 0.940 0.861 0.897 0.902 0.821 0.844 L1 0.908 0.881 0.893 0.882 0.872 0.849 L2 0.911 0.891 0.913 0.892 0.890 0.877 SVM 0.910 0.966 0.950 0.894 0.785 0.818 本文方法 0.997 0.986 0.987 0.994 0.973 0.974

表3 7种方法在2个数据集上的性能对比 Table3 Performance comparison of 7 methods on 2 datasets

由表3可看出, 相比表现最优的RF, 本文方法在Human数据集上AUC值提升5.7%, 在C.elegans数据集上AUC值提升3.2%.这可能是因为本文方法属于端到端的表示学习方法, 可基于数据进行表征学习, 自动学习到有意义的数据特征信息, 从而反映化合物-蛋白质交互性质.而上述机器学习方法以化合物的结构和蛋白质的结构作为输入特征, 利用化合物分子图和蛋白质序列的生物化学特性, 分类效果的优劣很大程度上依赖于人工设计的特征.

为了验证本文方法中各模块的有效性, 进行消融实验.对比方法如下.

1)GCN-ML-SRU.将本文方法中的GAT-GRU模块替换成普通GCN网络, 并与ML-SRU结合.

2)GAT-ML-SRU.将本文方法中的GAT-GRU模块替换成普通GAT网络, 并与ML-SRU相结合.

3)None-attn.在本文方法的基础上在计算蛋白质向量的过程中取消注意力机制, 直接将SRU特征提取后的特征向量作为蛋白质向量表示.

各方法的消融实验结果如表4所示.由表可见, 对比GAT-ML-SRU, 本文方法在Human、C.elegans数据集上的AUC值分别提升1.3%和0.6%, 表示在图注意力卷积操作中引入GRU单元有利于进一步提取高阶邻域信息, 从而验证GRU单元的有效性.对比GAT-ML-SRU、GCN-ML-SRU可发现, GAT-ML-SRU比GCN-ML-SRU的实验结果更优, 表示基于GAT的方法提取节点特征信息的能力更强, 而SRU能缓解CNN难以捕获长距离特征信息问题.对比本文方法与None-attn发现, 当不使用注意力对蛋白质进行表示时, 方法预测性能的评价指标均呈现不同程度的下降, 表示使用注意力机制能有效捕获蛋白质的特征信息, 说明蛋白质向量表示的改进可使方法性能得到提升.

表4

Table 4

表4(Table 4)

表4 各方法在2个数据集上的消融实验结果 Table 4 Ablation experiments of different methods on 2 datasets 方法 Human C.elegans AUC F1 AUC F1 GAT-ML-SRU 0.984 0.956 0.988 0.954 GCN-ML-SRU 0.979 0.943 0.981 0.942 None-attn 0.985 0.973 0.989 0.961 本文方法 0.997 0.986 0.994 0.973

表4 各方法在2个数据集上的消融实验结果 Table 4 Ablation experiments of different methods on 2 datasets

本节分析可能对CPIs预测性能产生影响的因素, 主要从SRU网络层数、嵌入维度、是否使用标准化操作及GAT-GRU中的模块数进行分析.

通常多层神经网络的非线性表征能力较强, 能学习到更复杂的变换, 可拟合更高级的特征, 使网络性能有所改进.SRU网络层数对评价指标的影响如图4所示.由图可看出, 本文方法在层数为5时表现最佳, 但随着层数的增加, 性能并不会一直提升.其原因可能是随着网络层数的增加, 参数增加, 训练难度增大, 导致通用性变差.

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 SRU的网络层数对评价指标的影响Fig.4 Effect of SRU network layers on evaluation metrics

不同嵌入维度对评价指标的影响如图5所示.由图可看出, 嵌入维度为10时AUC值低于嵌入维度为30的AUC值, 当嵌入维度为30时效果最佳, 而嵌入维度为50的AUC值、F₁值均有所下降.其原因可能是虽然嵌入维度越高能保留的特征信息也越多, 但引入的噪声也越多, 导致时间复杂度和空间复杂度也会随之升高.这在一定程度上表示CPIs预测问题适用于低维嵌入, 既能获取足够的特征信息, 也不会增加额外的时间和空间复杂度.

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 不同嵌入维度对评价指标的影响Fig.5 Effect of different embedding dimensions on evaluation metrics

通常随着网络层数的加深, 内部协变量偏移问题可能会愈加严重, 因此是否使用特征向量标准化的模型对AUC值和F₁值的影响如图6所示.由图可看出, 相比未进行标准化操作的GAT-GRU, 使用标准化操作的GAT-GRU的AUC值和F₁值更高, 这可能是由于输入数据经过网络处理后, 概率分布会发生偏移, 导致数据的表征能力下降, 影响网络收敛.但是, 对其进行标准化操作能使经网络处理后的数据分布与输入数据的分布基本保持一致, 有助于方法收敛.

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 特征向量标准化对评价指标的影响Fig.6 Effect of eigenvector normalization on evaluation metrics

GAT-GRU中不同模块数对AUC值和F₁值的影响如图7所示.由图可看出, 模块数量从3增加到4, F₁值增加, AUC值并无大幅提升, 当模块数增加到5时, AUC值和F₁值甚至略有下降.其原因可能是更深的网络带来的网络退化, 并且随着模块数的增加, 计算开销也会随之增加, 因此选择合适的模块数对实验结果至关重要.

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 不同模块数对评价指标的影响Fig.7 Effect of the number of different modules on evaluation metrics

为了进一步验证本文方法的泛化能力, 在2个数据集上进行评估实验.设置数据集上80%的样本为训练集, 20%的样本为测试集.对比方法采用性能较优的GCN-CNN^[16], 结果如表5所示.由表可看出, 本文方法略优于GCN-CNN, 泛化能力更优.从测试结果的方差可见, 本文方法的方差较小, 说明本文方法的性能较稳定.

表5

Table 5

表5(Table 5)

表5 两种方法的泛化性能对比 Table 5 Generalized performance comparison of 2 methods 方法 Human C.elegans ACC F1 ACC F1 GCN-CNN 0.973 ± 0.05 0.962 ± 0.06 0.974 ± 0.05 0.965 ± 0.07 本文方法 0.984 ± 0.03 0.983 ± 0.04 0.987 ± 0.03 0.986 ± 0.03

表5 两种方法的泛化性能对比 Table 5 Generalized performance comparison of 2 methods

在实际问题中, 药物发挥作用体现在小分子化合物与蛋白质活性位点发生交互, 使蛋白质的结构略有变化, 以此调节机体的生理活动.故不同位点相互作用强度有所不同, 通过使用注意机制可考虑到蛋白质中某些重要区域对交互的影响, 因此以得分较高的注意力权重对蛋白质的向量进行表示.

注意力机制对蛋白质向量表示的影响如图8所示, 由图可看出, 使用注意力机制的蛋白质向量表示对化合物与氨基酸相互作用越强, 位点颜色越深.通过减少对非重要位点的注意力, 突出显示蛋白质的重要区域, 达到增加模型准确率的效果.

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 是否使用注意力机制对蛋白质向量表示的影响Fig.8 Effect on protein vector representation using the attention mechanism or not