一种基于扩散模型的多目标分子生成方法和系统

本发明属于生物医药，更具体地，涉及一种基于扩散模型的多目标分子生成方法和系统。

背景技术：

1、药物从研发到上市需经过药物发现、预临床研究、临床试验等多个阶段。在药物发现阶段，科研人员通过实验筛选或计算机虚拟筛选从大量化合物中确定潜在的药物候选分子。这些候选分子还需经过进一步的药效实验验证。传统依赖实验筛选的随机高通量方法存在效率低的问题，这类方法需要研发人员逐一合成大量候选化合物，进行药效实验以找到潜在药物分子，整个筛选流程周期长、成本高。计算机辅助的虚拟筛选如分子动力学模拟也可用于药物发现。这类方法可以评估更多化合物，不需要逐一合成，但是对每个化合物进行精确的量子化学计算要求计算资源巨大，也难以处理复杂生物大分子。近年来，基于深度学习的分子生成方法提供了一种更高效的筛选思路。这类方法使用神经网络对已有药物数据进行训练，学习编码药物分子的化学结构特征。然后根据网络对药效的预测，有选择性地生成新分子。

2、目前基于深度学习的分子生成研究有两种方法，一种是基于变分自编码器或生成式对抗网络的模型，其直接建模数据分布，这类方法通过编码分子结构的语义特征，并加入分子性质标签，学习数据联合分布。然后可以根据采样或搜索这个学习到的联合分布，生成新的分子；另一类方法是基于语言模型的模型(例如transformer模型)，这类方法利用大规模分子数据进行预训练，学习编码分子smiles表示的语法及化学语义特征。然后根据这种语言模型采样生成新的分子smiles字符串。

3、然而，上述两种方法均存在一些不可忽略的缺陷：

4、第一，基于生成对抗网络的模型虽然能够生成看似高质量的样本，但是其难以训练，而且训练容易坍缩，从而导致生成样本唯一性很差；

5、第二、基于变分自编码器的模型建模复杂数据分布时，可能会过度平滑样本分布，导致生成样本多样性不足，这会导致生成样本的多样性不高；

6、第三，语言模型容易出现“数据饥饿”现象，也即需要更多的数据来训练模型，而在药物研发领域，高质量的小分子数据集样本数量很少，语言模型会对小样本数据集过拟合严重，进而导致生成的样本新颖性很低；

7、第四，现有基于变分自编码器或语言模型更注重单一药效或性质的优化。而高质量的药物候选物需要同时兼顾多个性能指标(如药效、毒性、合成难易度等)，因此该模型在处理这类多目标分子生成问题上还存在明显不足。

技术实现思路

1、针对现有技术的以上缺陷或改进需求，本发明提供了一种基于扩散模型的多目标分子生成方法和系统，其目的在于，解决现有基于生成对抗网络的模型虽然能够生成看似高质量的样本，但是其难以训练，而且训练容易坍缩，从而导致生成样本唯一性很差的技术问题，以及现有基于变分自编码器的模型建模复杂数据分布时，可能会过度平滑样本分布，导致生成样本多样性不足、生成样本的多样性不高的技术问题，以及现有语言模型需要更多的数据来训练模型，而在药物研发领域，高质量的小分子数据集样本数量很少，语言模型会对小样本数据集过拟合严重，进而导致生成的样本新颖性很低的技术问题，以及现有基于变分自编码器或语言模型更注重单一药效或性质的优化。而高质量的药物候选物需要同时兼顾多个性能指标，导致该模型在处理这类多目标分子生成问题时存在明显不足的技术问题c。

2、为实现上述目的，按照本发明的一个方面，提供了一种基于扩散模型的多目标分子生成方法，包括以下步骤：

3、(1)获取带有标签的活性分子数据集，将活性分子数据集输入到预先建立好的子结构搜索模型，以获取活性分子数据集的语义子结构s。

4、(2)将步骤(1)得到的语义子结构转换成含噪语义子结构简化分子线性输入规范smiles序列ynoise，并将该含噪语义子结构smiles序列ynoise进一步转换为含噪语义子结构词元y0。

5、(3)将步骤(2)获取的含噪语义子结构序列y0输入预先训练好的扩散模型，以获取多目标分子。

6、优选地，步骤(2)具体为，首先使用分子处理软件rdkit将步骤(2)获取的图格式的语义子结构s转换成语义子结构smiles序列y；随后，在语义子结构序列y的活性位点位置附近添加随机长度的随机原子，以获得含噪语义子结构smiles序列ynoise；最后，使用扩散模型中预先建立的词表将含噪语义子结构smiles序列ynoise进行编码，以获取独热编码形式的含噪语义子结构序列y0。

7、优选地，步骤(3)包括以下子步骤：

8、(3-1)设置时刻t＝t，从扩散模型预先建立的词表随机采样获取t时刻长度为l的含噪序列xt，其中t表示总采样步数，且有t＝100，l表示序列长度。

9、(3-2)判断是否有t＜1，如果是则转入步骤(3-6)，否则将含噪序列xt和时刻t输入到扩散模型中，以获取输出然后进入步骤(3-3)。

10、(3-3)将步骤(3-2)获取的输出输入扩散模型，以获取t-1时刻的含噪序列xt-1，然后进入步骤(3-4)；

11、(3-4)判断是否有t＜γ，如果是则进入步骤(3-5)，其中γ表示截断时间，且有γ＝10，否则获取步骤(2)获取的含噪语义子结构序列y0在t-1时刻的特征yt-1，将步骤(3-3)获取的t-1时刻的含噪序列xt-1与特征yt-1通过线性插值进行特征融合，以获取融合后的t-1时刻的含噪序列xt-1，然后进入步骤(3-5)；

12、(3-5)设置时刻t＝t-1，xt＝xt-1，并返回步骤(3-2)；

13、(3-6)获取t＝0时刻的序列x0，并进入步骤(3-7)；

14、(3-7)通过扩散模型预先建立的词表将步骤(3-6)获取的序列x0解码成smiles序列，以获取多目标分子。

15、优选地，子结构搜索模型是通过以下步骤训练得到的：

16、(4-1)获取gsk3活性分子数据集和jnk3活性分子数据集作为子结构搜索模型的输入，使用该gsk3活性分子数据集训练活性打分器，并利用训练好的活性打分器和蒙特卡洛树搜索mcts算法获取gsk3活性分子数据集中的活性打分器sg和语义子结构集合ssg；

17、(4-2)使用步骤(4-1)获取的jnk3活性分子数据集训练活性打分器，并利用训练好的活性打分器和蒙特卡洛树搜索mcts算法获取jnk3活性分子数据集的活性打分器sj以及语义子结构集合ssj；

18、(4-3)将步骤(4-1)得到的gsk3活性分子数据集中的语义子结构集合ssg和步骤(4-2)得到的jnk3活性分子数据集中的语义子结构集合ssj融合成一个多目标融合的语义子结构集合，作为子结构搜索模型的输出。

19、优选地，步骤(4-1)具体包括以下步骤：

20、(4-1-1)获取gsk3活性分子数据集，并在预先建立的预训练分子数据集中随机获取数量为该gsk3活性分子数据集100倍的多个非活性分子构成非gsk3活性分子数据集；

21、(4-1-2)将步骤(4-1-1)获取的gsk3活性分子数据集和非gsk3活性分子数据集合并成一个数据集，并对该数据集利用k-means算法进行聚类，以获取该合并数据集的聚类结果；

22、(4-1-3)根据步骤(4-1-2)获取的聚类结果，在步骤(4-1-2)获取的合并数据集中的gsk3活性分子数据集的聚类中心附近采样多个非活性分子，作为新的非gsk3活性分子数据集；

23、(4-1-4)为步骤(4-1-3)获取的所有新的非gsk3活性分子数据集中的所有非活性分子在[0，0.5]之间随机赋予标签，并为步骤(4-1-1)获取的gsk3活性分子数据集的中的所有活性分子在[0.7，1.0]之间随机赋予标签，所有被赋予标签的活性分子和非活性分子构成活性打分器数据集；

24、(4-1-5)利用分子处理软件rdkit将步骤(4-1-4)获取的活性打分器数据集转换成分子指纹数据集，并将分子指纹数据集按照8∶2的比例随机划分为训练集和测试集；

25、(4-1-6)使用步骤(4-1-5)获取的训练集和测试集对随机森林模型进行训练，以获取训练好的随机森林模型作为gsk3活性数据集的活性打分器sg，其中随机森林模型的输入是训练集和测试集的分子指纹，输出为预测的分子活性标签；

26、(4-1-7)使用(4-1-6)获取的活性打分器sg和mcts算法获取gsk3活性分子数据集中的语义子结构集合ssg，下标g用于标识当前活性分子数据集。

27、优选地，步骤(4-1-7)包括以下子步骤：

28、(4-1-7-1)设置计数器k＝1，初始化语义子图集合ssg为空集合，并转入步骤(4-1-7-2)；

29、(4-1-7-2)图g初始为gsk3活性数据集中第k个活性分子，设置计数器j＝1，子图g＝g，设置mcts算法的四个关键状态量，并转入步骤(4-1-7-3)：

30、c(g，a)＝0为从图g到子图g过程中，断开边操作的边数量总和；

31、v(g，a)＝0为从图g到子图g过程中，活性打分器的得分总和；

32、q(g，a)＝v(g，a)/c(g，a)＝0为从图g到子图g过程中，活性打分器的平均得分；

33、r(g，a)＝sg(g′)为从子图g执行断开边操作a得到的子图g′，利用活性打分器sg对子图g′进行打分，获得子图g′得分为sg(g′)；

34、(4-1-7-3)判断是否有j＞n，或者子图g中的节点数小于mmin，如果是，则转入步骤(4-1-7-10)，否则设置计数器i＝1，并转入步骤(4-1-7-4)，其中n为最大递归迭代次数，优选为n＝20，mmin为语义子结构最小节点数，优选为mmin＝15；

35、(4-1-7-4)判断是否有i＞第j次迭代时对应的子图g中的非环边总条数d，如果是，则转入步骤(4-1-7-7)，否则，对子图g的第i条边执行断开边操作ai，获得新的子图gi，并使用步骤(4-1-6)获得的活性打分器sg计算新的子图gi的得分sg(gi)，并转入步骤(4-1-7-5)；

36、(4-1-7-5)判断是否有sg(gi)≥δ，且新的子图gi中的节点总数小于mmax并大于mmin，如果是，则将新的子图gi加入到语义子结构集合ssg中，并转入步骤(4-1-7-6)，否则过程结束，其中δ是一个0到1之间的阈值，优选为δ＝0.5，mmax为语义子结构最大节点数，优选为mmax＝25；

37、(4-1-7-6)设置计数器i＝i+1，并返回步骤(4-1-7-4)；

38、(4-1-7-7)对子图g的所有d条非环边，执行断开边操作a1，...，ad，以获得多个新的子图g1，...，gd，并根据以下公式选择断开边操作an，以得到新的子图gn，然后转入步骤(4-1-7-8)：

39、an＝argmax q(g，a)+w(gi，ai)

40、

41、其中∑bc(gi，b)是新的子图gi中含有的非环边的总条数，b表示对子图gi的断开边操作，为了与断开子图g的断开边操作a区分。λ是用于平衡博弈论中探索与利用的参数，优选地λ＝10，这样就能够每次选择一个更好的断开边操作an以及子图gn，argmax q(g，a)+w(gi，ai)表示从所有可能的断开边操作ai中选择一个使得q(g，a)+w(gi，ai)值最大的断开边操作；

42、(4-1-7-8)利用步骤(4-1-7-7)获得的断开边操作an以及新的子图gn，对步骤(4-1-7-2)中设置的四个关键状态量进行更新，然后转入步骤(4-1-7-9)：

43、c(g，a)＝c(g，a)+1

44、v(g，a)＝v(g，a)+s(gn)

45、q(g，a)＝v(g，a)/c(g，a)

46、r(g，a)＝sg(gn)

47、(4-1-7-9)设置j＝j+1，g＝gn，a＝an，并返回步骤(4-1-7-3)；

48、(4-1-7-10)设置k＝k+1，并返回(4-1-7-2)。

49、优选地，步骤(4-3)包括以下子步骤：

50、(4-3-1)设置计数器x＝1，c＝1，sx为语义子结构集合ssg中的第x个子结构，sc语义子结构集合ssj中的第c个子结构，多目标融合的语义子结构集合ss＝{}，并转入步骤(4-3-2)；

51、(4-3-2)判断是否有x＞语义子结构集合ssg的大小kg，如果是则过程结束，否则转入步骤(4-3-3)；

52、(4-3-3)判断是否有c＞语义子结构集合ssj的大小lj，如果是，则转入步骤(4-3-6)，否则转入步骤(4-3-4)；

53、(4-3-4)使用分子处理软件rdkit获取语义子结构集合ssg中的第x个子结构sx和语义子结构集合ssj中的第c个子结构sc二者的公共子结构a，并通过公共子结构a将si与sj进行合并，以获得合并子结构s，判断是否有sg(s)≥δ且sj(s)≥δ，如果是则把合并子结构s加入到多目标融合的语义子结构集合ss中，并转入步骤(4-3-5)，否则过程结束；

54、(4-3-5)设置x＝x+1，并返回步骤(4-3-3)；

55、(4-3-6)设置c＝c+1，并返回步骤(4-3-2)。

56、优选地，扩散模型是通过以下步骤建立的：

57、(5-1)从chembl数据库获取一千万条无标签的smiles序列作为预训练数据。

58、(5-2)利用分子处理软件rdkit对步骤(5-1-1)获取的预训练数据进行数据处理操作，以获取预训练分子序列数据集；

59、(5-3)使用正则表达式对步骤(5-2)获得的预训练分子序列数据集进行分词，以获取多个词元，所有词元组成词表。

60、(5-4)将步骤(5-2)获取的预训练分子序列数据集按照8∶2的比例划分成训练集和测试集；

61、(5-5)为步骤(5-4)获取的训练集中的每一条序列seq随机采样一个时刻t，其中t是1到t的整数，与步骤(3-1)相同，t＝100；本步骤到步骤(5-16)的流程如图4所示；

62、(5-6)通过步骤(5-3)获取的词表将步骤(5-5)获取的每一条序列seq编码成长度为l的词元序列(如无特别指明，以下序列均表示词元序列)，其中l是1到512的整数；

63、(5-7)使用独热编码方式将步骤(5-6)获取的词元序列编码成序列x0，其中，序列x0中元素的取值为0或者1，序列的维度为l×k，k为步骤(5-3)获取的词表的大小；

64、(5-8)根据步骤(5-7)获取的序列x0获取t时刻的含噪序列xt：

65、(5-9)将步骤(5-8)获取的含噪序列xt输入嵌入层神经网络，以获取特征hx，其中嵌入层神经网络的嵌入维度为d；

66、(5-10)使用transformer的正余弦编码方法对步骤(5-5)获取的时刻t进行编码，以获取特征ht；

67、(5-11)将步骤(5-9)获取的特征hx和步骤(5-10)获取的特征hx相加，以获取特征hx；

68、(5-12)使用12层transformer编码器对步骤(5-11)获取的特征hx进行编码，以获取特征hout；

69、(5-13)使用矩阵wo将步骤(5-12)获取的特征hout编码成特征hfinal，其中矩阵wo的维度为d×k，：

70、

71、(5-14)使用归一化指数函数softmax将步骤(5-13)获取的特征hfinal转换为扩散模型预测t时刻含噪序列xt在t＝0时刻的概率分布：

72、(5-15)根据步骤(5-8)获取的t时刻的含噪序列xt获取t-1时刻的含噪序列xt-1，并根据步骤(5-14)获取的扩散模型预测的t时刻含噪序列xt在t＝0时刻的概率分布获取扩散模型预测的t时刻含噪序列xt在t-1时刻的概率分布

73、(5-16)根据步骤(5-15)获得的扩散模型预测的t时刻含噪序列xt在t-1时刻的概率分布和t-1时刻的含噪序列xt-1计算损失函数其中dkl表示kullback-leibler散度；

74、(5-17)重复上述步骤(5-5)到步骤(5-16)，直到损失函数的值最小为止，从而得到训练好的扩散模型。

75、优选地，步骤(5-2)具体为，首先，对预训练数据进行去重处理，以得到数据集d1；然后，去除数据集d1中的立体异构信息，以得到数据集d2；最后，去除数据集d2中序列长度超过512的所有序列，最终得到预训练分子序列数据集；

76、步骤(5-8)是采用以下公式：

77、xt＝q(xt|x0)＝c(xt；p＝x0q1q2…qt)

78、其中q表示序列x0的数据分布，q(xt|x0)是计算含噪序列xt关于x0的条件概率，p＝q1q2…qt表示序列x0的所有词元转变成其他词元的概率，c(xt；p)表示将含噪序列xt以概率p转变为均匀类别分布，qt是t时刻的概率转移矩阵，其取值在0到1之间，矩阵的维度为k×k，qt的具体形式为：

79、

80、其中k表示步骤(5-3)获取的词元的类型总数，i和j分别表示qt的行号和列号，βt为常量，其取值在0到1之间，随着t增大单调递增。

81、优选地，步骤(5-15)是采用以下公式：

82、xt-1＝q(xt-1|xt)＝c(xt-1；p＝xtqt-1)

83、

84、其中qt-1为t-1时刻的概率转移矩阵；

85、步骤(5-16)中损失函数的计算公式为：

86、

87、一种基于扩散模型的多目标分子生成系统，包括：

88、第一模块，用于获取带有标签的活性分子数据集，将活性分子数据集输入到预先建立好的子结构搜索模型，以获取活性分子数据集的语义子结构s。

89、第二模块，用于将第一模块得到的语义子结构转换成含噪语义子结构简化分子线性输入规范smiles序列ynoise，并将该含噪语义子结构smiles序列ynoise进一步转换为含噪语义子结构词元y0。

90、第三模块，用于将第二模块获取的含噪语义子结构序列y0输入预先训练好的扩散模型，以获取多目标分子。

91、总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：

92、(1)本发明由于采用了步骤(5-1)到(5-17)，其基于扩散模型的加噪去噪的数学框架，没有对抗式训练机制，使得训练容易且比较稳定不容易产生模式坍缩现象，能够很好的建模分子序列数据分布，生成的分子唯一性很高；

93、(2)本发明由于采用了步骤(5-7)，对样本加噪时采取渐进式的提升噪音比例(βt)，这可以让扩散模型学会更加复杂的数据分布，能够让生成的分子多样性更高；

94、(3)本发明使用语义子结构来引导扩散模型的生成，避免了直接将数据标签输入模型，只需通过调控语义子结构的引导系数控制条件生成过程，能够极大避免过拟合问题，提升生成样本新颖性；

95、(4)本发明基于蒙特卡洛树搜索算法搜索语义子结构，然后将多个语义子结构融合成一个含有多语义的子结构，通过多语义的子结构来引导扩散模型的生成，同时兼顾了多个性能指标；

96、(5)本发明使用语义子结构来引导扩散模型的生成也避免了在高质量标签数据较少情况下使用强化学习常导致过拟合的现象，能够更高效的利用小数据集。