人工智能技术已经被应用于化学、生物、材料发现等专业领域[1-3]。在使用传统的机器学习方法预测物质的性质时，通常会使用一些已知的性质来预测某些性质，这些特征通常是人工选择的[1,4,5]。后来，由于物质的性质与它的分子结构密切相关，更多的人选择利用图神经网络根据分子图或分子空间信息来预测它的性质，分子空间信息有助于提高预测结果[6-11]。
目前，获取分子几何结构的传统方法通常是通过实验方法测定或通过分子模拟计算。实验方法包括微波光谱、X射线衍射、电子衍射和中子衍射。但实验方法复杂，中子衍射等实验设备成本较高。传统的分子模拟计算方法主要是利用密度泛函理论(DFT)进行计算和优化，这些计算方法计算步骤多，计算成本高，往往以小时为单位。此外，虽然一些经验计算方法具有较高的计算效率，但计算结果往往较差。
近年来，一些研究通过预测距离矩阵来重构原子坐标，得到分子的几何结构[12-17]。这些研究大多采用概率生成模型，目的是生成多个低能量的构象，但仍存在一定的局限性。例如，文献[15]只能处理由固定原子组成的分子，而文献[12,14]需要先通过一个构象得到原子间的距离。
针对这种情况，本文提出了一种基于图卷积网络[18]的分子距离矩阵预测模型(DMGCN)。该模型以简化的分子输入线输入系统(SMILES)转换成的图形表示作为输入，预测分子中原子之间的成对距离，用于分子几何结构的重建或分子性质的预测等下游任务。本文的主要贡献如下：(1)提出了一种新的用于分子距离矩阵预测的图卷积网络模型；(2)进行了对比实验探索和分析最大原子数目不同数据集对我们的模型的影响，结果表明，我们的模型比DeeperGCN-DAGNN[17]和RDKit[21]有着更小的MAE，并且我们的模型预测的键长是接近DFT计算和实验测量的值的；(3)利用预测距离矩阵对分子性质进行预测，验证了该方法的有效性。
百度刚刚发表的ChemRL-GEM网络模型中，3D空间结构是用RDKit获取的，本文模型的精度高于RDKit。
2 方法
2.1 编码
本文将分子定义为一个图G = (V, E)，其中图节点V表示原子，边E表示键。原子的特征向量定义为v, v∈V，键的特征向量定义为e, e∈E，其中v被编码为19维的特征向量，e被编码为5维的特征向量。具体特征向量构成如表1所示。
表1 特征向量组成
2.2 模型概述
本文提出的DMGCN模型以SMILES为输入，在预处理阶段将其转换为图，模型整体框图如图1所示。由于需要预测分子距离矩阵，所以在构建输入图时，需要将它们之间没有化学键的原子连接起来，即构建一个完整的图。而完全图破坏了分子图的原始结构，因此，模型建立了两个分支，分别对完全图和分子图进行处理。分子图所在的分支需要完成两件事：一是预测原子之间有键时的距离，二是通过图卷积运算[18]将分子图送入更新层进行节点更新。更新完成后，将所有节点的特征向量分配给完整图中的节点，完全图使用更新后的节点更新原子之间无键时的边。经过3个更新层后，将每个更新层后对应的边连接起来，预测原子间无键时的距离。一个包含3个隐藏层的多层感知器(MLP)用于预测原子的距离[19]，其定义如下：
其中ei表示第i层的输出向量，W是可训练的权重参数，b是偏差，φ是非线性激活函数，以便更好地逼近真实值。这里使用ReLU作为激活函数。
图1 模型架构。(1)预处理：根据输入的SMILES构建分子图和完全图；(2)更新层：通过分子图上的图卷积更新节点，并根据更新的节点更新完整图中的边。(3)边拼接：将每一更新层的边特征向量与上一更新层的节点特征向量进行拼接；(4)MLP1以拼接边特征向量和分子图中对应节点特征向量为输入，MPL2以拼接边结果为输入。
2.3 更新层
更新层中主要目的是对结点和边进行更新，主要操作是图卷积操作，对分子图中的每一个结点，聚合它的邻居节点用于更新自身，聚合操作如下：
其中，N(i)表示结点i的所有邻居节点，l表示图卷积层数，aggregate表示聚合方式，本文使用了求平均值作为聚合方式，得到邻居的聚合信息之后，将其与结点i拼接起来，送入全连接层，输出更新后的i结点特征，公式如下：
其中||表示拼接操作，φ表示激活函数Relu，f表示全连接层。当分子图所有结点更新完毕后，将新的结点特征赋给完全图中的结点，然后更新完全图中原子之间没有键的边，操作如下：
其中，ei表示第i条边的特征向量，和表示形成边的两个结点的特征向量。当所有的边更新完毕后更新层结束，在经过3层这样的更新层后，将每一层得到的边特征向量和最后一层的形成该边的结点特征拼接起来，送入多层感知机中，进行边的长度预测，公式如下：
其中e表示所有原子之间没有键的预测距离，ei表示第i个更新层之后的边特征，和表示最后一次更新层中构成每条边的两个节点的特征向量。
值得一提的是，在每个更新层进行卷积运算后，使用批处理归一化层(BN层)来防止过拟合和梯度消失。
2.4 模型构建和训练
模型是基于PyTorch框架搭建的，为了模型能较快的收敛，将多个分子图组成一个批次，进行训练，采用的BatchSize为2048，一共训练200个epoch。模型使用Adam优化器优化梯度下降，学习率为0.01，后期为了能更好地收敛，将动态调整学习率的值。整个模型的总体损失函数定义为
其中loss1表示MLP1的损失，loss2表示MLP2的损失，经试验选择，α和β分别设置为0.4和0.6，这两部分损失都使用了平均绝对误差，定义为：
其中n为每个batch中边的条数ylabel表示标签，ypredict表示模型计算出的预测值。
模型是在Intel(R) Core(TM) i7-4710MQ CPU上进行训练的，图2是使用小于十一个原子的数据集训练时的损失函数下降曲线，DMGCN模型能在100个epoch内收敛，而随着数据量增大，模型能在更少的epoch内收敛。
图2 损失下降曲线
3 结果
3.1 数据集
在本文中，我们使用QM9数据集来训练和测试我们提出的模型。数据集中的分子由C、N、O、F、H等元素组成，并包含多达9个重原子。该数据集中的原子坐标先由CORINA解析分子的SMILES生成初始的笛卡尔坐标，随后使用MOPAC在PM7半经验理论水平上进行几何松弛，其结果用作高斯的B3LYP的几何松弛的输入[10]，得到优化的原子坐标。数据集中有133885个分子。去掉含有无法处理字符的错误数据的分子后，剩下的131,808个分子作为本文的数据集。
3.2 对比实验
为了研究训练集中分子的最大原子数对模型的影响，我们除了训练整个QM9数据集外，还从QM9中筛选出两个子集对模型进行训练。考虑到原子数的差异和子集中数据量的大小，分别选取原子数小于11和小于15的分子作为子集。我们将数据集按照训练集与测试集的比例为9:1进行划分，并在原子数小于11的子集上手动选择超参数。
此外，为了说明该方法的先进性，将我们的模型与RDKit化学计算包中计算原子几何结构[21]的方法和文献[17]中的DeeperGCN-DAGNN模型进行了比较。据我们所知，DeeperGCN-DAGNN是目前同类方法中最好的。结果如表2所示。
Table 2. RDKit、DeeperGCN-DAGNN和DMGCN在各个子集上的结果 (单位: Å)
RDkit中的方法主要有基于距离的方法和基于知识的方法(ETKDG)。基于距离的方法是通过生成和优化分子边界矩阵，根据边界矩阵随机产生距离矩阵，再将这个距离矩阵映射到三维空间，生成原子坐标，最后使用力场对原子坐标进行粗略的优化。这种方法生成的分子结构较为粗糙，需要再使用力场对其进行优化，而基于知识的方法是根据Riniker等人从晶体结构数据库的小分子结构中总结出的一些规则来修正基于距离的方法得到的结果。
从表2可以看出，随着数据集中原子数的差异越来越大，无论是RDKit中的方法还是我们的模型，产生的误差都在逐渐增加。这是因为随着原子数量的增加，原子之间无键时的距离差异也会增加，使得模型更难拟合。但是可以看出，DMGCN的误差普遍小于RDKit方法，并且随着原子数差异的增大，RDKit的误差增长范围大于本文模型，这也说明该模型比RDKit具有更好的鲁棒性。此外，DeeperGCN-DAGNN在数据量较小时表现较差，但在数据量充足时，MAE介于DMGCN和RDKit之间，RMSE略高于我们的模型。总的来说，DMGCN模型取得了不错的效果。
为了更直观地展示模型的效果，从小于11个原子的测试集中选取误差在平均误差附近的原子间距离预测结果，并与Ground Truth进行比较。如图3所示，原子之间有键时的距离(图3a)相对较小，接近于Ground Truth，而原子之间无键时的距离误差(图3b)相对较大，但趋势与ground truth基本一致。
图3 原子间距离预测；（a）原子间有键时的距离（b）原子间无键时的距离
此外，为了便于观察原子间有键时的距离预测误差分布，我们计算了原子数小于15的子集中原子间有键时距离的绝对误差。统计结果如图4所示。大多数误差集中在0 ~ 0.02Å之间，MAE为0.0208Å, RMSE为0.0301Å。可以看出，该模型在预测原子之间有键时的距离方面更为准确。
图 4. 原子间有键时的距离预测误差分布；红线表示误差的平均值。
此外，表4给出了DMGCN预测的键长、B3LYP/6-31 G (2df, p)计算的键长以及一些来自文献[24]中集合B和QM9数据集[20]交集中的分子上的实验测定的键长之间的误差。表4中，误差1表示模型预测的键长与B3LYP/ 6-31 G (2df, p)方法计算的键长差的绝对值(该方法用于计算QM9数据集[20]中的分子空间结构)，误差2表示B3LYP/ 6-31 G (2df, p)方法计算得到的键长与实验得到的键长之差的绝对值，误差3表示模型预测的键长与实验得到的键长之差的绝对值。在三种误差中，最小值的使用粗体字体，其次是绿色字体。由表4可以看出，B3LYP/ 6-31 G (2df, p)方法得到的键长与实验值非常接近，而且误差3中的小的值要多于误差1，表明多数DMGCN预测的键长比B3LYP/6 - 31 G (df, p)法计算的键长更接近实验值。
表4 DMGCN预测的键长与B3LYP/ 6-31G (2df, p)计算的键长与实验键长之间的两两误差(单位: Å)
3.3 消融实验
模型设计采用BatchNorm层(BN)和双分支结构。此外，为了提高模型的预测精度，我们在节点特征中加入全局信息(分子组成)，边特征采用五维特征，原子之间没有键的边特征是两个原子之间最短路径上的边特征之和。为了验证该模型各部分的有效性，我们在分子原子数小于15的数据集上进行了消融实验。结果如表5所示。可以看出，上面提到的每一部分，特别是双分支和路径最短的边特征，对模型有重要的影响。
表5 消融实验结果。
3.4 性质预测对比
为了表明所提方法的有效性，我们使用文献[25]中的方法，分别根据QM9数据集[20]计算出的分子距离矩阵和我们的方法预测的分子距离矩阵来预测分子的性质。我们将整个数据集按9:1的比例划分为训练集和测试集，并利用训练集对模型进行训练。然后我们在测试集上使用训练好模型来计算分子的距离矩阵,然后将测试集中的数据按照8:2分为训练集和测试集来对文献[25]中的模型进行训练，以预测分子的性质，如最高占据分子轨道的能量(ϵHOMO)、最低非占据分子轨道能(ϵLUMO)、零点振动能(ZPVE)和298.15 K (Cv)时的热容。结果如表6所示。可以看出，基于我们模型预测的分子距离矩阵得到的性质预测结果(Predicted)与基于QM9数据集计算的分子距离矩阵得到的性质预测结果比较接近。说明DMGCN预测的原子间距离可以应用于实际，误差在可接受范围内。
表6 使用QM9中的原子距离矩阵和预测的原子距离矩阵进行性质预测的误差对比
4 结论
本章提出了一种基于图卷积网络的双分支DMGCN模型，以解决传统计算方法在确定分子结构时实验成本高、计算成本高的问题。该模型以分子的SMILES表达式作为输入，通过数据预处理将其转换为图形表示，用于预测分子中任意两个原子之间的距离，达到确定分子结构的目的。
为了展示该方法的进步，将该模型的预测结果与DeeperGCN-DAGNN模型和RDKit化学计算包中的方法进行了比较，得到了更好的MAE。结果还表明，该模型预测的键长与DFT计算的键长和实验测得的键长比较接近。此外，随着数据集中原子数目差异的增大，我们模型的误差增加普遍小于RDKit，具有更好的鲁棒性。随后，为了揭示我们所提出的模型的有效性，我们还进行了分子性质预测实验。结果表明，基于DMGCN计算距离的分子性质预测误差与基于QM9数据集计算距离的分子性质预测误差接近，两者误差均在可接受范围内。
虽然DMGCN的预测结果普遍优于RDKit，但未来仍有很大的改进空间。下一步，我们将考虑如何让模型更全面地了解分子的特征，从而减小原子间无键时的距离预测误差，得到更精确的分子距离矩阵。
作者简介
林小惠（1997-），女，四川广安人，西南交通大学计算机与人工智能学院硕士研究生，主要研究方向为人工智能辅助药物发现。
杨燕（1964-），女，安徽合肥人，西南交通大学计算机与人工智能学院教授，博导，博士，主要研究方向为人工智能、大数据分析与挖掘等。