基于深度学习的SARS-CoV-2 RBD-ACE2结合亲和力预测方法

吴敏明，李维华**，王玉霜，丁海燕

(1.云南大学 信息学院，云南 昆明 650500；2.云南大学 云南生物资源保护和利用国家重点实验室，云南 昆明 650091)

新型冠状病毒感染疫情自2019 年底暴发以来，迅速波及全球多个国家和地区，引起全球大流行[1].截至2022 年4 月13 日，全球已向世界卫生组织(World Health Organization，WHO)报告了497 960 492 例确诊病例，其中6 181 850 例死亡，死亡率达1.24%[2].

SARS-CoV-2 属于巢病毒目（Nidovirales）冠状病毒科（Coronaviridae）的β属[3]，是一种具有单股正链RNA 基因组的包膜病毒.其Spike 蛋白（S）是一种Ⅰ型融合蛋白，在病毒粒子表面形成三聚体.SARS-CoV-2 在进入靶细胞(Target cell)时，弗林蛋白酶（Furin/TMPRSS2）会把S 蛋白水解成两个功能亚基[4]：S1 和S2 亚基.S1 亚基上的受体结合域（Receptor Binding Domain,RBD）主要负责与人体细胞表面受体—血管紧张素转化酶2（Angiotensin Converting Enzyme，ACE2）结合[5].因此，S 蛋白在一定程度上决定了病毒的传染性及其在宿主中的传播性[6].在RBD 与ACE2 结合后，S1 和S2亚基发生一系列构象变化，并介导病毒和宿主细胞膜融合，使病毒基因组进入细胞，具体过程如图1所示.

图1 新型冠状病毒融合过程示意Fig.1 Schematic diagram of the fusion process of SARS-CoV-2

对SARS-CoV 的流行病学和生物化学研究表明，不同SARS-CoV 在宿主细胞中的感染性与每种毒株的RBD 与宿主细胞的ACE2 之间的结合亲和力成正比[7].因此，获取RBD-ACE2 结合亲和力对评估SARS-CoV-2 毒株的侵染能力、突变后的传染性变化以及未来感染趋势都具有重要意义.

结合亲和力，也称为结合自由能，可以通过生物学试验或者自由能计算模型获得.测定RBD-ACE2结合亲和力的生物学试验方法包括免疫荧光共振（Surface Plasmon Resonance）[8]、酶联免疫吸附测定（Enzyme-Linked Immuno Sorbent Assay，ELISA）[9]和定量深度突变扫描（Quantitative Deep Mutational Scanning Approach）[10]等.然而，生物试验方法费时费力，并且测定一般都在病毒流行一段时间后.生物试验方法还无法满足对突变株快速预测和监控的需求.自由能计算模型一般采用基于蛋白质物理化学性质的分子动力学进行模拟，例如MM/PB(GB)SA[11]方法有助于在分子层面上认识RBD-ACE2的结合机理，但突变株晶体结构很难获得，同时可信的结合自由能往往需要较长的模拟时间.实验数据的积累为机器学习方法直接模拟突变与结合亲和力变化之间的内在关系提供了前所未有的机会.例如Geng 等[12]利用界面结构、进化和近似能量得出的预测特征，基于随机森林预测突变后的亲和力变化.但是，这些方法所依赖的特征往往需要大量的模拟时间获得[13].

近年来，神经网络在预测蛋白质结合亲和力上取得了显著的成果.例如，Li 等[14]利用图神经网络学习原子层蛋白质结构特征，并利用这些特征评估蛋白质的结合亲和力，但是这个方法依赖于蛋白质原子层面的复杂特征；Chen 等[15]利用卷积神经网络抽取特征，再利用梯度提升树评估SARS-CoV-2的S 蛋白和宿主ACE2 受体在突变后的结合亲和力变化，但是这个方法依赖蛋白质的物理化学特征.测序技术可以低成本且快速地获得毒株的RBD序列信息，这为探索基于序列的结合亲和力计算方法提供了基础.

因此，本文结合卷积神经网络、循环神经网络和注意力机制，设计了一个基于序列的RBD-ACE2结合亲和力的计算模型.该模型首先对RBD 序列进行特征表示，其次通过深度神经网络提取关键的序列特征，并用这些特征预测RBD-ACE2 的结合亲和力.本文主要贡献如下：

（1）本文模型仅仅基于序列，预测多点突变毒株的RBD-ACE2 结合亲和力，避免复杂特征工程以及长时间模拟；

（2）在真实数据集上的实验结果表明，模型在预测RBD 多点突变毒株的RBD-ACE2 结合亲和力是有效的.

1 本文模型框架

基于深度神经网络的RBD-ACE2 结合亲和力预测模型的主要框架如图2 所示.RBD-ACE2 结合亲和力预测的任务是根据突变前后的RBD 序列，预测RBD-ACE2 结合亲和力的变化.该模型主要包括输入模块、卷积模块、循环模块、注意力模块和预测模块.输入模块利用ProtVec 方法[16]对RBD序列进行特征表示；卷积模块用于提取RBD 序列的局部信息；循环模块进一步提取序列的全局信息；注意力模块用于区分RBD 特征矩阵的重要性；预测模块用于输出RBD-ACE2 的结合亲和力.

图2 新型冠状病毒RBD-ACE2 结合亲和力的预测模型Fig.2 Prediction model of SARS-CoV-2 RBD-ACE2 binding affinity

1.1 输入模块近年来，分布式表示在自然语言处理(Natural Language Processing，NLP)训练中的单词嵌入、单词到数字向量空间的映射方面取得了巨大成功，已在生物信息学中有广泛的应用，如蛋白质的分类和结构的预测.本文对输入长度为n的RBD 序列按照k=3 且步长s=1 的k-mer[17]将RBD 序列划分为n-2 个氨基酸三元组；对突变前后的RBD 序列分别采用ProtVec[17]方法将其映射为特征矩阵Xi和Xj.进一步通过按位差运算获得毒株对的特征表示：

1.2 卷积模块卷积模块采用多尺度卷积神经网络提取上的多个局部特征，并且包括两个卷积层.其中，第一层卷积使用卷积核大小分别为f(f=3,5,7)的卷积神经网络提取上的局部特征，每个卷积核按照下面的公式进行映射：

式中：σ是激活函数ReLU，Wc是 权重矩阵，bc是偏置项，⊙是矩阵点乘操作.将卷积核得到的特征Cf进行最大池化操作=max{Cf}，分别将池化操作得到的特征Cf进行拼接，形成第一个卷积层的特征表示：

式中：||表示特征矩阵的连接.第二个卷积层接收第一个卷积层的输出Ct，并使用两个卷积核大小为f=3 的滤波器分别进行卷积操作：

1.3 循环模块长短期记忆网络(Long Short-Term Memory，LSTM)、门控循环单元(Gated Recurrent Unit，GRU)可以克服传统的递归神经网络中梯度消失和梯度爆炸问题.和LSTM 相比，GRU 可以用较少的参数获得较好的性能，因此，在循环神经模块中使用双向GRU(Bidirectional GRU,BGRU)更新模型状态并计算每个氨基酸三元组的全局特征.具体来说，将卷积模块学习到的局部特征向量表示为对于当前的输入，隐状态ht通过从正向获得的特征和反向获得的特征计算得到，计算公式如下：

式中：ht-1表示上一个时刻的隐状态.

对于BGRU 提取的RBD 序列全局特征表示为h=[h1,h2,···,ht].为了获得RBD 序列的整体表示，将卷积模块获取的局部特征和循环模块获取的全局特征拼接起来：

式中：Uf是由卷积模块提取的局部特征，h是由循环模型提取的全局特征.

1.4 注意力模块注意力模块输入RBD序列的整体表示g=由于不是所有的特征都对预测任务有同等的贡献，所以在训练中采用注意机制[18]自动学习和分配RBD 序列特征的权重，计算公式如下：

式中：Wc是 权重矩阵，bc是偏置项.注意力模块在RBD 序列特征的时间步长t处为每个k-mer 分配一个权重，并形成一个时间步长t处新的特征向量st，计算公式如下：

通过注意力模块得到的特征矩阵表示为s=[s1,s2,···,sn].

1.5 预测模块预测模块由全连接层组成，接收注意力模块学习的特征s，并使用非线性激活函数SeLU 计算预测值：

式中：γ 和 α是归一化后得到指定均值和方差.

模型参数通过训练样本上最小均方误差获得，预测模块进一步根据训练好的模型进行预测.均方误差的定义如下：

式中：yi和分别表示样本真实值和预测值，N表示训练样本.

2 实验设置

2.1 实验环境和数据本文实验环境的主要参数为：处理器Intel(R) Core(TM) i5-8250U CPU@1.60 GHz，内存8 GB，操作系统Windows 10；采用深度学习框架TensorFlow 2.2 以及内置的Keras 构建神经网络并进行模型的训练和测试.

这就好比要搞清楚自己接手的是什么样的工程，自己又该有哪些具体的实施构想，然后才能有的放矢对着图纸付诸施工，至于“梦想照进现实”，那只是时间的问题。

本文使用的RBD-ACE2 结合亲和力试验数据集[19]，包含4 003 条RBD 突变体氨基酸序列及其对应的结合亲和力.其中，3 388 条突变体的结合亲和力小于野生型毒株；615 条突变体的结合亲和力大于野生型毒株.具体统计数据信息如图3 所示.其中横坐标表示结合亲和力的分布区间，纵坐标表示各区间突变体的数量.

图3 RBD 突变体结合亲和力的区间分布Fig.3 Interval distribution of binding affinity of RBD mutants

2.2 模型参数设置卷积模块共有两个卷积层，第一层包括3 种类型的卷积运算，卷积核长度分别为3、5、7，每个卷积层之后的Max-pooling 运算，pooling 长度为2；第二层的卷积核长度为3.循环模块使用双向GRU，其中每个BGRU 的输出维度设置为128，共有256 个隐藏单元.预测模块由全连接层组成，全连接层的神经单元分别为128、64、1.

本文按照0.8∶0.2 的比例将数据集划分为训练集和测试集，使用随机梯度下降Adam 算法作为模型的优化器；初始学习率设为0.000 1；卷积模块和循环模块使用Dropout 策略和正则化参数L2 来避免训练过程中的过拟合问题，且Dropout 和L2分别设为0.5 和0.001.

2.3 评价标准为了全面地评估模型的性能，本文采用4 种不同的评价指标评估模型性能.

均方误差(mean squared error，EMS)：

平均绝对误差(mean absolute error，EMA)：

解释回归模型的方差(explained variance score，EVS)：

3 实验和结果分析

3.1 不同基线方法对比为了评估本文模型的性能，本文使用几种回归算法作为基线，包括SVR[20]、CART[21]、KNN[22]、Random Forest[23]、XGBoost[24]和MLP Regressor[25].表1 给出本文模型与这些基线方法在相同测试集上的RBD-ACE2 结合亲和力的预测性能.

表1 不同基线方法性能对比Tab.1 Performance comparison of different baseline methods

从表1 可以看出，在相同数据集的情况下，本文设计的模型获得的预测结果明显优于基线模型.例如，本文模型的EMS、EMA评价指标分别为0.53、0.42，分别比基线模型XGBoost 的EMS(0.89)和KNN的EMS( 0.69)低0.36 和0.17；本文模型的EVS、R2评价指标分别为0.76、0.76，分别比基线模型XGBoost的EVS( 0.57)和R2(0.57)高0.19.因此，实验结果验证了本文模型在RBD-ACE2 结合亲和力的预测上的有效性.

表2 嵌入向量结果对比Tab.2 Comparison of embedding vector results

从表2 可以看到，模型在预测RBD-ACE2 结合亲和力上的实验结果表明ProtVec 方法优于Word2Vec 和One-hot.因此，本文使用ProtVec 方法对RBD 序列进行特征进行表示是合理的.

3.3 不同循环模块的性能为了分析不同循环模块和不同损失函数值对预测RBD-ACE2 结合亲和力的影响，使用LSTM 模块替换本文模型中的GRU模块，并使用两种不同的损失函数MSE 和MAE对模型分别进行训练.得到RBD-ACE2 结合亲和力的预测效果如表3 所示.

表3 不同循环模块的性能对比Tab.3 Performance comparison of different cycle modules

从表3 中可以看出，当损失函数为MSE，循环模块采用GRU 时，本文模型的评价指标比使用LSTM 模块的效果更优；当损失函数为MAE，循环模块采用GRU 时，本文模型的评价指标同样比使用LSTM 模块的效果更优.因此，本文的模型使用GRU 作为循环单元模块，而且采用MSE 作为损失函数值评价模型.

3.4 预测模型对RBD 多点突变的预测性能为了验证本文模型在具有多点突变的毒株RBD-ACE2结合亲和力上的预测效果.由于沙贝病毒RBDs 最大似然发育树中相近的病毒结构相近[19].因此，本文以具有多点突变的SARS-CoV-2(沙贝病毒亚属)的不同亲缘关系毒株RaTG13、GD-Pangolin 和SARS-CoV-1 为例[19],作为对比毒株.如图4 所示.

图4 沙贝病毒 RBDs 的最大似然系统发育树Fig.4 Maximum likelihood phylogenetic tree of sarbecovirus RBDs

通过比较最大似然发育树上的不同亲缘毒株的模型预测值、加权平均法[19]和生物试验计算值，结果如表4 所示.

表4 不同亲缘关系毒株结合亲和力预测Tab.4 Predicted binding affinity of different strains with different affinities

从表4 可以看出，使用本文模型预测的结合亲和力比加权平均法[19]更接近实验值.对RaTG13，本文模型的预测值（-2.23）相比较加权平均（-4.31）更接近实验值（-2）；对GD-Pangoin，本文模型和加权平均法的预测值都和实验值非常接近.然而，对SARS-CoV-1，本文模型和加权平均法的预测值都离实验值较远，但本文模型的预测值更接近实验值.其原因可能是SARS-CoV-1 与SARS-CoV-2 属于沙贝病毒亚属，本文模型可以学习到部分RBD序列特征，从而使得本文模型的预测效果更接近试验值.

Starr 等[19]认为，当不同亲缘毒株的RBD 区域氨基酸和SARS-CoV-2 的差别不大时，其结构变化也不大，那么少量突变是在相同的背景下发生，多点突变产生的亲和力差别可以用单点突变的加权平均获得.当不同亲缘毒株的RBD 区域氨基酸和SARS-CoV-2 的差别较大时，其结构可能发生较大变化，多点突变产生的亲和力差别与单点突变的加权平均不相符.从上面的分析可以看出，本文模型可以成功地预测系统发育树上较近(相应地RBD突变较少)毒株的结合亲和力，但对预测系统发育树上较远毒株的结合亲和力有一定局限性.从图4可以看出，所有SARS-CoV-2 突变株在发育树上都离野生型很近，RBD 的突变也比较少.因此，我们模型具有预测RBD 多点突变毒株的RBD-ACE2结合亲和力的能力.

3.5 关切变异株中的应用分析WHO 公布了5种关切变异株Alpha、Beta、Gamma、Delta、Omicron.其中，Alpha 毒株与最初在武汉发现的野生株相比，其RBD 与ACE2 受体的结合亲和力提高了1.98 倍[26]，原因是N501Y 突变诱导中性电荷和官能团的疏水性增加，通过疏水相互作用稳定了与ACE2 受体结合；Beta 毒株与ACE2 的结合亲和力是Alpha 株的2.7 倍[27]；Gamma 毒株与早期流行毒株相比，传播能力增强约40%[28]，原因是该毒株的K417T 突变将S 蛋白中的赖氨酸转化为苏氨酸，从而形成中性电荷、氢键，并减少空间位阻.同时，K417T 和H655Y 的协同作用也可能增强病毒的传染性.Delta 株的传染性比Alpha 株高64%[26]，提高的原因可能是L452R、T478K 和E484Q 突变都改变了S 蛋白的表面电荷，病毒与带负电的细胞膜相互作用，增强细胞对病毒的摄取；Omicron 毒株的传染性可能是野生株的10 倍以上，该病毒通过突变增强RBD-ACE2 结合亲和力或通过突变逃避抗体保护来增强其在RBD 上的进化优势[29].

为了验证模型在新型冠状病毒关切变异株上的预测有效性，本文分别对这5 种关切变异株的RBD-ACE2 结合亲和力进行预测，将模型预测的值使用相对解离常数Δlg(KD,app)[19]表示.实验结果如图5 所示.

图5 关切变异株预测结果Fig.5 Predicted results of VOC

从图5 可以看到，本文模型对Alpha、Beta、Gamma、Delta 和Omicron 进行预测的相对解离常数Δlg(KD,app)分别为0.10、-0.12、0.06、0.17 和0.19.相对于原始SARS-CoV-2 毒株来说，除了Beta，其他变异株的预测值和现有文献公开结果一致[28-29]，结合亲和力都高于原始SARS-CoV-2 毒株.Beta 预测值和现有文献公开结果有出入.这一点可能和Starr 等[19]通过深度突变扫描获取的突变体K417N、E484K 和N501Y 的结合亲和力相关，因为K417N、E484K 和N501Y 分别是-0.45、0.06、0.24，K417N 上的扫描值支配了模型的预测值.同时，造成现实中Beta 株传播力高于野生毒株的原因可能是复杂的，例如，Beta 株Sipke 三聚体更容易采用站立态，利于和ACE2 结合[30]，或是K417N突变消除了RBD 和中和抗体CB6 之间埋藏的界面盐桥，降低了结合能，从而促进了该变体能有效逃逸[31].

从上面的分析可以看出，本文模型在具有RBD多点突变的关切变异株的RBD-ACE2 结合亲和力上具有良好的预测性能.

4 结束语

为了获取RBD-ACE2 结合亲和力，了解SARSCoV-2 在现有突变后的传染性变化和侵染能力，本文利用基于深度神经网络，设计了一个基于RBD序列的RBD-ACE2 结合亲和力的预测模型.该模型充分利用RBD 序列并进行特征表示，结合卷积神经网络、循环神经网络和注意力机制挖掘序列中影响亲和力的关键特征，并用于预测.在真实数据集对模型进行训练和评估，实验结果表明，5 种关切变异株中有4 个预测值和现有文献公开值一致.和现有的回归模型相比，本文模型能更有效地预测多点RBD 突变变异株的RBD-ACE2 结合亲和力.虽然本文模型具有良好的预测性能，但是本文仅基于序列对RBD-ACE2 结合亲和力进行预测，未考虑到其他辅助特征，未来可以考虑引入相关的蛋白质结构和理化性质以提高模型预测效果.