CN116072234B - 基于深度学习的分子动力学轨迹分析方法、存储介质和装置 - Google Patents

Abstract

本发明涉及基于深度学习的分子动力学轨迹分析方法、存储介质和装置，属于生物信息学领域，所述方法首先构建一个包含蛋白质结构以及化合物信息的初始数据集；然后将初始数据集进行分子动力学模拟，之后以20ps的间隔从每个200ns轨迹中提取构象快照，并转化为像素图作为后续机器学习分析的构象数据集；然后，提出一种新颖的分子动力学轨迹分析模型，其融合Swin‑Transformer和全连接神经网络，从而提高分子动力学轨迹分析的准确性；还融合了SHAP算法，能够实现对模型的有效解释并确定不同活性构象状态的重要残基。本发明还提供了执行所述方法的计算机可读存储介质和设备。

Description

基于深度学习的分子动力学轨迹分析方法、存储介质和装置

技术领域

本发明属于生物信息学领域，结合生物化学和分子动力学以及深度神经网络的知识设计并实现了一种基于深度学习的分子动力学轨迹分析方法。

背景技术

早期的研究认为，蛋白质的天然状态对应自由能最小的结构唯一状态。但是近年来越来越多的研究表明，蛋白质天然状态应该是自由能底部的一组构象的集合，蛋白质可以在这组构想中互相转变，这就是蛋白质的功能运动。生物体内蛋白质的功能运动往往需要蛋白质发生尺度较大的构象变化，这与蛋白质功能动力学有关，需要在原子层次上完整的理解蛋白的多个尺度的构象运动之间的耦合关系。简而言之，蛋白质构象变化在蛋白质功能中起重要作用。然而，涉及构象的原子细节对于实验技术来说通常是一项艰巨的任务，特别是对于复杂环境中的大型生物分子系统。

分子动力学(MD)模拟可以在原子水平上提供构象变化，从而为揭示许多生物系统的结构和功能机制做出了巨大贡献。但分子动力学模拟产生的蛋白质构象数量繁多且是高维的，随着系统的复杂度以及模拟时间的增加，这些构象结果数量也会随之激增。目前，手动分析仍存在忽略重要结构变化的风险，并且湿实验极其耗费人力物力。因此，如何智能的从广阔的构象空间中捕捉到功能特性的重要结构变化是目前领域内要面临的一个紧迫问题。

机器学习(ML)可以挖掘复杂数据的重要信息，因此在各个领域都取得了巨大的成功，包括分子动力学领域，其中机器学习可以用于分析蛋白质口袋动力学，增强采样，并生成新的构象。但样本的特征表示是ML模型的关键，ML模型一般通过计算一些特定的描述符来表示构象结构，这一方面会带来数据准备的压力，另一方面很难避免结构信息的丢失，因此必须探索一种简单又充分的构象表示方法。事实上，像素图表示的方法能够有效的解决这一问题，虽然近些年已经存在使用CNN或RNN对图像进行特征提取的例子，但融合Swin-Transformer的深度神经网络与构象像素表示能够更好的结合，并已经表现出了图像分类的高精度。这得益于Swin-Transformer在图像识别、像素级语义分割、图像分类等领域中强大的学习能力，以及通用性强、感受野巨大带来的比分子描述符更少的像素表示信息损失。然而，在使用深度学习方法时，我们必须面对另一个挑战，即其本质的黑盒。换句话说，在许多情况下，深度学习可以实现与人类相当甚至更好的预测准确性，但它的黑盒阻碍了我们对结果的理解，从而限制了其在生物分子系统中的应用。

事实上，近几年越来越多的研究者尝试将模型中加入模型解释算法，例如：LIME算法、shapley value算法等。然而，LIME算法只能通过线性模型的局部近似来进一步解释模型，仍然存在很大的误差，而shapley value算法的穷举以及权重分配极大的耗费运算资源。显然，要真正提高模型的解释准确性，就需要探索一种更加简单且充分的算法。SHAP算法融合了LIME算法的线性局部近似的思想以及shapley value算法的权重分配方法，能够更准确的反映出每一个样本中各特征的影响力的正负和大小。

综上，目前尚未有人尝试将构象像素表示、Swin-Transformer融合到深度神经网络中对分子动力学轨迹进行分析，以提升模型预测准确度，并使用SHAP算法对模型进行解释，提升模型可读性，从本质上发现科学规律。

发明内容

本发明要解决的技术问题在于提供一种基于深度学习的分子动力学轨迹分析算法。所述算法首先构建一个包含蛋白质结构以及化合物信息的初始数据集；然后将初始数据集进行分子动力学模拟，之后以20ps的间隔从每个200ns轨迹中提取构象快照，并转化为像素图作为后续机器学习分析的构象数据集；然后，提出一种新颖的分子动力学轨迹分析模型，其融合Swin-Transformer和全连接神经网络，从而提高分子动力学轨迹分析的准确性。除此之外，还融合了SHAP算法，以实现对模型的有效解释并确定不同活性构象状态的重要残基。

本发明是通过如下技术方案来实现的：

一种基于深度学习的分子动力学轨迹分析算法，所述算法的具体步骤如下所示：

步骤一、构建初始数据集；根据现有技术公开信息，从蛋白质数据库(PDB)获得所有正在研究的系统的初始晶体结构以及相关信息，设计并构建一个包含蛋白质结构以及化合物信息的初始数据集；

步骤二、对初始数据集进行分子动力学模拟；进行分子动力学模拟时，在晶体结构中保留了配体和换能器；模拟系统中，所有蛋白质残基在pH＝7时均设置为标准CHARMM质子化状态；然后，将这些晶体结构浸入4：1质量比的1-棕榈酰基-2-油酰基磷脂酰胆碱：胆固醇扩展膜中，并用TIP3P水模型溶剂化，且通过CHARMM-GUI平台设置0.15M的NaCl对其进行中和，然后进行5000步能量最小化；之后，每个系统在250ps内从0加热到310K，并在310K下通过5ns NVT(微正则系综)预平衡；最后进行了3次初始随机速度差异的并行200ns仿真，以增强构象采样；

步骤三、MD构象的像素表示；进行分子动力学模拟之后，从轨迹中提取快照，作为后续机器学习分析的构象数据集，其中仅通过从分子动力学快照中移除配体和换能器来考虑受体构象；使用像素贴图来表示数据集中的每个构象，其中每个像素点对应于一个原子；在转换为像素贴图之前，所有分子动力学坐标都已对齐以删除平移和旋转，且H原子没有被考虑在内；RGB色彩空间用于像素表示；

步骤四、融合Swin-Transformer和全连接神经网络的分子动力学轨迹分析新模型；模型融合Swin-Transformer，采取层次化设计；通过包含的4个隐藏层对特征图进行缩小分辨率操作，扩大感受野；

进一步，特征图依次通过patch partition层、linear Embedding层、SwinTransformer Block层以及三个Patch Merging和Swin Transformer Block复合层；为增强模型学习准确率并增加计算效率，Patch Merging层通过降采样以及全连接层改变特征图维度，并且特征图在第二个复合层中需要循环通过三个Swin Transformer Block层，使模型在特征图为的维度上进一步学习；每一个Swin Transformer Block层都包含两个子层，两个子层依次运行，且每个子层都是由Layer Normalization层、WindowAttention层、残差连接以及一个Layer Normalization和MLP复合层组成的；唯一区别在于第一个子层的Window Attention层为W-MSA层，第二个子层为Shifted Window Attention层SW-MSA，所以每一个Swin Transformer Block层都经过两次self-attention计算，以此来获取全局的视野；通过Swin-Transformer后构建了一个三层全连接的神经网络，最后输出神经元数为两个，并使用激活函数softmax分类；

步骤五、使用步骤三中处理后的数据对步骤四中所述的新模型进行训练和评估；使用五重交叉验证，将每条轨迹按照时间顺序分组，每组又有一个以上的折叠；然后，从每组中取出一个折叠来构成验证集，而将每组中剩余的其余折叠组合起来作为训练集；依次重复五次，构建五重交叉验证集；并使用Accuracy评估新模型的性能；

步骤六、可解释模型的SHAP算法

对特征图进行分割，并为分割区域进行编号，通过对特征图部分分割区域的遮掩进行扰动，将扰动后的所有特征图重新进入模型进行分析预测，得到带有小扰动的扰动数据集，构建样本空间，最后通过局部线性近似拟合一次函数，每一个未知数对应分割区域，未知数的系数为起作用大小；SHAP算法加入扰动权重，权重大小通过保留未遮掩的区域多少决定。

本发明还提供一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序适用于由处理器加载并执行所述分子动力学轨迹分析算法。

本发明还提供一种计算机设备，所述设备包括存储器和处理器，所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器执行所述分子动力学轨迹分析算法步骤。

本发明与现有技术相比的有益效果：

1、使用像素图的方式表示分子动力学轨迹，构建像素表示的构象数据集。数据集摈弃现有轨迹表示方法，通过使用像素图的形式表示了分子动力学的轨迹，减少了蛋白质构象的结构信息丢失，大大提升了后续新模型的学习准确率。

2、设计并实现了融合Swin-Transformer和全连接神经网络的分子动力学轨迹分析新模型。模型采用图像分类中相对新颖的Swin-Transformer方法对表示分子动力学轨迹的像素图进行学习。其W-MSA层和SW-MSA层通过self-attention算法使模型能够有全局的视野进行学习，增强了预测的准确性。并且通过Patch Merging的降采样以及全连接神经网络对特征图进行降维，使模型能够在不同的维度对特征图进行学习，在合适的维度进行深入学习，使模型训练时间短，可用性强。与传统的CNN图卷积相比有更强的学习能力以及更高的预测准确度。

3、引入并实现了可解释模型的SHAP算法；算法通过遮掩扰动特征图获取扰动数据集，并对数据集再次通过模型进行分析预测得到样本空间。最后通过加入扰动权重和局部线性近似拟合一次函数得到的分类结果，快速识别不同活动状态下的重要残基，增强模型的可解释性。

附图说明

图1为本发明的整体框架图；

图2为本发明MD构象与像素图转化原理示意图；

图3为SHAP解释器流程图。

具体实施方式

下面通过实施例结合附图来对本发明的技术方案做进一步解释，但本发明的保护范围不受实施例任何形式上限制。

实施例1

一种基于深度学习的分子动力学轨迹分析算法，如图1所示，所述算法的具体步骤如下所示：

步骤一、构建初始数据集

从蛋白质数据库(PDB)获得所有正在研究的系统的初始晶体结构以及相关信息，构建包含蛋白质结构以及化合物信息的初始数据集。

步骤二、对初始数据集进行分子动力学模拟；

在模拟中，将配体和换能器保留在晶体结构中。在分子动力学系统中，所有蛋白质残基在pH＝7时设置为标准CHARMM质子化状态。然后，将这些蛋白质浸入4：1质量比的1-棕榈酰基-2-油酰基磷脂酰胆碱：胆固醇扩展膜中，并用TIP3P水模型进行溶剂化，通过CHARMM-GUI在0.15M NaCl中中和，然后进行5000步能量最小化。之后，每个系统在250ps内从0加热到310K，并在310K下通过5ns NVT进行预平衡。为了增强构象采样，在初始随机速度不同的三个平行200ns模拟中进行每个系统的温度为310K，压力为1个大气压。对于所有模拟，采用2fs步长，并使用SHAKE算法来约束所有涉及氢的化学键。使用粒子网格Ewald方法计算静电相互作用。CHARMM36力场用于蛋白质、脂质和盐离子，而所有配体的参数都是使用CHARMM一般力场生成。

步骤三、MD构象的像素表示

在MD之后，以20ps的间隔从每个200ns轨迹中提取10000个快照，作为后续ML分析的构象数据集，通过从MD快照中去除配体和换能器仅考虑受体构象。为了尽量减少结构信息损失，使用像素图来表示数据集中的每个构象，其中每个像素点对应一个原子，转化原理如图2所示。需要注意的是，所有MD坐标都对齐以在转换为像素图之前移除平移和旋转且不考虑H原子。最后通过矩阵变换将每个构象的XYZ坐标转换为RGB坐标，得到特征像素图数据集。具体转化原理如图2所示。RGB色彩空间用于像素表示，因为它是通过更改和叠加三个颜色通道来表示颜色的标准：红色(R)、绿色(G)和蓝色(B)。RGB空间中的每个点都对应一种独特的颜色；特别是，每个构象的XYZ坐标通过矩阵变换传输到RGB坐标；

步骤四、融合Swin-Transformer和全连接神经网络的分子动力学轨迹分析新模型；

模型融合Swin-Transformer，采取层次化设计；通过包含的4个隐藏层对特征图进行缩小分辨率操作，扩大感受野；

特征图依次通过patch partition层、linear Embedding层、Swin TransformerBlock层以及三个Patch Merging和Swin Transformer Block复合层；为增强模型学习准确率并增加计算效率，Patch Merging层通过降采样以及全连接层改变特征图维度，并且特征图在第二个复合层中需要循环通过三个Swin Transformer Block层，使模型在特征图为的维度上进一步学习；每一个Swin Transformer Block层都包含两个子层，两个子层依次运行，且每个子层都是由Layer Normalization层、Window Attention层、残差连接以及一个Layer Normalization和MLP复合层组成的；唯一区别在于第一个子层的Window Attention层为W-MSA层，第二个子层为Shifted Window Attention层SW-MSA，所以每一个Swin Transformer Block层都经过两次self-attention计算，以此来获取全局的视野；通过Swin-Transformer后构建了一个三层全连接的神经网络，最后输出神经元数为两个，并使用激活函数softmax分类。

步骤五、使用步骤3)中处理后的数据对步骤4)中所述的新模型进行训练和评估；使用五重交叉验证，将每条轨迹按照时间顺序分为10组，每组又分为5个折叠；然后，从每组中取出一个折叠来构成验证集，而将每组中剩余的四个折叠组合起来作为训练集；依次重复五次，构建五重交叉验证集；并使用Accuracy评估新模型的性能；

对于分子动力学轨迹分析模型，其包含的4个隐藏层对特征图进行缩小分辨率操作，扩大感受野。输入一个H×W×3的图片，首先将图片达成patch，patch size设置成4×4。然后经过patch partition后将图片尺寸变为之后通过linear Embedding层以及Swin Transformer Block层，其中linear embedding层会设置超参数C＝96。之后在特征图在尺寸为/>和

时通过一个Swin Transformer Block层，在尺寸为/>时通过三个Swin TransformerBlock层。最后通过三层全连接的神经网络分别包含512、64和2个神经元，以及0.5的dropout防止过拟合。对于512和64个神经元的两个全连接层，使用ReLU作为激活函数，而在最后一层使用softmax激活函数进行分类。

为了验证模型的鲁棒性，在一条平行MD轨迹上使用训练好的模型来预测其他两条平行轨迹。结果如表1所示。验证集的预测准确度达到100％，独立测试集的预测准确度超过99％，显示出高鲁棒性。

表1.模型在三种类型的G蛋白偶联受体独立测试集上的表现

步骤六、基于SHAP的解释器

对特征图进行分割，并为分割区域进行编号，通过对特征图部分分割区域的遮掩进行扰动，将扰动后的所有特征图重新进入模型进行分析预测，得到带有小扰动的扰动数据集，构建样本空间，最后通过局部线性近似拟合一次函数，每一个未知数对应分割区域，未知数的系数为起作用大小；SHAP算法加入扰动权重，权重大小通过保留未遮掩的区域多少决定，如图3所示。

对于选择的特定构象，首先添加一些小扰动以获得包含与原始构象相似的构象的扰动数据集。πx用于定义所选构象x周围的局部性。使用训练的新模型在扰动数据集上进行预测。由于扰动数据集非常小且与所选构象相似，分类边界变得简单，接近线性。因此，我们可以根据扰动大小确定权重，并使用直观且可解释的线性模型除以权重来拟合局部边界。具体而言，我们在扰动数据集上训练线性模型，使新模型的预测值与线性模型的预测值之间的误差尽可能小。对于每个构象，SHAP解释器生成一个SHAP矩阵，该矩阵评估每个像素在分类结果中的重要性得分。SHAP矩阵与图像具有相同的维度，其中每个元素对应一个代表原子的像素点。SHAP矩阵中每个元素的取值为0或1。0值表示该元素对分类选择的影响很小，而1值表示该元素对分类决策的影响很大。我们将每个构象的所有SHAP矩阵求和，并对它们进行平均以获得0到1值之间的分数。值越大，分类结果中的原子越重要。然后，对残基中所有原子的得分进行平均，以表示残基的重要分数。

Claims (3)

1.一种基于深度学习的分子动力学轨迹分析方法，其特征在于，所述方法的具体步骤如下所示：

步骤一、构建初始数据集；从蛋白质数据库获得所有正在研究的系统的初始晶体结构以及相关信息，设计并构建一个包含蛋白质结构以及化合物信息的初始数据集；

步骤二、对初始数据集进行分子动力学模拟；进行分子动力学模拟时，在晶体结构中保留了配体和换能器；模拟系统中，所有蛋白质残基在pH＝7时均设置为标准CHARMM质子化状态；然后，将这些晶体结构浸入4：1质量比的1-棕榈酰基-2-油酰基磷脂酰胆碱：胆固醇扩展膜中，并用TIP3P水模型溶剂化，且通过CHARMM-GUI平台设置0.15M的NaCl对其进行中和，然后进行5000步能量最小化；之后，每个系统在250ps内从0加热到310K，并在310K下通过5ns微正则系综预平衡；最后进行了3次初始随机速度差异的并行200ns仿真，以增强构象采样；

步骤三、MD构象的像素表示；进行分子动力学模拟之后，从轨迹中提取快照，作为后续机器学习分析的构象数据集，其中仅通过从分子动力学快照中移除配体和换能器来考虑受体构象；使用像素贴图来表示数据集中的每个构象，其中每个像素点对应于一个原子；在转换为像素贴图之前，所有分子动力学坐标都已对齐以删除平移和旋转，且H原子没有被考虑在内；RGB色彩空间用于像素表示；

步骤四、融合Swin-Transformer和全连接神经网络的分子动力学轨迹分析新模型；模型融合Swin-Transformer，采取层次化设计；通过包含的4个隐藏层对特征图进行缩小分辨率操作，扩大感受野；

步骤五、使用步骤三中处理后的数据对步骤四中所述的新模型进行训练和评估；使用五重交叉验证，将每条轨迹按照时间顺序分组，每组又有一个以上的折叠；然后，从每组中取出一个折叠来构成验证集，而将每组中剩余的其余折叠组合起来作为训练集；依次重复五次，构建五重交叉验证集；并使用Accuracy评估新模型的性能；

步骤六、可解释模型的SHAP算法；对特征图进行分割，并为分割区域进行编号，通过对特征图部分分割区域的遮掩进行扰动，将扰动后的所有特征图重新进入模型进行分析预测，得到带有小扰动的扰动数据集，构建样本空间，最后通过局部线性近似拟合一次函数，每一个未知数对应分割区域，未知数的系数为起作用大小；SHAP算法加入扰动权重，权重大小通过保留未遮掩的区域多少决定；

所述的步骤四的具体方法为特征图依次通过patch partition层、linear Embedding层、Swin Transformer Block层以及三个Patch Merging和Swin Transformer Block复合层；为增强模型学习准确率并增加计算效率，Patch Merging层通过降采样以及全连接层改变特征图维度，并且特征图在第二个复合层中需要循环通过三个Swin Transformer Block层，使模型在特征图为的维度上进一步学习；每一个Swin Transformer Block层都包含两个子层，两个子层依次运行，且每个子层都是由LayerNormalization层、WindowAttention层、残差连接以及一个LayerNormalization和MLP复合层组成的；唯一区别在于第一个子层的Window Attention层为W-MSA层，第二个子层为ShiftedWindowAttention层SW-MSA，所以每一个Swin TransformerBlock层都经过两次self-attention计算，以此来获取全局的视野；通过Swin-Transformer后构建了一个三层全连接的神经网络，最后输出神经元数为两个，并使用激活函数softmax分类。

2.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有计算机程序，所述计算机程序适用于由处理器加载并执行权利要求1所述分子动力学轨迹分析方法。

3.一种计算机设备，其特征在于，所述设备包括存储器和处理器，所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器执行权利要求1所述分子动力学轨迹分析方法步骤。