基于深度学习的冠状病毒刺突蛋白拉曼特征峰的理论研究

倪 爽， 温家星， 周民杰， 黄景林， 乐 玮， 陈 果，何智兵， 李 波， 赵松楠， 赵宗清， 杜 凯

中国工程物理研究院激光聚变研究中心， 四川 绵阳 621900

引 言

新冠疫情对全球的经济造成了巨大破坏。 为了有效控制新冠疫情， 快速检测新冠病毒(严重急性呼吸综合症冠状病毒2, severe acute respiratory syndrome coronavirus 2, SARS-CoV-2)是一个急需解决的问题。 新冠病毒的刺突蛋白是病毒攻击人体的钥匙， 且在病毒表面大量分布， 因此成为拉曼光谱技术检测新冠病毒的检测点[1-3]。 要实现刺突蛋白的拉曼技术检测， 关键之一在于构建刺突蛋白的拉曼特征峰。 鉴于纯刺突蛋白的拉曼谱图较难获得且成本较高， 需要从理论上快速构建刺突蛋白拉曼特征峰。 此外， 人类已知可以感染的七种冠状病毒刺突蛋白的结构[4-9]相近[图1(a)]， 是否可以通过拉曼技术区分他们对于新冠病毒的准确检测是一个十分重要的问题。 基于理论构建的刺突蛋白拉曼特征峰， 可以分析七种冠状病毒刺突蛋白拉曼特征峰的不同， 为实验提供指导。

(1)

式(1)中：εi为局域模频率，βij为局域模之间的相互作用， 激子模型哈密顿量可以写成矩阵形式

(2)

式(2)中：n为局域模的数量， 对角化H矩阵可以获得激子的频率及其对应特征矢量(局域模展开为简正模的线性展开系数)。 对于蛋白酰胺Ⅰ峰而言， 其局域模近似为组成蛋白骨架的每个酰胺单元的酰胺Ⅰ振动模。

本文关注的重点是定性分析七种冠状病毒刺突蛋白拉曼特征峰的差异而非绝对值且在激子模型中H矩阵对角元的值远大于非对角元的值。 因此， 本文采用简化的激子模型， 仅考虑结构变化对H矩阵对角元的校正而将非对角元舍去。

为了校正对角元， 需要对七种刺突蛋白结构中的每个酰胺单元计算拉曼谱图， 每个刺突蛋白有数千个酰胺单元， 逐个计算拉曼谱图费时且不利于统一误差。 江俊[15-16]等近期提出了基于机器学习构建分子结构和拉曼性质(频率、 强度)之间映射关系的策略。 其先通过分子动力学构建分子的动力学结构， 然后计算结构的拉曼频率以及强度， 最后使用机器学习拟合结构和性质的映射关系。 本文基于此策略， 构建酰胺单元结构和拉曼峰频率以及强度之间的模型。

文献中多数研究集中在酰胺Ⅰ峰， 对酰胺Ⅲ拉曼峰研究较少。 本文基于深度学习技术， 从理论上构建蛋白酰胺单元结构和酰胺Ⅰ、 Ⅲ拉曼特征峰的映射关系， 然后统计七种冠状病毒刺突蛋白的结构差异， 带入模型中获得拉曼特征峰。 最后通过洛伦兹线型展开谱线获得谱图， 并比较谱图的差异。

1 模型的构建和计算方法

使用N-甲基乙酰胺(NMA)分子来模拟蛋白的酰胺单元[图1(b， c)]， 构建酰胺Ⅰ、 Ⅲ峰的深度学习模型。 首先， 通过VASP软件， 采用分子动力学的方法构建10 000个NMA分子随时间演化的结构， 未考虑H2O的影响。 然后通过Gaussian09软件计算这些结构对应的拉曼谱峰， 最后通过深度学习技术构建NMA分子结构特征和酰胺Ⅰ、 Ⅲ峰频率以及强度的映射关系。 VASP计算采用GGA-PBE泛函， 周期性的边界条件为14.6×14.6×14.6 Å3， 平面波的截断能设置为400 eV， 能量收敛标准为1×10-5eV， 分子动力学模拟选择正则系综(NVT)， 时间步长1 fs， 温度为300 K， 总共模拟10 000步。 Gaussian09计算拉曼谱峰基于B3LYP杂化泛函水平， 基组使用6-311++G(d, p)。

深度学习模型由1个输入层， 3个隐藏层以及1个输出层组成， 对于每一个隐藏层， 使用线性修正单元(Rectified Linear Unit)激活函数。 对于NMA分子， 选用10个结构特征来描述(4个键长， 4个键角， 2个二面角， 如图2)。 为了增强模型的鲁棒性， 对输入特征以及输出结果进行了归一化处理。

图2 用来预测NMA分子拉曼位移以及强度的描述符

2 结果与讨论

2.1 酰胺Ⅰ、 Ⅲ峰模型

基于分子动力学构建了10 000个NMA分子结构， 计算了其均方根误差， 为1.06 Å(图3)， 这说明通过动力学演化得到的NMA分子结构变化较大且相关性较小。

图3 NMA分子结构的均方根误差

图4 NMA描述符间的皮尔逊相关系数

图5 酰胺Ⅰ峰的拉曼位移(a)、 强度 (b)以及酰胺Ⅲ峰的拉曼位移(c)、 强度(d)通过深度学习模型和DFT计算的结果比较

2.2 七种冠状病毒刺突蛋白骨架的结构特征

基于实验上的冠状病毒刺突蛋白结构[4-9](PDB code: 6u7h, 5i08, 5szs, 6ohw, 5x59, 5x58, 6vsb)， 统计其骨架特征键长、 键角、 二面角的分布。

图6 七种冠状病毒刺突蛋白骨架的键长特征统计分布图

从键角[图7(a—d)]的分布可以看出： 键角的分布较为广泛且均匀， 说明键角变化的力常数较小， 七种冠状病毒的键角变化相差不大， 键角对酰胺特征峰的影响较小。

图7 七种冠状病毒刺突蛋白骨架的键角和二面角特征统计分布图

从二面角的分布[图7(e,f)]中可以看出， 七种冠状病毒刺突蛋白的二面角均在180°(3.14弧度)附近， 这是由于酰胺平面共轭导致的。

2.3 七种冠状病毒刺突蛋白的拉曼特征峰比较

根据前面获得的酰胺Ⅰ、 Ⅲ峰的模型(酰胺单元的10个特征和酰胺Ⅰ、 Ⅲ振动峰的频率以及强度的映射关系)， 从七种冠状病毒刺突蛋白的实验结构中计算出每个酰胺单元的10个特征， 带入到模型中获得每个酰胺单元的酰胺Ⅰ、 Ⅲ拉曼峰(振动频率和强度)， 最后通过洛伦兹线型将每个酰胺单元的Ⅰ、 Ⅲ拉曼峰展开获得七种冠状病毒刺突蛋白的酰胺Ⅰ、 Ⅲ谱带(图8)。 从图中可以看出， 七种冠状病毒刺突蛋白的酰胺Ⅰ、 Ⅲ谱带根据最高峰频率可以各分成三个组。 对于酰胺Ⅰ峰， SARS-CoV-2， SARS-CoV， MERS-CoV刺突蛋白频率相近， 其频率值在1 636～1 637 cm-1区间； HCoV-HKU1， HCoV-NL63刺突蛋白频率相近， 其频率值在1 657～1 658 cm-1区间； HCoV-229E， HCoV-OC43刺突蛋白频率相近， 其频率值在1 673～1 674 cm-1区间。 对于酰胺Ⅲ峰， SARS-CoV-2， SARS-CoV， MERS-CoV刺突蛋白频率相近， 其频率值在1 263～1 265 cm-1； HCoV-229E， HCoV-HKU1， HCoV-NL63刺突蛋白频率相近， 其频率值在1 272～1 275 cm-1； HCoV-OC43刺突蛋白单独一个频率， 其频率值为1 285 cm-1。

图8 七种冠状病毒刺突蛋白的酰胺Ⅰ、 Ⅲ拉曼特征峰比较

根据上面的分析， 可以根据酰胺Ⅰ、 Ⅲ峰的频率划分七种冠状病毒， 如图9所示， 七种冠状病毒分为四组： SARS-CoV-2， SARS-CoV， MERS-CoV在同一个组； HCoV-HKU1， HCoV-NL63为一组； HCoV-229E一组； HCoV-OC43一组。 不同组之间其刺突蛋白特征峰的数值差异较大， 可以区分开来。 同组中的刺突蛋白特征峰的数值差异较小， 较难区分。

图9 通过酰胺Ⅰ、 Ⅲ特征峰频率区分七种冠状病毒

3 结 论

基于深度学习的技术构建了刺突蛋白酰胺Ⅰ、 Ⅲ特征拉曼峰模型， 结合实验上的冠状病毒刺突蛋白结构， 获得了七种冠状病毒刺突蛋白的酰胺Ⅰ、 Ⅲ拉曼特征峰， 通过比较七种冠状病毒刺突蛋白的拉曼特征峰， 可以把七种冠状病毒分为四组： SARS-CoV-2， SARS-CoV， MERS-CoV一组； HCoV-HKU1， HCoV-NL63一组； HCoV-229E一组； HCoV-OC43一组。 不同组的冠状病毒特征峰差异较大， 可以区分开来； 同一组的冠状病毒特征峰差异较小， 较难区分。