Datawhale 春训营材料设计赛道篇二

Pythonliu7

已于 2025-05-02 10:16:09 修改

阅读量316

点赞数 5

文章标签：西风胡杨

于 2025-04-27 09:30:52 首次发布

本文链接：https://blog.csdn.net/Pythonliu7/article/details/147547188

版权

概述

随着深度学习技术的发展和定量构效关系(Quantitative structure–activity relationship, QSAR)数据的剧增, 多种深度生成模型联合强化学习或贝叶斯优化算法在生成具有理想属性、新颖且多样的分子上显示出较大潜力, 这为学术界和制药行业大幅缩减了新药研发的时间和成本. 近年来基于深度学习的分子生成方法主要包括基于SMILES和基于分子图表示的方法. 尽管这些分子表示方法在各种应用中都很有效, 但更本质和更富含信息的分子表示是3D几何, 其中原子由其笛卡尔坐标表示. 在化学信息学和计算化学等领域处理3D分子构象至关重要, 这是因为它决定了分子的生物、化学和物理性质, 如电荷分布、势能、分子对接、药效团搜索等. 比如在药物设计中, 理解分子如何与特定的靶蛋白结合至关重要, 而这一过程很大程度上取决于这两者的3D构象, 包括几何(形状匹配)和化学(疏水、亲水)相互作用[2]. 预测3D分子构象的传统方法包括实验测定和基于力场或量子力学的第一性原理预测, 这些方法十分耗时且成本高昂. 随着几何深度学习在处理几何图形数据上的快速发展[3], 近年来面向3D分子构象的深度生成模型逐渐成为研究热点.

分子生成(molecular generation)利用深度学习技术来学习分子数据集的概率分布, 并生成具有理想生物活性和药代动力学属性的新颖分子. 分子生成通过在类药化学空间中搜索和采样多目标属性约束的新分子, 有效减小了搜索空间, 进而加速上述药物发现的进程. 现有的绝大多数基于字符串或基于图的分子生成方法由于缺乏空间信息而无法识别出空间异构体、原子间非键相互作用、以及分子和靶标蛋白在三维空间中的相互作用, 因而在实际药物设计场景中受到限制. 3D分子生成(3D molecular generation)有助于解决上述挑战, 旨在设计出具有理想生物化学属性和几何属性的新颖分子. 比如在基于受体结构的药物设计中, 不同于2D分子图的简化结构表示, 3D分子生成算法生成的3D构象可以在生成过程中考虑分子内或分子间相互作用, 以及来自蛋白口袋的局部约束, 因而能更好地保证预测的结构具有生物活性或药物属性. 根据生成算法的输入将现有的3D分子生成方法总结为三大类, 第1类基于隐变量的生成是指在没有给出参考分子的情况下生成3D分子构象, 即从零开始生成3D分子构象; 第2类基于2D分子图的生成是指给定2D分子图来生成或优化3D分子构象; 第3类基于3D分子构象的生成是指给定3D分子构象来生成或优化3D分子构象. 根据生成算法的输出将3D分子生成方法划分为定向生成和非定向生成. 非定向生成是指利用生成模型学习训练数据集的分布并从中采样出具有新颖、多样、化学有效以及几何有效(稳定态)等通用属性的3D分子构象. 定向生成是指在符合这些通用属性的基础上还需满足多种特定的生物活性和药物属性, 如对特定靶蛋白具有高结合亲和力、靶标选择性、化学可合成性、安全性等. 具体分类以及代表模型如图1和表1所示, 各方法的创新点和优劣性将在第3~5节着重讨论.

3D分子生成是药物发现领域一个极具挑战的任务, 其存在的问题和难点主要包括以下几个方面: (1)搜索空间巨大: 由键长、键角和扭转角组成的三维结构空间巨大, 虽然分子图对可能的三维构象施加了特定的约束(例如键长范围取决于键类型, 四面体中心决定了局部原子的空间排列), 然而可能的构象空间随着图尺寸和可旋转键数量的增加呈指数增长, 这阻碍了穷尽探索的可能性. (2)对称性归纳偏置: 3D分子生成需要处理几何图形式的数据. 不同于一般的图数据, 几何图不仅为每个节点分配一个特征, 还为每个节点分配一个几何向量. 分子可以用几何图来表示, 其中原子的3D坐标是几何向量. 分子的3D构象在空间中具有旋转、平移的几何对称性, 这是因为分子在3D空间中发生旋转、平移后其3D坐标会发生变化, 但控制原子动力学的物理定律是不变的. 在处理分子的3D几何图数据时, 将对称性这种归纳偏置纳入模型设计是至关重要的, 这也推动了几何等变图神经网络的发展. (3)构象分布多模态: 一个分子体系的总能量随其结构的变化而改变的方式可以通过势能面来描述, 势能面上存在多个能量高点和低点, 其中最低处称为全局能量最低点, 其他能量较高的极小点则称为局部能量最低点. 真实的分子倾向于以能量较低的稳定态构象存在, 对应于势能面上的能量低点, 这个数量在一到几千之间变化, 如何有效搜索这些势能面上的局部能量最低点是一大挑战. 最近, 扩散模型凭借其强大的生成能力在建模复杂多模态构象分布上显示出优势, 能更好地生成高质量和多样性的结构. 我们将在下文介绍扩散模型在3D分子生成建模上的出色表现. 现有的机器学习模型正试图从不同角度解决上述挑战, 目前已有一些综述对此领域进行了报道[4~7]. 为了促进更多研究者的参与, 我们对3D分子生成模型从算法到应用的最新进展进行系统回顾. 本文主要内容安排如下: 首先, 概述了3D分子生成的基本概念; 其次, 梳理了基于隐变量的生成、基于2D分子图的生成和基于3D分子构象的生成这三种方法的研究进展, 每一种方法又根据3D分子生成算法的输出分为定向生成和非定向生成进行讨论; 随后, 总结了相关生成模型在主要的公开数据集上的性能, 来说明各种生成模型的优缺点; 最后, 对3D分子生成领域进行总结和展望.

2 分子生成的基本概念

2.1 分子表示

神经网络能将复杂的输入表示转换成解决特定任务所需的潜在空间的表示, 因此输入表示的选择对于模型如何学习分子信息至关重要. 在本节中, 我们将介绍分子表示(molecular representation)的三类主流方法: 1D字符串, 2D分子图和3D分子几何. 图2给出了一个特定分子(青霉烷)的三种表示方法示例.

2.1.1 1D字符串表示

最常见的1D字符串表示类型为简化分子线性输入规范(Simplified Molecular Input Line Entry System, SMILES), 由David Weininger和Arthur Weininger于1986年开发, 并由日光化学信息系统(Daylight Chemical Information Systems)改进和扩展[8]. SMILES是一种简单的基于ASCII字符的表示, 根据预定义的原子排序规则将一个分子转换成一个字符序列, 能够明确地描述原子、键、环、芳香性和立体化学等信息. 同一个分子由于起点原子差异可以表示成多个不同的SMILES, 实际使用时选择统一的规范化SMILES (Canonical SMILES)表示法, 而通用SMILES可用于数据增强. 由于SMILES基于序列表示的特性, 自然语言处理算法可以自然地被移植到分子生成领域. 先前的研究工作[9~12]利用循环神经网络(Recurrent Neural Network, RNN)、长短期记忆网络(Long Short-Term Memory, LSTM)等自回归生成模型生成分子的SMILES字符串. 这些方法忽略了SMILES语法的复杂性, 容易生成无效的SMILES, 即生成的SMILES不符合化合价规则导致不能转换成分子结构. 为了解决这个问题, Krenn等人[13]设计了一种改进的字符串表示方法, 称为自引用嵌入式字符串(Self-Referencing Embedded Strings, SELFIES), 基于SELFIES的RNN模型可以生成100%有效的分子. 另外, 在化学反应预测任务中, 上述针对单一分子设计的方法忽略了反应物和生成物SMILES之间的关系, 我们课题组[14]提出了一种用于化学反应预测的紧凑的分子字符串表示R-SMILES, 通过将输入和输出SMILES的根原子进行对齐的方式, 避免模型学习复杂的SMILES语法, 从而专注于学习化学反应相关的化学知识. 此外, 1D字符串表示还包括国际化学物标识(International Chemical Identifier, InChI)[15]和SMARTS(SMiles ARbitrary Target Specification)[16].

import pickle

# 替换为你的 .pkl 文件路径
pkl_file_path =path to competition/data/competition/competition_smiles.pickle'

# 加载 .pkl 文件
with open(pkl_file_path, 'rb') as file:
    data = pickle.load(file)

# 显示内容
data

['O=C1NCCC1c1ccc(Cl)cc1',
 'C1C2CC12',
 'O=S1(=O)CCN(S(=O)(=O)CC2CCCCC2)CC1',
 'O=C(Cc1cnco1)N1CCNCC1',
 'c1csc(C(c2nnn[nH]2)N2CCCCC2)c1',
 'CC(=O)OC1CNC1=O',
 'O=C(O)CC1NC=NC1=O',
 'O=c1[nH]c(=O)n2ccnc2[nH]1',
 'Cc1cccc(=O)o1',
 'CCOC1C=NOC1',
 'O=C1CCCC/C=C/CCO1',
 'O=C1NCCN1CCc1ccccc1',
 'N=C(Nc1ccc(Cl)cn1)Nc1nccs1',
 'CN(C)/N=C/c1c[nH]c(=O)[nH]c1=O',
 'O=Cc1cc(F)c(F)c(F)c1',
 'CC(O)OC(N)=O',
 'CCN/C(O)=C1/N=CC=CC1=O',
 'NC(=O)Oc1c(F)c(F)c(F)c(F)c1F',
 'NC(=O)N1CCOCC1',
 'CN(C1CCNCC1)S(=O)(=O)c1cccc(Cl)c1',
 'O=c1ocnc2ccc(Cl)cc12',
 'O=[SH](=O)N1CCCNC1',
 'C1=C(N2CCOCC2)CCCC1',
 'COC(=O)NC1=NNC(=O)C1',
 'O=C1N=CSC1F',
...
 'O=C(Cn1cnnn1)N1CCOCC1',
 'O=S(=O)(NCc1ccccc1F)c1ccc(F)cc1F',
 'O=C(c1ccc(Cl)nc1)N(C(O)=S)c1ccc(Cl)cc1',
 'NS(=O)(=O)CC1CCCCC1',
 ...]
Output is truncated. View as a scrollable element or open in a text editor. Adjust cell output settings...

2.1.2 2D分子图表示

表示分子更直观的方法是根据分子的路易斯结构将其转化成一个分子图(molecular graph), 其中原子表示成节点, 化学键表示成边, 并根据原子类型和化学键类型为节点和边赋值对应的标签. 与1D字符串表示相比, 2D分子图表示具有如下优势: 一方面, 图表示能直接捕获原子之间的连接, 而该信息在1D字符串表示中需要模型学习其复杂的语法规则; 另一方面, 如上文所述, 基于SMILES的分子生成方法无法确保生成的分子具有100%化学有效性, 而基于2D分子图的生成方法很容易确保生成的分子是100%化学有效的. 近年来已有许多研究工作致力于开发能生成真实分子图的神经网络架构. 2D分子生成根据生成策略可分为三类: (1)一次输出图的节点类型和邻接矩阵来生成分子图, 有许多深度生成模型可以一次性生成整个分子, 如基于变分自编码器(Variational Autoencoder, VAE)的VGAE[17]、RVAE[18]、GraphVAE[19], 基于生成对抗网络(Generative Adversarial Network, GAN)的MolGAN[20], 基于流(Flow)模型的GraphNVP[21]等. (2)基于子结构生成分子图, 即将子结构或官能团作为构建块来生成高质量的分子. 代表性模型有基于VAE的JT-VAE[22]、HierVAE[23], 基于强化学习(Reinforcement Learning, RL)的RationaleRL[24]、MolEvol[25], 基于采样的MARS[26]、MIMOSA[27]等. (3)通过顺序添加节点和边生成分子图, 代表性模型有基于VAE的CGVAE[28], 基于RL的GCPN[29]、DeepGraphMolGen[30], 基于Flow的GraphAF[31]、GraphDF[32]等.

2.1.3 3D分子几何表示

现有的基于1D字符串或2D分子图的生成模型没有考虑化合物分子的空间构象, 例如同分异构体具有一样的拓扑结构, 但其3D空间结构是不同的. 捕获3D空间结构特征对于准确预测分子的物理化学性质如量子性质、分子构象生成至关重要. 因此, 更本质和更有信息的分子表示是3D几何(3D geometry), 又称3D分子构象(3D molecular conformation), 它包含了3D几何信息(如位置)和重要的物理规律(如对称性). 3D几何将原子坐标加入到它们的表示中, 并能够描述三维空间中原子之间的相对位置. 这里我们考虑处于平衡态的构象, 对应于势能面上的局部最低点. 与1D字符串和2D分子图表示相比, 3D几何表示显著扩展了有待探索的分子空间, 同时也增加了难度, 主要挑战在于模型需要在3D空间中对旋转和平移操作保持等变性. 最近, 3D分子几何表示的研究工作在探索等变性上主要通过以下两种策略, 一种是等变图神经网络, 代表性工作有SE(3)-Transformers[33], EGNN[34]等; 另一种是借助中间几何元素如原子间距离和角度信息来使用E(3)不变的三维信息, 代表性工作有SchNet[2], PhysNet[35], DimNet[36]等. 下文将详细介绍该领域现有的研究进展并对其进行分类整理. 此外, 分子的3D表示还包括3D密度图(3D density map)和3D表面网格(3D surface mesh). 分子的3D密度图由组成分子的所有原子的体积组成[37], 常用3D网格(3D grid)表示, 可用3D卷积神经网络(3D Convolutional Neural Network, 3DCNN)对其进行编码[38]. 3D表面网格由定义分子表面形状的节点、边和面组成, 也常用3D网格表示, 可用测地线卷积神经网络(Geodesic Convolutional Neural Network, GCNN)[39]或3D图神经网络(3D Graph Neural Network, 3DGNN)[2]对其进行编码. 这些分子表示方法在分子生成任务中应用较少, 在这里仅作简单介绍以便读者了解.

2.2 等变图神经网络

2.2.1 等变性

物体(或系统)的对称性(symmetry)是指对其变换操作后, 物体的某一特性保持等变(equivariant)或不变(invariant)的性质[40]. 对称性在物理系统中无处不在, 例如动量守恒中的空间平移不变性, 分子构象、蛋白质或点云的旋转对称性. 在建模物理系统时, 赋予深度学习模型这种归纳偏置对于模型的训练和泛化能力都至关重要. 下面给出等变性的定义和相关符号说明.

E(n)等变(E(n)-equivariance)是指模型在n维欧式空间中对平移、旋转、反射具有等变性[41]. 设有一个函数ϕ:X→Yϕ:X→Y, XX是输入, YY是输出, 定义某变换群GG, 对于任意群元素g∈Gg∈G, 若有

ϕ(Tg(x))=Sg(ϕ(x))(1)

则称函数ϕϕ对变换群GG是等变的, 其中TgTg指变换gg在输入空间中的群表示, SgSg是指同一个变换gg在输出空间中的群表示, 上式说明了对输入的变换操作等同于对输出施加相同的变换. 当Sg=ISg=I, 即输出的群表示是单位矩阵时, 式(1)表示不变性(invariance), 即对输入做群中的任何变换, 输出保持不变. 神经网络ϕ:y=ϕ(x)ϕ:y=ϕ(x)的E(n)E(n)等变性包括以下三种对称性变换, 图3展示了其中的旋转等变性.

(1)平移等变性(translation equivariance): 对输入施加一段平移g∈Rng∈ℝn, 则输出也发生了同样的平移, 有y+g=ϕ(x+g)y+g=ϕ(x+g).

(2)旋转等变性(rotation equivariance): 对输入进行旋转导致输出也发生同样的旋转, 用正交矩阵Q∈Rn×nQ∈ℝn×n表示旋转变换, 则有Qy=ϕ(Qx)Qy=ϕ(Qx).

(3)反射等变性(reflection equivariance): 对输入关于平面(或直线)作反射(翻转)后得到输出, 即输入和输出互为镜像, 用正交矩阵P∈Rn×nP∈ℝn×n表示反射变换, 则有Py=ϕ(Px)Py=ϕ(Px).

2.2.2 图神经网络

图神经网络(Graph Neural Network, GNN)是近年来出现的一种利用神经网络来学习图结构数据, 提取发掘图结构数据中特征和模式的框架, 其优异的性能引起了学者的高度关注和深入探索[42]. 通过在图中的节点和边上制定一定的策略, GNN将图结构数据转化为规范而标准的表示, 并输入到多种不同的神经网络中进行训练, 在节点分类、边信息传播和图分类等任务上取得优良的效果. 当应用于分子时, GNN通过逐步聚合原子及分子环境的信息来提取特征. 下面给出 GNN 中相关符号的说明和定义.

图(Graph)是由节点(Vertex)和连接节点的边(Edge)构成, 记为G=(V,E)G=(V,E). 其中V={v1,v2,...,vn}V={v1,v2,...,vn}表示节点集合, E={e1,e2,...,em}E={e1,e2,...,em}表示边集合. 每个节点ii的节点特征表示为hihi, 连接节点i和j的边特征表示为eijeij. 近年来多种GNN框架被提出[43], 如图卷积神经网络(Graph Convolutional Network, GCN)、图注意力网络(Graph Attention Network, GAT)和图同构网络(Graph Isomorphism Network, GIN)等, 其中最流行的是Gilmer等人[44]提出的消息传递神经网络(Message Passing Neural Network, MPNN). MPNN将GNN统一到了一个通用的框架中, 它迭代地对每个节点或边进行消息传递和邻域聚合. MPNN的消息传播机制表示如下:

mij=ϕm(hi,hj,eij)hi=ϕh⎛⎜⎝∑j∈N(i)mij,hi⎞⎟⎠(2)mij=ϕm(hi,hj,eij)hi=ϕh(∑j∈N(i)mij,hi)(2)

其中mijmij是节点ii经消息传递后生成的信息, ϕmϕm和ϕhϕh分别是消息传递函数和邻域聚合函数, N(i)N(i)是节点i的邻居. 邻域聚合操作ϕhϕh满足置换等变性, 即改变输入顺序, 输出顺序也发生相应改变. 有几项研究工作利用该置换等变性在降低节点排列组合复杂度上的优势, 将GNN应用于具有几何输入的复杂系统, 例如粒子的三维坐标. 如Li等人[45]提出的动态粒子交互网络DPI-Net和Mrowca等人[46]提出的层次化关系网络HRN通过构建动态交互图来模拟粒子动力学. Kipf等人[47]提出的神经关系推断模型NRI通过自动推断潜在交互图, 在多体粒子系统动力学模拟和人体运动捕获方面取得了优异性能. 此外, 在量子化学任务中, GNN通常借助中间几何元素如原子间距离和角度信息来使用E(3)E(3)不变的三维信息, 进而提高平衡态和非平衡态分子构象的力和能量的预测精度, 代表性方法有Schütt等人提出的通用量子化学深度学习架构SchNet[2]和极化原子交互网络PaiNN[3].

当节点特征除了标量还包含实体坐标, 如分子图的节点特征包括原子类型及其3D空间坐标信息时, 还需要考虑原子的3D空间坐标会随着分子在空间中变换而发生变换. 也就是说除了上述讨论的GNN标准置换等变性, 还需要模型满足平移、旋转和反射等变性. 关于E(3)E(3)等变性在分子3D空间中如何发挥作用, 我们将在下文介绍等变图神经网络及其在3D分子构象生成任务中的应用.

2.2.3 等变图神经网络

基于图的方法中, 最近发展的一个领域是E(n)E(n)等变图神经网络(Equivariant Graph Neural Network, EGNN), 它们处理n维图的绝对坐标系. 因此, 这些网络可能适合应用于分子的3D表示. 根据2.2.2对图的定义, 将图表示为G=(V,E)G=(V,E), 具有节点vi∈Vvi∈V和边eij∈Eeij∈E. 在分子3D几何图中, 对于每个节点的特征, 除了之前的节点特征hi∈RHhi∈RH, 还包括节点的三维坐标xi∈R3xi∈ℝ3. Victor等人[34]提出的EGNN不仅能在节点坐标xixi上保持平移、旋转和反射等变性, 还能在节点上保持前文所述的GNN标准置换等变性. 具体来说, EGNN采用了基于标量化建模等变性质的方法, 通过将几何特征转化为不变的标量, 利用多层感知机(Multilayer Perceptron, MLP)得到标量变化, 然后将此标量变化加回到原有的几何特征上从而得到等变性. EGNN的通用框架表示如下:

mij=ϕm(hi,hj,∥xj−xi∥2,eij)xi=xi+C∑j≠i(xi−xj)ϕx(mij)hi=ϕh⎛⎜⎝hi,∑j∈N(i)mij⎞⎟⎠(3)mij=ϕm(hi,hj,‖xj−xi‖2,eij)xi=xi+C∑j≠i(xi−xj)ϕx(mij)hi=ϕh(hi,∑j∈N(i)mij)(3)

其中∥xj−xi∥‖xj−xi‖是几何特征的标量化, ϕϕ函数为不同的MLP, 式(3)的邻域聚合步骤将几何信息和非几何信息进行了关联, 保证了几何特征和非几何特征在消息传播过程中的等变性.

2.3 模型框架

分子生成的模型框架主要有两类. 第一类是深度生成模型, 它旨在利用深度学习模型对数据分布进行建模. 第二类是利用启发式策略来实现目标的组合优化方法. 两者在利用结构化数据上的主要区别在于, 深度生成模型利用的是连续隐空间的表示, 而组合优化方法直接在结构化数据空间中进行搜索. 在这里我们简要描述这些模型框架的原理, 以及它们在3D分子生成上的应用.

2.3.1 深度生成模型

深度生成模型通过学习可观测数据的概率密度来生成新的合成数据, 包括图像、视频和文本等. 这里主要关注药物发现中3D分子生成的深度生成模型. 深度生成模型的目标是最小化生成数据分布与真实数据分布之间的距离, 可以用极大似然法进行求解. 从处理极大似然函数方法的角度, 可将生成模型分为三类: 第一类方法通过变分或抽样的方法求似然函数的近似分布, 主要包括直接优化似然函数变分下界的变分自编码器和采用抽样方法近似计算似然函数的玻尔兹曼机; 第二类方法是避开求极大似然过程的隐式方法, 其代表模型是通过生成器和判别器之间的对抗行为来优化模型参数从而巧妙避开求解似然函数的生成对抗网络; 第三类方法是对似然函数进行适当变形以简化计算, 包括流模型和自回归模型, 流模型利用可逆函数构造似然函数后直接优化模型参数, 自回归模型将目标函数分解为条件概率乘积的形式. 下面将分别介绍这三类模型的代表性模型框架: 变分自编码器、生成对抗网络、流模型、自回归模型、扩散模型, 如图4所示.

(1) 变分自编码器

变分自编码器(Variational Autoencoder, VAE)[48]的基本结构包括一个编码器和一个解码器, 其中编码器学习将分子映射到连续低维空间的隐变量, 然后解码器将隐变量重建为新分子. 训练VAE时最小化的损失函数包括重构项和正则化项两部分, 重构项迫使解码器从隐变量中重建出正确的分子, 正则化项则使用Kullback–Leibler散度(KL散度)来衡量重建的样本分布与数据分布之间的距离. 在3D分子生成中, Mansimov等人[49]提出的条件变分图自编码器(CVGAE)首次以2D分子图为条件, 使用VAE学习原子表征并直接生成原子坐标, 但不能保证构象的旋转平移等变性, 需将生成后的构象与真实构象进行比对, 性能较差. 为了解决这个问题, Simm等人[50]提出了图距离几何模型(GraphDG)利用中间几何元素(键长)来生成构象. 具体来说, GraphDG使用VAE预测键长, 再通过距离几何算法利用先前计算得到的键长来求解构象. Xu等人[51]提出的ConfVAE模型使用条件VAE和距离几何算法将GraphDG中分离的两个步骤结合起来, 改进成端到端的双层优化框架.

(2) 生成对抗网络

生成对抗网络(Generative Adversarial Network, GAN)[52]是当前深度学习领域最热门的生成模型, 广泛应用于图像生成领域. GAN本质上是将难以求解的似然函数转化成神经网络, 让模型自己训练出合适的参数来拟合似然函数, 这个神经网络就是GAN的判别器. GAN由一个生成器和一个判别器组成, 生成器的目标是从高斯噪声中生成逼真的分子让判别器无法判别出真伪, 判别器的目标是正确区分分子是真实的(来自训练数据集)还是虚假的(由生成器生成). 这两个网络在训练过程中互相博弈, 不断优化自身的生成能力和判别能力, 最终达到两者之间的纳什均衡, 即当判别器的判别能力达到一定程度却无法判断数据来源时, 就获得了一个学习到真实数据分布的生成器. 先前已有一些研究工作使用GAN来生成分子的SMILES字符串[53]或分子图[18,20]的邻接矩阵. 最近, 研究人员采用这种方法来生成3D分子构象, 如Hoffmann等人[54]提出的欧几里得距离矩阵网络(EDMNet)使用GAN来生成原子间距离矩阵, 进而生成分子3D结构. 然而, 使用GAN生成序列、分子图或3D构象仍然具有挑战性, 主要体现在以下几点: (1) 模型难以训练, 构建序列、图或3D坐标需要离散选择反向传播梯度, 容易出现梯度消失; (2) 出现模式崩溃, 表现为生成的分子多样性低; (3) 目标函数的形式导致模型在训练过程没有可以指示训练进度的指标.

(3) 流模型

主流深度生成模型中, VAE推导出了似然函数的变分下界来代替真实的数据分布, 这种近似模型并不能生成最精确的分子; GAN使用对抗和交替训练的方法避免了优化似然函数, 保留了模型的精确性, 但在训练过程中会出现训练困难、模式崩塌等问题, 因此研究一种既能保证模型精度又容易训练的深度生成模型具有重要意义. 标准化流(Normalizing Flow)[55]通过学习先验分布(如高斯分布)和真实世界的高维数据(如分子)之间一系列可逆变换来生成分子, 同时该可逆变换函数可计算得到. 流模型相对于VAE的主要优势在于其可逆映射能够计算精确的数据似然, 有望生成质量足够好的分子. 得益于该优势, 近年来已有一些研究将流模型应用于分子图的生成[21,31,56]. 如Shi等人[31]提出了基于流的自回归模型GraphAF来生成分子图, 其性能优于VAE. 最近, 研究人员采用流模型来生成3D分子构象, 如Satorras等人[57]提出的E(n)E(n)等变标准化流模型(E-NFs)将EGNN集成到标准化流模型, 能一次性同时生成所有原子的特征和3D坐标. Luo等人[58]提出了生成3D分子几何结构的自回归流模型G-SphereNet, 将分子生成建模成一个序列决策过程, 根据条件信息逐步生成原子的类型和相对位置信息(距离、角度和扭转角), 保证了模型的不变性和等变性.

(4) 自回归模型

自回归(Autoregressive, AR)模型可以对序列数据进行密度估计, 用同一变量在之前各个时刻的观测值预测该变量在当前时刻的观测值. AR将目标函数分解为条件概率乘积, 即首先用条件概率表示可见层数据相邻元素之间的关系, 然后以条件概率乘积来表示联合概率分布. 分子图数据可以转化成一系列节点和边的序列数据, 这些节点和边作为分子图的子组件存在潜在的依赖关系. AR将分子图数据的联合分布分解为一系列子组件的概率乘积, 可通过前面的子组件预测下一个子组件. AR最大的优势是可以对序列数据进行很好的密度估计, 在图像、文本和语音等领域的生成任务中得到成功应用, 包括循环神经网络(RNN)[59], PixelCNN[60], WaveNet[61], Transformer[62]和BERT[63]等. 近年来, AR也被广泛应用于3D分子生成领域, 如Gebauer等人[64]提出的G-SchNet模型将欧氏空间约束和原子分布的旋转不变性作为先验知识, 使用基于SchNet[65]的自回归神经网络模型依次生成新原子, 并将其放置在焦点原子的局部网格点. Joshi等人[66]提出的3D-Scaffold模型建立在G-SchNet基础上, 从现有的3D骨架中自回归地生成新原子, 生成的新候选分子保证具有所需的骨架. 这种方法需要分子骨架和目标性质之间关系的先验知识, 并且可能会因为对特殊分子的采样而影响模型性能. 该论文的作者在此基础上又提出了条件G-SchNet (cG-SchNet)[67], 能够以多种结构和化学特性的任意组合为条件, 有针对性地生成3D分子结构. 然而, 使用AR生成序列、图或3D构象仍然具有挑战性, 由于AR的每一步输出都需要按顺序等待前一步的输出, 这种按顺序生成的方式会阻碍AR利用并行运算, 导致训练和生成样本所需要的计算量远大于VAE和GAN等其他通用模型. 实际使用时可采用数据并行来利用并行计算资源从而加速生成.

(5) 扩散模型

扩散模型(Diffusion models)[68]的提出受到物理学中非平衡热力学的启发, 它定义了一个扩散步骤的马尔科夫链, 可分为正向扩散过程和逆向扩散过程. 扩散模型在正向扩散过程中逐步对数据添加随机噪声, 直至最后一步达到各向同性高斯分布; 逆向扩散过程则反转这一过程, 从噪声分布中构建所需的数据样本. 得益于生成真实样本的强大能力, 扩散模型在建模复杂的多模态分子构象分布上显示出优势, 能更好地生成高质量和多样性的结构. 如Shi等人[69]提出的ConfGF是基于扩散模型的分子构象生成早期工作, 在保证旋转和平移等变性的同时, 利用朗之万动力学方法对分子构象进行采样. 然而, ConfGF只建模了一阶、二阶和三阶邻居之间的局部距离, 忽略了分子内部的长程非键相互作用力. 由此, Luo等人[70]提出的动态图评分匹配(DGSM)模型通过动态地建模三维空间中原子对之间的相互作用力, 包括化学键和非键相互作用, 将ConfGF扩展到了长程作用力的建模. 然而, DGSM在扩散学习过程中输入的扰动距离矩阵可能违反三个原子之间距离的三角不等式. 为了解决这一挑战, Xu等人[71]提出的几何扩散模型(GeoDiff)引入了旋转–平移不变马尔可夫过程来对坐标分布施加约束. Hoogeboom等人[72]提出的E(3)E(3)等变扩散模型(EDM)进一步扩展了上述方法, 通过一个等变网络学习去噪扩散过程, 能同时对连续特征(原子坐标)和分类特征(原子类型)进行操作. 由于构象生成的难点主要在于扭转角的自由度, Jing等人[73]提出的扭转扩散模型(Torsional Diffusion)将超环面的扩散过程和一个从外向内的评分模型作用在扭转角空间来生成分子构象, 显著提升了预测性能和速度.

2.3.2 组合优化方法

深度生成模型需要大量的标注数据进行训练, 但药物发现场景难以满足这种数据要求, 因为标注数据需要资源和时间密集型的生化实验. 在这种情况下, 组合优化方法可以通过将分子生成问题转化成离散优化问题来帮助缓解这种数据稀疏的现状, 但难点在于需要在巨大的化学空间中找到探索和利用的平衡. 在这里我们主要讨论5种组合优化方法: (1) 强化学习; (2) 贝叶斯优化; (3) 遗传算法; (4) 蒙特卡洛树搜索; (5) 最优传输.

(1) 强化学习

强化学习(Reinforcement Learning, RL)通常被描述为马尔可夫决策过程(Markov decision process, MDP), 用来建模智能体与环境的交互过程. MDP的公式MDP=(S,A,P,R)MDP=(S,A,P,R)描述了RL中的几个基本概念: 状态集合SS、动作集合AA、状态转移概率PP和奖励函数RR. 状态集合SS表示所求解问题中所有可能出现的状态所构成的有限或无限的集合; 动作集合AA表示所求解问题中智能体能够采取的所有动作所构成的有限或无限的集合; 状态转移概率P(St+1|St,At)P(St+1|St,At)表示在状态StSt下采取了动作AA后进入下一个状态St+1的概率, 满足马尔可夫性; 奖励函数R(St,At,St+1)R(St,At,St+1)表示在状态StSt执行了动作AA后到达状态St+1St+1时, 智能体能够得到的奖励. RL的目标在于学习一种能使奖励最大化的最优策略. 在2D分子生成任务中, 状态通常被定义为部分生成的分子图, 动作可以是在某个位置添加一个原子、化学键或子结构, 奖励可以是部分生成的分子的某种属性. 如You等人[29]提出的图卷积策略网络(GCPN)和Zhou等人[74]提出的分子深度Q网络(MolDQN)使用强化学习在分子图上迭代地添加或删除原子和化学键以最大化类药性(Quantitative Estimate of Drug-likeness, QED)分数, 从而生成具有高QED的类药化合物. Olivecrona等人[75]在SMILES字符串上操作, 使用RL来微调先前预训练的生成模型, 以最大化生成分子与多巴胺受体的结合亲和力. 近期, 利用强化学习进行3D分子结构从头设计也取得了很大进展. Simm等人[76]提出的MolGym模型使用强化学习在笛卡尔坐标下从头生成分子. MolGym利用用量子力学性质能量作为奖励函数, 使用基于平移和旋转不变状态–动作表示的行动者–批评神经网络来保证模型的平移旋转等变性质. 该模型在MolGym RL环境的一系列分子设计任务中能生成有效的分子. 在此基础上, 作者进一步提出了一种行动者–批评神经网络COVARIANT[77], 通过基于球面谐波级数展开的旋转协变状态–行动者表示保证了分子设计过程的对称性, 结果表明利用分子设计过程的对称性显著提升了生成分子的质量和模型泛化能力.

(2) 贝叶斯优化

贝叶斯优化(Bayesian Optimization, BO)是一种调整超参的技术, 常用于优化计算成本高或不可导的目标函数. BO通过代理优化的方式, 根据新的信息来更新前信念以产生一个更新的后信念. BO有两个关键步骤: 1) BO为目标函数建立代理函数, 常用高斯过程表示的概率分布来量化代理函数的不确定性; 2) 采集函数在由代理函数确定的可能最小值附近做更多采样, 并随之迭代更新代理函数, 最终达到全局最小值. BO[8,22,78]通过学习分子的连续嵌入, 将离散优化问题转化为连续优化问题. 这些方法包括首先训练一个VAE将离散的分子映射到连续的嵌入空间, 然后训练另一个神经网络从其连续嵌入向量中预测原分子的化学属性. 然后, 他们在连续嵌入空间中使用BO来找到一个关联属性得分最好的嵌入, 该嵌入进而被解码器解码成一个离散的分子. 在3D分子构象采样中, 一个核心挑战是使用较少的能量评价指标来找到低能稳定构象, Chan等人[79]提出了一种叫做贝叶斯优化算法(BOA)的随机搜索方法来寻找分子的低能稳定构象. 与穷举搜索或随机搜索方法相比, BOA由于较好平衡了利用和探索之间的关系, 能使用更少的能量评价指标来找到更低的能量构象. 该论文作者在此基础上又提出了基于先验知识的贝叶斯优化模型BOKEI[80], 使用相关扭转角知识作为约束条件, 将低能构象搜索过程定义为一个约束BO问题, 与基准相比该方法能在减少搜索空间的同时找到更低能量构象.

(3) 遗传算法

遗传算法(Genetic Algorithm, GA)是借鉴生物的自然选择和遗传进化机制而开发出的一种全局优化自适应概率搜索算法. GA使用群体搜索技术, 通过对当前群体施加选择、交叉、变异等一系列遗传操作, 从而产生出新一代的群体, 并逐步使群体进化到包含或接近最优解的状态. 由于其具有思想简单、易于实现、应用效果明显等优点而被众多应用领域所接受, 并在自适应控制、组合优化、模式识别、机器学习、人工生命、管理决策等领域得到广泛应用. 在分子生成任务中, GA通过分子突变来寻找具有理想属性的化合物分子. GA的核心由两部分组成: 一组变异规则和一个适应度函数. 适应度函数包括一个加权插值的属性预测分数和一个用于惩罚长期存活分子的惩罚分数, 其中惩罚分数鼓励模型探索多样的分子. 对现有分子应用突变规则得到新分子, 在每次迭代中, 适应度分数较低的旧分子会被移除. 如Nigam等人[81]通过神经网络参数化适应度函数, 使用遗传算法设计出高油水分配系数(logP)的分子. 不同于Nigam等人使用一组固定的突变规则, Ahn等人[82]提出的遗传专家指导学习(GEGL)框架使用策略函数执行分子突变来生成高回报的分子, 进一步拓展了GA在分子生成上的应用. GA也被广泛应用于3D分子生成领域, 如Spiegel等人[83]提出的半自动化计算机辅助药物设计开源程序AutoGrow4利用遗传算法和反应库来对种群中的分子进行突变, 迭代应用虚拟反应来最大化期望适应度函数来生成或优化药物分子. Fu等人[84]提出的增强遗传算法(RGA)将遗传算法中的进化过程建模成MDP, 用神经网络代替随机游走方法来指导遗传算法中的交叉和变异操作, 以更智能地探索化学空间. 此外, 该网络使用多种配体–靶标复合物进行预训练来学习不同靶标共享的结合亲和力物理特性, 然后在优化过程中用强化学习进行微调.

(4) 蒙特卡洛树搜索

蒙特卡洛树搜索(Monte Carlo Tree Search, MCTS)方法在分子优化方面显示出优越性. 不同于GA, MCTS可以从头生成分子. MCTS包括四个步骤: 选择、扩展、模拟和反向传播. 具体来说, MCTS通过迭代以上四个步骤来构建搜索树, 树中每个节点代表分子生成过程的一个中间状态. 每次迭代中, 模型首先从搜索树中选择一个有希望的状态, 枚举该状态的可能操作, 并模拟生成原子、化学键或子结构来得到下一个中间状态, 通过收集奖励信息并反向传播. 如Li等人[85]提出的DeepLigBuilder模型能直接在靶蛋白的3D结合口袋内构建和优化配体的3D分子结构. 该模型使用配体神经网络(L-Net)生成有效的3D类药分子, 用MCTS搜索高结合亲和力的分子, 可以直接输出类药物分子的3D结构和拓扑结构, 而不需要额外的原子位置或键序推断.

(5) 最优传输

最优传输(Optimal Transport, OT)是一种比较概率分布的几何工具, 即利用几何方法来建模概率分布和衡量概率分布之间的距离, 是连接几何和概率的桥梁. OT的核心是研究如何用最经济的方式将一个概率分布变换到另外一个概率分布. 由于深度神经网络本质上是在训练概率分布间的传输映射, 常需比较两个概率分布间的距离, OT为此提供了解决思路. 近年来, OT被广泛应用于深度学习领域, 如Arjovsky等人[86]提出的WGAN利用生成器计算OT映射, 将隐空间的均匀分布映射成数据分布, 判别器计算数据分布和生成分布间的Wasserstein距离. 最近OT被应用于3D分子生成, 如Ganea等人[87]提出的分子3D构象扭转几何生成模型GeoMol是一种端到端、非自回归和SE(3)SE(3)不变的机器学习方法, 用于生成低能构象异构体分布. GeoMol利用MPNN捕获局部和全局图信息, 预测局部原子3D结构和扭转角来组装相邻的局部3D结构对, 避免由于几何自由度导致不必要的过度参数化. 模型在局部预测中使用基于非对抗性OT的损失函数来促进多样化构象异构体的生成, 最后测试时将所有预测的相邻局部3D结构对组合成完整构象. Stärk等人[88]提出的SE(3)SE(3)等变几何深度学习模型EquiBind利用E(3)E(3)协变网络生成配体和受体上的结合关键点, 而后在三维空间中使用OT损失对这两组点进行匹配从而计算其差异. 利用上述得到的关键点, 模型将输出的小分子构象叠合至蛋白口袋中来确定结合模式. 同时, 该网络还能够对配体的三维结构进行变换实现小分子柔性对接.

3 基于隐变量的生成

基于隐变量的生成是指在没有给出参考分子的情况下在3D空间中生成分子, 即从零开始生成3D分子构象, 生成的构象需同时满足化学有效与几何有效, 以及多种理想属性, 包括生物活性、靶标选择性、结构新颖性、化学可合成性、成药性和安全性等. 根据3D分子生成算法的输出, 可将基于隐变量的生成任务分为定向生成和非定向生成. 表1从各方面对这些方法进行了比较.

3.1 非定向生成

该方法的目标是在没有给出参考分子的情况下从给定数据集中学习3D分子构象分布, 并从中采样出新颖、多样、化学有效以及几何有效的3D分子构象.

Gebauer等人[64]提出的G-SchNet模型将欧氏空间约束和原子分布的旋转不变性作为先验知识, 使用基于SchNet[65]的自回归神经网络模型依次生成新原子, 并将其放置在焦点原子的局部网格点. 作者将生成模型的分布偏向于具有小HOMO–LUMO间隙的分子, 表明模型能捕捉到3D几何结构与电子特性之间的复杂关系. Hoffmann等人[54]提出的欧几里得距离矩阵网络(EDMNet)和Nesterov等人[89]提出的3D分子结构生成模型3DMolNet分别使用GAN和VAE来生成原子间距离, 进而生成分子3D结构. Satorras等人[57]提出的E(n)E(n)等变标准化流模型(E-NFs)将EGNN集成到Flow模型来保证生成模型的对称性. 具体来说, 通过将EGNN集成到连续时间可微的神经常微分方程获得Flow的可逆等变函数. 该模型能一次性同时生成所有原子的特征和3D坐标, 生成的分子结构在准确性和稳定性上大大优于现有的Flow模型. Luo等人[58]提出的生成3D分子几何结构的自回归流模型G-SphereNet, 将分子生成建模成一个序列决策过程, 根据条件信息逐步生成原子的类型和相对位置信息, 保证了模型的不变性和等变性. 此外, 在序列生成过程中, 使用球形消息传递网络SphereNet和注意机制来提取条件信息. G-SphereNet在随机分子几何结构生成和以量子性质为目标的分子发现任务方面优于E-NFs和G-SchNet. 近期, 利用RL进行3D分子结构从头设计也取得了很大进展. Simm等人[76]提出的MolGym模型使用RL在笛卡尔坐标下从头生成分子. MolGym将量子力学性质能量作为奖励函数, 使用基于平移和旋转不变状态–动作表示的行动者–批评神经网络来保证模型的平移旋转等变性质. 该模型在MolGym RL环境的一系列分子设计任务中能生成有效的分子. 在此基础上, 作者进一步提出了一种行动者–批评神经网络COVARIANT[77], 通过基于球面谐波级数展开的旋转协变状态–行动者表示保证了分子设计过程的对称性, 结果表明利用分子设计过程的对称性显著提升了生成分子的质量和模型泛化能力.

3.2 定向生成

该方法的目标是在没有给出参考分子的情况下从头生成出同时满足化学有效和几何有效, 以及具有多种理想药物属性或能与特定蛋白口袋结合的3D分子构象.

Masuda等人[90]提出的一种基于条件变分自编码器(CVAE)与GAN的分子生成模型LiGAN, 能够在蛋白口袋中生成分子, 主要的策略是先生成分子的原子密度图, 再将其转换成分子结构. 该实验室的研究者们在此基础上又提出了一步生成基于受体结合位点的3D分子结构[91]. 具体来说, 首先使用CVAE来学习蛋白质–配体复合物的原子密度网格表示, 然后用原子拟合和键推断的方法从生成的原子密度构建出有效的分子构象. Li等人[85]提出的DeepLigBuilder模型使用深度生成模型直接在靶蛋白的3D结合口袋内构建和优化配体的3D分子结构. 该模型提出了一个配体神经网络(L-Net)来生成有效的3D类药物分子, 用MCTS搜索强结合分子, 可以直接输出类药分子的3D和拓扑结构, 而不需要额外的原子位置或键序推断. Luo等人[92]提出的DSBDD模型将给定的蛋白结合位点作为3D上下文, 估计原子在3D空间中出现的概率密度, 即更有可能出现原子的位置将被赋予更高的概率, 然后使用自回归采样方案从学习的分布中顺序地采样原子, 直到没有空间容纳新的原子. 结合该采样方案, 该模型可以生成有效且多样的分子. 实验结果表明, 即使模型没有明确地被优化, 从该模型中采样的分子也表现出与特定靶点的高亲和力和良好类药性.

4 基于2D分子图的生成

不同于上述基于隐变量的方法, 基于2D分子图的生成是指从给定的2D分子图生成或优化出同时满足化学有效和几何有效, 以及具有多种理想的生物活性和药物属性的3D分子构象. 类似于基于隐变量的生成方法的子分类, 根据3D分子生成算法的输出, 上述每一大类生成方法可二次划分为定向生成和非定向生成. 表1按照不同的算法和生成方法进行了组织和分类.

4.1 非定向生成

该方法的目标是从给定的2D分子图生成新颖、多样、化学有效以及几何有效的3D分子构象. 从2D分子图生成3D分子构象对许多下游任务如受体–配体结合亲和力预测、基于受体结构的药物设计等具有重要意义.

Mansimov等人[49]提出的条件变分图自编码器(CVGAE)首次以2D分子图为条件, 使用VAE模型学习原子表征并直接生成3D原子坐标, 但不能保持构象的旋转平移等变性, 需将生成后的构象与真实构象进行比对, 性能较差. 为了解决这个问题, 大多数后续的模型都是基于中间几何元素(如原子间距离和扭转角)来生成构象. 这些元素在分子平移、旋转时能保持不变, 这种平移旋转不变性已被证明是分子几何建模中的一个重要的归纳偏置. 如Simm等人[50]提出的图距离几何模型(GraphDG)和Xu等人[93]提出的条件图连续流模型(CGCF)分别使用VAE和Flow模型预测原子间距离, 再通过距离几何算法利用原子间距离来求解构象. Xu等人[51]提出的条件变分自编码器框架ConfVAE使用条件VAE和距离几何算法将GraphDG、CGCF中分离的两个步骤结合起来, 改进成端到端的双层优化框架. 但这些方法都存在累积误差问题, 即预测原子间距离过程的噪声会传导到后续的坐标搜索过程, 影响生成的3D结构的准确性. 为了克服上述缺陷, Shi等人[69]提出的ConfGF模型受到经典计算化学中力场的启发, 使用去噪得分匹配算法来学习原子间距离的对数密度梯度场, 并通过链式法则计算出原子坐标, 然后利用朗之万动力学方法对分子构象进行采样. 然而, ConfGF忽略了分子内部的长程非键相互作用力, 比如分子图中相距较远的两个原子可能在三维结构中形成氢键. Luo等人[70]提出的动态图评分匹配(DGSM)模型通过动态地建模三维空间中原子对之间的相互作用力, 包括化学键和非键相互作用, 将ConfGF扩展到了长程作用力的建模. 这种基于中间几何元素方法的优势是能保证旋转和平移的不变性. 然而, 预测的原子间距离可能会相互冲突, 违反三个原子之间距离的三角不等式, 例如有研究发现GraphDG在GEOM–Drugs测试集上对8.65%的分子产生了违反三角不等式的距离矩阵[94]. 此外, 最新最前沿的方法是最近提出的GeoMol[87], GeoDiff[71]和Torsional Diffusion[73]. GeoMol首先学习预测几何元素的最小集合(即键长和键角)得到局部3D结构, 然后通过预测扭转角来组装相邻的局部3D结构对, 最后在测试时将所有预测的相邻局部3D结构对组合成完整构象. 为了生成不同的构象, 随机高斯噪声被附加到每个初始节点, 并使用基于OT的损失函数对模型进行训练. GeoDiff模拟扩散过程将分子坐标扰乱成噪声分布, 再用去噪扩散模型模拟逆向扩散过程来恢复真实构象. 作者将去噪扩散模型与几何表示相结合, 将逆向扩散过程建模成马尔可夫链, 并将平移旋转等变性引入马尔可夫核, 以保证模型的平移旋转等变性质. 由于构象生成的难点主要在于扭转角的自由度, Jing等人[73]提出的扭转扩散模型(Torsional Diffusion)将扩散过程只作用于扭转角, 而键长和键角则通过标准的化学信息学方法进行计算. Torsional Diffusion成为首个在GEOM–Drugs数据集上效果超越商业软件OMEGA[95]的深度生成模型, 并且在训练速度上比GeoDiff快了两个数量级.

4.2 定向生成

该方法的目标是从给定的2D分子图生成或优化出同时满足化学有效和几何有效, 以及具有多种理想药物属性或能与特定蛋白口袋结合的3D分子构象.

Méndez-Lucio等人[96]提出的DeepDock模型利用几何深度学习方法来预测配体和靶蛋白结合构象. 该模型通过学习每个配体–靶体对基于距离似然的统计力学势, 使用能量模型(EBM)作为全局优化算法对其进行优化, 从而构建出配体的实验结合构象. 同已有的用于对接和筛选任务中的评分函数相比, 此处基于距离似然的统计势有类似或更好的效果, 这种方法代表了人工智能用于改进基于结构的药物设计的一个案例. Stärk等人[88]提出的SE(3)SE(3)等变几何深度学习模型EquiBind利用E(3)E(3)协变网络生成配体和受体上的结合关键点, 而后在三维空间中对这些点进行匹配来确定结合模式. 同时, 该网络还能够对配体的三维结构进行变换实现小分子柔性对接.

5 基于3D分子构象的生成

基于3D分子构象的生成是指从给定的3D分子生成或优化出同时满足化学有效和几何有效, 以及具有多种理想的生物活性和药物属性的3D分子构象. 根据生成算法的输出, 二次分类为定向生成和非定向生成. 表1按照不同的算法和生成方法进行了分类比较.

5.1 非定向生成

该方法的目标是从给定的3D分子生成或优化出新颖、多样、化学有效以及几何有效的3D分子构象.

针对两个分子片段的连接问题, Huang等人[97]提出了一种基于条件变分自编码器的生成模型3DLinker, 利用EGNN来结合分子的不变特征与等变特征, 将分子图的生成与分子三维信息的表征进行整合, 从而有效地预测连接器的锚点与边的具体连接状态. 该模型在分子片段连接器的恢复上表现出良好的性能优势, 并能够准确预测分子的3D坐标信息.

5.2 定向生成

该方法的目标是从给定的3D分子生成或优化出具有理想生物活性或药物属性的3D分子构象, 如低能分子构象优化、基于骨架的分子构象优化、由粗到细的分子构象优化等.

Chan等人[79]提出了一种叫做贝叶斯优化算法(BOA)的随机搜索方法来搜索分子的低能稳定构象. Wang等人[98]提出的粗粒化变分自编码器(Coarse-Graining VAE)使用VAE模型学习从粗粒度分子构象生成细粒度分子构象. Joshi等人[66]提出的3D-Scaffold模型建立在G-SchNet基础上, 从现有的3D骨架中自回归地生成新原子, 生成的新候选分子保证具有所需的骨架. 这种方法需要分子骨架和目标性质之间关系的先验知识, 并且可能会因为对特殊分子的采样而影响模型性能. 该论文的作者在此基础上又提出了条件G-SchNet (cG-SchNet)[67], 能够以多种结构和化学特性的任意组合为条件, 有针对性地生成3D分子结构. 与以前的方法相比, 该模型不需要在特定目标属性的数据集上微调, 显式条件使cG-SchNet能够从所有可用的参考数据集中有效地学习. Fu等人[84]提出的增强遗传算法(RGA)将遗传算法中的进化过程建模成MDP, 用神经网络代替随机游走方法来指导遗传算法中的交叉和变异操作, 以更智能地探索化学空间. 此外, 该网络使用多种配体–靶标复合物进行预训练来学习不同靶标共享的结合亲和力物理特性, 然后在优化过程中用强化学习进行微调. Gogineni等人[99]提出了一种基于RL的顺序构象搜索技术TorsionNet, 将3D构象生成问题建模成MDP, 使用基于热力学的吉布斯分数作为奖励, 使用近端策略优化(PPO)算法和课程学习策略(curriculum learning strategy)来训练模型. TorsionNet是第一个用于大分子构象搜索的深度学习模型, 其在大支链烷烃和木质素聚合物数据集上的性能明显优于经典的化学信息学和分子动力学采样方法.

6 3D分子生成算法评价

6.1 常用数据集

近几十年来, 3D分子数据库得到不断积累. 这里我们主要关注基于深度学习的3D分子设计任务相关的公开数据集, 并将其分为以下两类: 第一类化学空间数据集是指以类似于列举化学空间的方式包含特定原子组成的化合物; 第二类任务数据集用于药物设计中的特定任务, 如分子动力学模拟, 能量、力等量子属性预测和亲和力预测. 第3、4、5节中讨论的大量工作已经证明了这些数据集的有效性, 在准确预测或生成目标的同时大大降低了传统方法密度函数理论(DFT)所需的计算量. 下面列出这些数据集, 并在表2展示其统计信息.

6.1.1 化学空间数据集

(1) GEOM[100]: 一个用于3D分子生成的基准数据集, 包含了超过430000个分子的数百万个构象, 由第一性原理模拟得到. GEOM被拆分成了GEOM-QM9和GEOM-Drugs这两个子集. GEOM-QM9是QM9数据集的扩展, 它仅限于含有9个重原子(总共29个原子)的小分子, 分子质量小, 可旋转键少. GEOM-Drugs数据集由中等大小的药物分子组成, 最多181个原子(91个重原子), 平均包含44.2个原子, 每个分子平均有6.5个可旋转键.

(2) Molecule3D[101]: 一个用于3D分子生成的基准数据集, 包含大约400万个分子的精确基态几何构象, 由DFT计算得到. 该基准数据集提供了一套软件工具用于数据处理、拆分、训练和评估等操作.

(3) QM9[102]: 一个用于量子属性预测的基准数据集, 来自GDB17数据集, 包含133885个有机小分子, 由不超过9个重原子(C, O, N和F)组成. 每个分子都具有几何、能量、电子和热力学等12种量子属性, 由DFT计算得到.

6.1.2 任务数据集

(1) ISO17[2]: 一个分子动力学轨迹数据集, 包含了从QM9中提取的原子组成为C7O2H10的129个同分异构体的分子动力学轨迹, 每条轨迹包含5000个构象, 数据集由64.5万个带标记的实例组成. 其中的分子反映了平衡态周围的构象分布, 可用于评估模型分布与平衡态周围数据分布之间的相似性.

(2) MD17[103]: 一个分子动力学轨迹数据集, 包含8个小分子的分子动力学轨迹及其对应状态的能量和相互作用力信息. 每个数据都由单个分子的轨迹组成, 覆盖了各种构象.

(3) Atom3D[104]: 一个涵盖小分子、蛋白质和核酸三维结构的综合数据集, 包含了8个预测任务涉及分子结构、分子之间相互作用、分子功能和设计等.

(4) PDBbind[105]: 是蛋白质数据库(PDB)的一个子集. PDBBind v2020包含19443个蛋白–配体复合物结构及其各自的势能(半峰抑制浓度、解离常数等), 其中包含3890个独特受体和15193个独特配体.

(5) CrossDock2020[38]: 一个利用Pocketome数据库和smina对接程序对PDBbind进行扩增得到的数据集, 包含总计22500000个蛋白-配体的结合构象.

(6) scPDB[106]: 一个蛋白–配体结合位点数据库, 共16034个条目, 包含由3678个独特蛋白质和5608个独特配体组成的9283个结合位点. 该数据库提供了蛋白质、配体、结合位点和结合模式的全原子描述.

(7) DUD-E[107]: 一个广泛用于分子对接和虚拟筛选的数据集, 包含22886种活性化合物及其对102个靶点的结合亲和力分数.

6.2 3D分子生成评价指标

3D分子生成的评价指标包括均方根偏差, 覆盖分数和匹配分数, 下面分别介绍.

(1) 均方根偏差(Root Mean Squared Deviations, RMSD): 衡量两个3D分子构象的相似性. 这里使用了RMSD的一种类型Kabsch-RMSD, 利用Kabsch算法[108]旋转和平移两个构象使其距离最小, 即将两个3D坐标矩阵叠加对齐后来计算RMSD. 公式为:

RMSD(R,ˆR)=(1nn∑i=1∥∥Ri−ˆRi∥∥2)12(4)RMSD(R,R^)=(1n∑i=1n‖Ri−R^i‖2)12(4)

RMSD越小, 两个构象差异越小; RMSD为零时两个构象相同.

(2)覆盖分数(Coverage score, COV)和匹配分数(Matching score, MAT)[93]: COV用来衡量生成构象与稳定态构象的相似性, 表示生成构象的准确性. MAT用来衡量模型生成多样性构象的能力, 这对于发现新的构象至关重要. 上述Kabsch-RMSD在计算中考虑了对称性, 因此这两个基于构象之间Kabsch-RMSD的指标COV和MAT也满足对称性. 公式为:

COV(Sg,Sr)=1|Sr|∣∣{R∈Sr∣∣RMSD(R,ˆR)<δ,ˆR∈Sg}∣∣COV(Sg,Sr)=1|Sr||{R∈Sr|RMSD(R,R^)<δ,R^∈Sg}|

(5) (5)

MAT(Sg,Sr)=1|Sr|∑R∈SrminˆR∈SgRMSD(R,ˆR)(6)MAT(Sg,Sr)=1|Sr|∑R∈Sr minR^∈SgRMSD(R,R^)(6)

其中SgSg表示生成构象, SrSr表示参考的稳定态构象, δδ是给定的RMSD阈值. 一个好的生成方法具有高COV和低MAT.

6.3 3D分子生成效果对比

由于第3节讨论的生成模型除了3.2.1以2D分子图为条件的非定向分子生成方法采用了GEOM数据集并使用了相同的实验设置, 其他模型大多没有使用统一的基准难以进行比较, 因此这里仅对3.2.1中以2D分子图为条件的非定向分子生成任务的实验结果进行总结. 实验设计沿用了Shi等人[69]的设置, 使用GEOM-QM9和GEOM-Drug基准数据集来进行3D分子构象生成, 两个数据集的训练集、验证集和测试集分别为200000、2500和22408 (GEOM-QM9)/14324 (GEOM-Drug)分子图–构象对. 采用COV和MAT对生成的构象进行质量评估, GEOM-QM9和GEOM-Drug的RMSD阈值δδ分别被设定为0.5和1.25.

实验结果均来自相关论文, 如表3所示. RDKit[109]是一款广泛应用的传统开源软件, 基于距离几何来生成构象. 其余基于深度学习的方法可以分为三类, 第一类是基于原子间距离的方法, 首先生成原子间距离或关于距离的梯度场, 然后根据距离或梯度生成原子坐标. 其中基于距离的方法有: GraphDG[50]和CGCF[93]分别通过VAE和Flow预测原子间距离矩阵, 然后通过距离几何技术预测原子坐标; ConfVAE[51]通过二层优化将上述流程改进成一个端到端的框架; EGNN[34]是一个等变模型, 将非线性变换作用在节点距离上, 并利用邻居节点坐标为标架来表示向量信息. 基于梯度场的方法有: ConfGF[69]通过学习原子间距离矩阵的梯度场来生成构象并保证等变性, DGSM[70]在ConfGF的基础上提出了一种动态图构建策略来建模长程相互作用, ClofNet[110]采用了与ConfGF相同的基于分数的框架, 通过一种具有完备局部标架的SE(3)SE(3)等变图网络来估计梯度场. 第二类是直接生成分子中所有原子的坐标: CVGAE[49]是生成原子坐标的早期尝试, 不能保持构象的旋转平移等变性; GeoMol[87]首先预测分子中原子的局部3D结构, 然后在一个复杂的过程中组装全局结构; GeoDiff[71]模拟扩散过程将分子坐标扰乱成噪声分布, 再用去噪扩散模型模拟逆向扩散过程来恢复真实构象; DMCG[94]通过设计专门的损失函数来保证生成构象的旋转平移和置换不变性, 其模型包括多个堆叠块, 每个块细化其前一个块生成的构象. 第三类是预训练模型: Zhu等人[111]提出了一种新的联合2D和3D表示的预训练方法, 设计了3个预训练任务: 掩码原子和坐标的重建、基于二维图生成三维构象、基于三维构象生成二维图. Uni-Mol[112]是一个大规模3D分子表征学习框架, 由具有SE(3)SE(3)等变transformer框架的分子预训练模型和口袋预训练模型组成, 可用于多种下游任务, 如分子性质预测、构象生成等.

从表3可以看出, 第三类的两种预训练模型在生成构象的多样性指标COV和准确性指标MAT上基本全面超越了现有的基准模型, 表明了预训练模型在构象生成任务上的有效性. 第一类基于原子间距离的方法中, 总体上基于梯度场的模型在生成效果上优于基于距离的模型, 在GEOM-Drug数据集上优势更大, 说明基于梯度场的模型能更好地建模多模态分布并生成准确而多样的构象, 在GEOM-Drug数据集上优势说明其能较好地模拟包含更多原子和官能团的复杂分布的梯度场. 第二类直接生成原子坐标的方法中, GeoDiff和DMCG在构象多样性和准确性上的表现接近或超越了预训练模型, 展现出直接生成方法的巨大潜力. GeoDiff和DMCG也超越了第一类基于原子间距离的方法中表现最优的DGSM和ClofNet, 可能是因为第一类方法中预测的原子间距离矩阵存在违反三角不等式约束的风险, 而直接预测原子坐标能避免各种距离违反情况因而在结构预测上表现出色.

7 总结与展望

基于深度学习的3D分子生成旨在利用深度学习技术设计出安全有效且新颖的3D分子构象, 在新药研发的实际应用中发挥着越来越重要的作用. 本文对3D分子结构设计进行了系统的梳理分析: 首先, 根据3D分子生成算法的输入将3D分子结构设计方法分为基于隐变量的生成、基于2D分子图的生成和基于3D分子构象的生成三大类, 接着根据算法输出的不同将其分类为定向生成和非定向生成; 总结了相关生成模型在主要公开数据集上的算法性能. 希望通过以上介绍, 读者能够对3D分子生成有一个较为全面的了解, 并在相关基于深度模型的实际任务中加以利用.

虽然基于深度学习的3D分子结构设计发展迅速并取得了很大的成功, 但由于分子结构的复杂性, 仍存在一些挑战有待解决. 我们列出了几个可能的未来发展方向, 以供参考.

(1) 更好地建模3D分子表征

在过去几年里, 生成模型取得了许多进展, 能从1D 字符串、2D分子图和3D分子几何表示中学习. 生成模型从3D表示中能学习到更丰富的信息, 其中包含更多关于分子潜在的物理信息, 但这些框架仍存在技术上的挑战. 目前建模3D分子表征主要通过EGNN或者将几何信息编码进MPNN这两种方法. 然而, 不同的几何信息对分子表征学习的影响还缺乏严格论证, 尚未建立起标准的空间信息学习方法, 这将是3D分子表示学习领域中一个很有潜力的研究方向. 此外, 大多数3D分子生成模型以自回归的方式生成3D构象, 而分子结构实际上是无序的, 这种在原子之间强加顺序的规则会引入额外偏差, 进而影响模型性能. 最后, 化合物在结构上呈现出层次结构, 如小分子由多个官能团组成, 蛋白质由氨基酸组成. 利用这些层次结构, EGNN可以尝试在多个粒度上对化合物进行建模, 以提高模型的效率和通用性.

(2) 药物设计大模型

近年来, 利用GNN进行分子表示学习取得了显著进展. 为了充分利用未标记的分子数据, 研究人员首先在大规模的分子数据库上预训练GNN以获得一个可以处理各种分子相关任务的通用编码器, 然后在下游任务中微调这些预训练GNN. 隐式编码在模型参数中的知识可以赋能多种下游任务, 有助于缓解药物设计中标记数据不足的问题. 在2D分子图上, GROVER[113]和MPG[114]等具有巨量参数的预训练GNN在学习通用分子图表示上表现出强大的表达能力. 目前, 已有一些研究将预训练模型用于3D分子表征学习, 如Zhu等人[115]提出的基于多视角对比学习的分子预训练框架MEMO, 能够自适应地集成来自2D拓扑图、3D几何、SMILES字符串和分子指纹这四种不同分子特征空间的信息, 在下游分子性质预测基准任务上性能优于SOTA, 并对学习到的模型权重提供了符合先验化学知识的解释. 近期深势科技提出的大规模3D分子表征学习框架Uni-Mol[112]包含两个相同架构的模型, 一个是由2.09亿个3D分子构象数据训练的分子预训练模型, 另一个是由320万个蛋白口袋数据训练的口袋预训练模型. Uni-Mol在14/15分子性质预测任务中表现优于SOTA, 在蛋白质–配体结合位点预测、分子构象生成等3D相关下游任务中表现尤其出色. 此外, 开发生物多模态预训练模型用于生物医药领域中的靶点发现、药物开发等环节, 将成为AI制药领域的基础设施.

(3) 干湿结合

干湿结合是指由AI计算模拟的干实验, 和实际中进行生化操作的湿实验相结合的研发方式. 首先, 由计算机在化学空间中对化合物进行模拟生成、预测和筛选. 随后, 由干实验筛选得到的化合物将提交到实验室中, 由研究人员进行合成、测试、检验等实际操作来进行优化筛选. 同时, 湿实验数据将被反馈给AI, 用于优化迭代干实验模型. 基于这种循环, 只需少量针对性的湿实验数据, AI就能够不断地生成、迭代, 直至得到具有理想活性的先导化合物. 如瑞典的阿斯利康AstraZeneca iLab对药物发现的设计–制造–测试–分析(DMTA)循环自动化进行了探索和尝试. DMTA包括设计新分子、化学合成、实验测试和分析改进四个环节, 然后再开始下一轮设计, 其中AI主要参与了DMTA循环中的设计和分析环节, 湿实验则参与了制造和测试环节. 他们提出的条件递归神经网络使化学家能够与计算机交互工作, 以加快化学空间的探索和潜在新药分子的设计[9]. 除了提升效率和节约成本, 干湿结合也是新药研发领域解决数据不足问题的一种方式.

(4) 更多真实世界的应用和数据集

虽然各种深度生成模型在理论和计算上发展迅速, 但在实际药物研发中使用AI生成模型的例子鲜有耳闻. 这是因为实际场景中仍存在一系列困难有待攻克, 比如实际分子设计中的多目标优化问题, 设计出的分子必须是可合成的, 用于筛选或估计分子属性的预测模型本身容易出错等. 另外, 目前使用的公开数据集主要关注低能的真空、水相模拟构象, 缺少实际中和靶标蛋白结合的药效构象数据集, 针对该问题除了上文提到的预训练大模型和干湿结合方法可以显著减少样本量的需求, 在算法方面还可以使用迁移学习、多任务学习、小样本学习、数据增强等方法来解决该领域数据不足的问题. 此外, 最近涌现了许多关于蛋白质折叠结构的工作, 如AlphaFold2[116]和RoseTTAFold[117]在蛋白质结构预测上取得了突破性进展, 推进了几何等变性在蛋白质等大分子结构领域的研究. 蛋白质主要是线性骨架结构, 小分子结构高度分支且常包含多种环, 这使得蛋白质折叠结构预测方法与本文讨论的小分子结构预测方法设置上存在差异. 本文讨论的方法虽然只在有限规模和复杂性的小分子数据集上得到了评估, 也可能有利于其他具有更复杂数据结构、更多样化约束和更复杂交互的应用.

实战

1.下载官方数据集与baseline

读取初始数据

import pickle

# 替换为你的 .pkl 文件路径
pkl_file_path = path to data_all.pkl'

# 加载 .pkl 文件
with open(pkl_file_path, 'rb') as file:
    data = pickle.load(file)

# 显示内容
data

[{'natoms': 23,
  'elements': ['O',
   'C',
   'N',
   'C',
   'C',
   'C',
   'C',
   'C',
   'C',
   'C',
   'Cl',
   'C',
   'C',
   'H',
   'H',
   'H',
   'H',
   'H',
   'H',
   'H',
   'H',
   'H',
   'H'],
  'coordinates': [[-1.1114834526, 1.6094932194, -1.9331164053],
   [-1.5420758704, 0.9264053946, -1.0295958406],
   [-2.7998947492, 0.3901707568, -0.9499828031],
   [-3.0798225155, -0.3450547044, 0.2701855393],
   [-1.6644553062, -0.6583011355, 0.7908986734],
   [-0.8085545958, 0.5176051925, 0.2691921396],
   [0.6697544643, 0.2462745422, 0.1225187619],
   [1.1817921402, -0.4470233574, -0.9860740901],
   [2.5424282061, -0.7382418619, -1.0923697091],
   [3.4130793875, -0.333418246, -0.076572089],
   [5.1208508658, -0.6952288215, -0.2025351142],
   [2.9327020409, 0.3627707447, 1.0345976298],
   [1.5679371316, 0.6468035778, 1.1226404623],
   [-3.5048404036, 0.6422995872, -1.6337841006],
   [-3.651866459, 0.2657274597, 0.9958601638],
   [-3.6661706337, -1.2563590389, 0.07023373],
   [-1.3017637484, -1.596493753, 0.3412719265],
   [-1.6275282242, -0.7734336107, 1.8835060808],
   [-0.9357233997, 1.3802088926, 0.949278011],
   [0.5135106936, -0.7512109932, -1.7948906059],
   [2.9327126095, -1.271115929, -1.9610179711],
   [3.6227349699, 0.6811371164, 1.8177653429],
   [1.1967768489, 1.1972849675, 1.9919902676]]},

字段中文名数据类型说明

natoms 原子数量 int 表示当前样本总共的原子数量

elements 原子元素种类 List[str] 描述当前样本的组成原子种类，第i个元素表示第i个原子的元素种类

coordinates 原子坐标 List[List[float]] 描述当前样本每个原子的3D位置，顺序与elements的value对应，第i个元素表示第i个原子在笛卡尔坐标系中的坐标

properties 化学性质 List[float] (复赛提供) 描述当前样本的绝对能量，数值经过脱敏处理。

2.训练数据

3Dmol/competition/configs/datasets_config.py

# 定义一个包含QM9数据集信息的字典，其中包含氢原子（with_h）
competition_with_h = {
    'name': 'qm9',  # 数据集名称
    'atom_encoder': {  # 原子编码器，将原子符号映射为整数
        'H': 0, 'C': 1, 'N': 2, 'O': 3, 'F': 4, 'P': 5, 'S': 6, 'Cl': 7, 'Br': 8
    },
    'atom_decoder': [  # 原子解码器，将整数映射回原子符号
        'H', 'C', 'N', 'O', 'F', 'P', 'S', 'Cl', 'Br'
    ],
    'n_nodes': {  # 不同节点数（原子数）的分子数量
        21: 2705, 25: 2464, 36: 584, 27: 2145, 9: 265, 20: 2631, 32: 1143, 14: 1519, 29: 1777, 10: 399,
        24: 2612, 22: 2742, 23: 2657, 19: 2595, 28: 2017, 30: 1560, 6: 55, 26: 2379, 35: 667, 15: 1785,
        11: 571, 31: 1333, 16: 2070, 13: 1168, 17: 2276, 12: 917, 34: 797, 37: 486, 18: 2504, 42: 167,
        39: 291, 8: 158, 40: 235, 33: 942, 4: 9, 47: 66, 38: 363, 41: 184, 51: 19, 45: 89,
        46: 92, 43: 140, 44: 99, 7: 79, 50: 26, 59: 4, 49: 32, 48: 54, 55: 4, 52: 20, 57: 7,
        53: 7, 5: 28, 56: 9, 60: 3, 1: 1, 54: 5, 3: 4, 58: 5, 2: 1
    },
    'max_n_nodes': 60,  # 数据集中分子的最大原子数
    'atom_types': {  # 不同原子类型的数量
        0: 143905848, 1: 129988623, 2: 20266722, 3: 21669359, 4: 1481844, 5: 36290, 6: 3999872, 7: 1224394, 8: 298702
    },
    'distances': [  # 分子中原子间距离的分布（可能用于分析或可视化）
        903054, 307308, 111994, 57474, 40384, 29170, 47152, 414344, 2202212, 573726,
        1490786, 2970978, 756818, 969276, 489242, 1265402, 4587994, 3187130, 2454868, 2647422,
        2098884, 2001974, 1625206, 1754172, 1620830, 1710042, 2133746, 1852492, 1415318, 1421064,
        1223156, 1322256, 1380656, 1239244, 1084358, 981076, 896904, 762008, 659298, 604676,
        523580, 437464, 413974, 352372, 291886, 271948, 231328, 188484, 160026, 136322,
        117850, 103546, 87192, 76562, 61840, 49666, 43100, 33876, 26686, 22402,
        18358, 15518, 13600, 12128, 9480, 7458, 5088, 4726, 3696, 3362, 3396, 2484,
        1988, 1490, 984, 734, 600, 456, 482, 378, 362, 168, 124, 94, 88, 52, 44, 40,
        18, 16, 8, 6, 2, 0, 0, 0, 0, 0, 0
    ],
    'colors_dic': [  # 不同原子类型的颜色编码（可能用于可视化）
        '#FFFFFF99', 'C7', 'C0', 'C3', 'C1', 'C4', 'C8', 'C9', 'C11'
    ],
    'radius_dic': [  # 不同原子类型的半径编码（可能用于可视化）
        0.3, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6
    ],
    'with_h': True  # 是否包含氢原子
}

# 定义一个函数，根据数据集名称和是否移除氢原子，返回相应的数据集信息
def get_dataset_info(dataset_name, remove_h):
    if dataset_name == 'qm9':  # 如果是QM9数据集
        if not remove_h:  # 如果不移除氢原子
            return competition_with_h  # 返回包含氢原子的QM9数据集信息
        else:
            return qm9_without_h  # 返回移除氢原子的QM9数据集信息（未定义）
    elif dataset_name == 'geom':  # 如果是GEOM数据集
        if not remove_h:  # 如果不移除氢原子
            return geom_with_h  # 返回包含氢原子的GEOM数据集信息（未定义）
        else:
            raise Exception('Missing config for %s without hydrogens' % dataset_name)  # 抛出异常
    elif dataset_name == 'qm9_second_half':  # 如果是QM9后半部分数据集
        if not remove_h:  # 如果不移除氢原子
            return qm9_second_half  # 返回包含氢原子的QM9后半部分数据集信息（未定义）

3Dmol/competition/configs/qm9_config.yaml

# 实验设置
exp_name: debug_10  # 实验名称，用于标识不同的实验
dataset: 'qm9'      # 使用的数据集，这里是QM9数据集
filter_n_atoms: null  # 筛选分子中原子数量的参数，null表示不筛选
n_report_steps: 1  # 每多少步报告一次训练进度
wandb_usr: cvignac  # Weights & Biases 用户名，用于记录实验
no_cuda: False  # 是否禁用CUDA（GPU加速）
wandb: False  # 是否使用Weighs & Biases进行实验跟踪
online: True  # 是否在线模式（True）还是离线模式（False）
data_dir: 'data'  # 数据存储的目录路径

# 训练设置
n_epochs: 100  # 训练的总轮数
batch_size: 128  # 每个批次的样本数量
lr: 0.0002  # 学习率
brute_force: False  # 是否使用暴力方法（通常用于某些特定的优化策略）
break_train_epoch: False  # 是否在训练过程中中断训练
dp: True  # 是否使用数据并行（Data Parallelism）
condition_time: True  # 是否对时间进行条件化处理
clip_grad: True  # 是否对梯度进行裁剪，防止梯度爆炸
save_model: True  # 是否保存训练好的模型
generate_epochs: 1  # 每多少轮生成一次样本
num_workers: 0  # 数据加载器使用的线程数
test_epochs: 1  # 每多少轮进行一次测试
data_augmentation: False  # 是否使用数据增强
resume: null  # 恢复训练的路径，null表示从头开始训练
start_epoch: 0  # 开始训练的轮数，通常用于恢复训练
ema_decay: 0.  # EMA（指数移动平均）衰减率，0表示不使用EMA
augment_noise: 0  # 添加的噪声量，用于增强数据

# 模型设置
model: 'egnn_dynamics'  # 使用的模型类型，这里是EGNN动力学模型
probabilistic_model: 'diffusion'  # 概率模型类型，这里是扩散模型
diffusion_steps: 200  # 扩散过程的步数
diffusion_noise_schedule: 'learned'  # 扩散噪声的时间表，learned表示学习得到的
diffusion_loss_type: 'vlb'  # 扩散模型的损失函数类型
n_layers: 6  # 模型的层数
nf: 64  # 每层的神经元数量（层大小）
ode_regularization: 0.001  # ODE（常微分方程）的正则化参数
trace: 'hutch'  # 追踪方法，hutch表示Hutchinson估计
dequantization: 'argmax_variational'  # 去量化方法
tanh: True  # 是否在坐标MLP中使用tanh激活函数
attention: True  # 是否在EGNN中使用注意力机制
x_aggregation: 'sum'  # 聚合方法，这里是求和
conditioning: []  # 条件化参数列表，可以包含多个条件
actnorm: True  # 是否使用激活归一化
norm_constant: 1  # 归一化常数，用于防止除以零

# 占位符
context_node_nf: 0  # 上下文节点的特征维度
data_root_dir: null  # 数据根目录，null表示使用默认路径

import numpy as np

# 加载.npz文件
data = np.load('path to /competition/data/competition/train.npz')

# 查看文件中包含的所有数组的名称
print(data.files)

['num_atoms', 'charges', 'positions']

['原子数', '电荷', '位置']

import numpy as np

# 加载.npz文件
data = np.load('/group_share/3Dmol/competition/data/competition/test.npz')

# 获取各个数组的数据
num_atoms = data['num_atoms']
charges = data['charges']
positions = data['positions']

# 打印数据信息
print("原子数量 (num_atoms):")
print(num_atoms)
print("\n形状:", num_atoms.shape)
print("-" * 50)

print("\n电荷 (charges):")
print(charges)
print("\n形状:", charges.shape)
print("-" * 50)

print("\n位置坐标 (positions):")
print(positions)
print("\n形状:", positions.shape)

# 关闭文件
data.close()

原子数量 (num_atoms):
[17 21 19 ... 25 32 18]

形状: (4585,)
--------------------------------------------------

电荷 (charges):
[[6 6 8 ... 0 0 0]
 [8 6 6 ... 0 0 0]
 [6 6 6 ... 0 0 0]
 ...
 [6 8 6 ... 0 0 0]
 [6 6 6 ... 0 0 0]
 [6 6 8 ... 0 0 0]]

形状: (4585, 60)
--------------------------------------------------

位置坐标 (positions):
[[[-3.21805733  0.71001875  0.25017957]
  [-2.23032949 -0.405445    0.01498888]
  [-2.50381737 -1.5580496  -0.20794194]
  ...
  [ 0.          0.          0.        ]
  [ 0.          0.          0.        ]
...
  [ 0.          0.          0.        ]
  [ 0.          0.          0.        ]]]

形状: (4585, 60, 3)

[ 0.          0.          0.        ]
[ 0.          0.          0.        ]]

形状: (4585, 60, 3)
Output is truncated. View as a scrollable element or open in a text editor. Adjust cell output settings...
输出被截断。将其作为可滚动元素查看或打开在文本编辑器中。调整单元格输出设置...

import numpy as np

# 加载.npz文件
data = np.load('/group_share/3Dmol/competition/data/competition/train.npz')

# 获取各个数组的数据
num_atoms = data['num_atoms']
charges = data['charges']
positions = data['positions']

# 打印数据信息
print("原子数量 (num_atoms):")
print(num_atoms)
print("\n原子分布:", num_atoms.shape)
print("-" * 50)

print("\n电荷 (charges):")
print(charges)
print("\n电荷分布状态:", charges.shape)
print("-" * 50)

print("\n位置坐标 (positions):")
print(positions)
print("\n坐标分布:", positions.shape)

# 关闭文件
data.close()

原子数量 (num_atoms):
[12 15 23 ... 15 38 29]

原子分布: (36683,)
--------------------------------------------------

电荷 (charges):
[[ 6  7  8 ...  0  0  0]
 [ 6  6  6 ...  0  0  0]
 [ 8  6  7 ...  0  0  0]
 ...
 [ 9  6  6 ...  0  0  0]
 [17  6  6 ...  0  0  0]
 [ 8  6  6 ...  0  0  0]]

电荷分布状态: (36683, 60)
--------------------------------------------------

位置坐标 (positions):
[[[ 2.56339586 -0.12884328 -0.26679122]
  [ 2.53272274  1.16219749 -0.06282373]
  [ 1.18434732  1.42929415  0.11404887]
  ...
  [ 0.          0.          0.        ]
  [ 0.          0.          0.        ]
...
  [ 0.          0.          0.        ]
  [ 0.          0.          0.        ]]]

 坐标分布: (36683, 60, 3)

[ 0.          0.          0.        ]
[ 0.          0.          0.        ]]

形状: (36683, 60, 3)
Output is truncated. View as a scrollable element or open in a text editor. Adjust cell output settings...
输出被截断。将其作为可滚动元素查看或打开在文本编辑器中。调整单元格输出设置...

baseline EDM 生成实战

## 配置环境
```
conda create -n com python==3.12.0
pip install -r requirements.txt
conda activate com
```
## 训练EDM模型
```
bash train.sh
```
## 采样生成数据
```
python sample.py --model_path "./outputs/competition" \
    --n_samples
```

path to competition/train.sh

python main.py --n_epochs 196 \
    # --exp_name edm_competition \
    --n_stability_samples 1000 \
    --diffusion_noise_schedule cosine \
    --diffusion_noise_precision 1e-6 \
    --diffusion_steps 500 \
    --diffusion_loss_type vlb \
    --batch_size 48 \
    --nf 256 \
    --n_layers 9 \
    --lr 1e-4 \
    --normalize_factors '[1,4,10]' \
    --test_epochs 1 \
    --ema_decay 0.9999 \
    --dataset "competition" \
    --datadir "./data/competition"\
    --test_valid_epoch 0
    --exp_name edm_competition_20250329_225950 \
    --resume outputs/edm_competition_20250329_225950

GeoDiff 尝试 GeoDiff: A Geometric Diffusion Model for Molecular Conformation Generation | OpenReview

有一知 GeoDiff 基于2D分子图的生成是指给定2D分子图来生成或优化3D分子构象;

## Environments

### Install via Conda (Recommended)

```bash
# Clone the environment
conda env create -f env.yml
# Activate the environment
conda activate geodiff
# Install PyG
conda install pytorch-geometric=1.7.2=py37_torch_1.8.0_cu102 -c rusty1s -c conda-forge
```

由于赛事要求

我们放弃使用预训练模型与增加数据集的策略，转而从算法层面上研究使模型拥有更好的输出能力与泛化专业性。

听从骆佬的建议，比较离谱的提分思路 😀
我们找到最新的图神经网络（GotenNet）找到它，然后把相关的特征数据输进去，有两条同样的特征数据输进去，让它做一个自回归。有点像对比学习一样的叫什么 simple1样的字回归这回归之后，他能得到一些比较好。跑的这个映射，拿到这个映射之后再去接一个模型，然后去做一些扩散，就像现在大家看到这样做一些扩散diffusion，也可能会有一些好的分手。
GotenNet++自回归AR

Mucunshuo/GotenNet: [ICLR 2025] GotenNet: Rethinking Efficient 3D Equivariant Graph Neural Networks

核心三个要素图神经网络。这个已经是被广泛证实说在材料领域很好用的东西，因为材料本身是一种图结构，另一个事情是。等变性，刚刚也有同学提到等变性，因为物理上面有等变性，相当于是引入了一些物理先验知识进去。第三个事情是你可以用DFT （密度泛函理论）去做验算，因为 DFT 是一种稳定的解释，解析的找到这个材料是否可靠的一种方法