1.本发明属于药物筛选技术领域，特别涉及一种结合机器学习和构象计算提高药-靶活性预测精度的方法。
背景技术：2.在药物研究过程中，巨额的人工费用和漫长研发周期始终掣肘新型药物的问世，直到使用计算机的虚拟筛选技术应用于药物化学领域，一定程度上加快了前期药物研发的进程。传统的分子对接是基于计算的方法，通过靶点特征以及药物分子和靶点之间的相互作用来预测其结合模式和亲和力，例如使用rosetta、ledock、autodockvina等分子对接软件。这种对接方法目前依然广泛应用于虚拟筛选领域。但此类方法受限于需要精准的靶点蛋白结构，众多蛋白质的结构被未被解析，而且此类方法计算效率较低，通常依赖于高性能的计算资源。近年来的研究表明，将特征学习能力集成到机器学习模型中可以提高预测性能，使用人工智能模型针对目标蛋白筛选潜在活性药物也逐步被应用于此领域。预测分子活性在指导药物发现方面起着决定性作用。
3.在生物医学领域，数据集的规模正逐步扩大。近年来，众多学术机构针对数据集规模问题做出了极大努力，研究报道了众多具备实验活性的药物-靶复合物数据。公开可用的生物活性数据的数量每年都在不断增加，并已成为众多生命科学研究的宝贵资源。同时，在过去的几十年里，由于x射线晶体学、核磁共振光谱和电子显微镜等技术在精度和吞吐量上显著提高，研究者也通过实验手段发表了大量的活性结构数据。
4.命名实体识别方法作为自然语言处理领域的初始步骤，已经有相当的研究高度。与通用领域相比，新报道的活性化合物为生物医学领域专有实体，具有领域性强、特征性复杂的特点，将通用领域的命名实体识别模型应用于该领域会出现将专有实体切分的情况；同时，生物医学领域对于化合物实体的命名没有统一的规则，存在大量含特殊字符、模糊缩写的情况，导致通用领域的命名实体识模型无法识别该领域的专有名词；且传统的命名实体识别模型大多针对于句子实现，当应用于文献级别的实体识别时，会出现序列标注不一致的情况。
5.随着机器学习和深度学习技术的逐步成熟，药物发现领域的机器学习模型也逐步迭代。在初始阶段，有研究者从简单的化学子结构和序列信息中自动提取药物和靶点的特征，然后使用深度置信网络(dbn)构建分类模型；在之前的基础上，有研究将药物和靶点的分子指纹或分子描述符作为输入，不考虑原子之间局部联系和氨基酸局部化学结构，构建了半监督的深度学习模型。这些方法都是通过药物靶标对的特征来预测药物靶标对是否相互作用，使用药物和蛋白质序列的一维表示，并使用强大的高级深度学习模型来提取复杂的局部化学信息和序列中局部结构之间的上下文关系，最后将药物与靶标的特征信息进行拼接，输入神经网络进行预测。但这些方法并没有像对接方法一样充分应用到靶蛋白的结构信息，很多药物分子呈现活性与靶蛋白的结构密切相关。
6.另一类被研究者广泛应用的深度学习模型的是图神经网络(gnn)。基于图神经网
络的方法首先将生物分子的原子及原子间的相互关系分别抽象成图数据节点和边，构建一个包含药物和靶点的网络，然后学习分子图特征，根据已知节点和边预测未知边，其主要思想是药物倾向于与相似的目标结合，反之亦然。在过去的几年中，图神经网络的概念、操作和模型一直在不断地演进和发展，并且图卷积神经网络在提取嵌入图形方面表现出很强的性能，极大地促进了药物-靶点相互作用预测的发展。尽管图方法应用广泛，但这类方法主要是表征固有异构网络的节点，难以实现异构网络外的药物目标预测。此外，节点的特征表示主要基于异构网络中的拓扑信息，没有深入考虑药物和靶点的生物结构信息，其准确性仍存在较大优化空间。有研究表明，在预测精度和泛化性上，很多传统的机器学习算法并不低于一些图方法，甚至在某些任务上表现还为优越。
7.基于深度学习的模型依赖于大量高质量的活性数据集，即使数据集规模有所扩大，但现有的公共数据集仍然缺少足够的活性分子和相关靶点。正样本缺失、数据集质量不高使得模型很难从中挖掘到深层次的关联关系和作用模式，以致于模型缺少足够的泛化性。另外，对靶点结构信息的忽略也同样限制了深度学习方法进一步突破性能瓶颈，从而在药物研发领域实际应用。
8.该领域的重点仍然是寻找更先进、性能更好和更合理的方法。这类方法的特点是基于生物学特性，能够充分利用靶蛋白和分子结合信息，可扩展性比较强。在模型泛化能力强的情况下，可以预测数据集外的药物靶点相互作用，取得更好的预测精度，并能兼顾计算效率。
技术实现要素：9.本发明提供一种结合机器学习和构象计算提高药-靶活性预测精度的方法，以提高药-靶结合活性预测的精度，使得药物筛选更高效、成本更低，且结果的准确度和可靠性提高。
10.本发明是通过如下技术方案来实现的
11.一种结合机器学习和构象计算提高药-靶活性预测精度的方法，所述方法包括命名实体识别获取文献中活性数据集，训练神经网络以输出权重参数，样本聚类以及构建多分类器；具体步骤如下：
12.(一)命名实体识别获取文献中活性数据以获得活性数据集，包括以下步骤：
13.1)从现有论文数据库中获取活性化合物和靶点的文献信息；
14.2)使用pubmedbert预训练模型将输入文献信息向量化；
15.3)使用bilstm获取包含上下文信息的状态向量；
16.4)利用attention机制实现对不同词汇的不同关注；
17.5)使用crf层获取最优标签序列；
18.6)根据序列标注结果，将标注为化合物、靶点及其详细信息的实体抽取保存；
19.(二)构建能量参数生成模型，优化模型输出权重参数，形成多个准确的能量函数模型
20.包括以下步骤：
21.1)从(一)中获得的活性数据集抽取一部分作为小样本集，将小样本集使用加入参数生成器的神经机器翻译模型(nmt)进行训练，生成一组权重系数，通过反复迭代训练此参
数生成神经网络，输出最终的权重系数组合，所述的参数生成神经网络是指所述加入参数生成器的神经机器翻译模型，应用此参数生成神经网络输出的权重系数组合与rosetta原有的能量项计算公式组成新的优化后的能量函数模型进行活性预测，公式如下：
22.δe
total
＝∑ω
iei
(θi，aai)
23.能量函数近似于生物分子和靶蛋白耦合构象的能量，这个能量称为
△etotal
，根据能量项的线性组合计算得出ei，ei为几何自由度θi、化学恒等式aai的函数，并按每个项ω的权重缩放；
24.2)将(一)中获得的完整的活性数据集输入步骤1)所得的能量函数模型，将输出结果与文献报导的实际活性值进行比较，设置活性阈值，在能量函数模型中验证，将判定为活性的样本取出，把高于活性阈值而被判定为活性的样本按照8：1：1的比例划分为训练集、测试集和验证集，形成新的样本集，使用训练集再次对参数生成神经网络进行迭代训练，输出新的权重系数组合，从而形成新的能量函数模型；经过迭代训练后的新能量函数模型使用测试集用以测试能量函数模型的泛化误差，验证集用以评估模型准确度以及设定相关超参数；
25.3)重复步骤2)，进行多次迭代，不断优化参数生成神经网络，验证准确度时设置阈值，达到阈值时终止步骤2)，最终形成一个相对较为准确的能量函数模型；
26.重复步骤1)、2)、3)，从完整样本集中随机选取小样本集随时抽取多次，获得多个对具有某些特征有偏好识别能力、更敏感的能量函数模型，所获得的能量函数模型数量与随机抽取样本集的次数一致；这些模型对迭代训练过程中其最后一次训练所输入的数据集具有相对较好的活性预测准确度，在此依据活性表现将最后一次训练所输入的数据完成样本聚类；
27.(三)构建多分类器
28.构建多分类器以区分新化合物分子，以匹配(二)中获得的相对应的能量函数模型预测活性，具体构建方法如下：
29.1)使用自动编码器进行分子、靶点编码
30.2)使用boosting方法和softmax方法生成基分类器
31.3)选择、组合分类器。
32.进一步的，所述(一)中步骤2)的具体操作是：采用基于pubmedbert的预训练模型将包含m个语句的输入文献d＝(x1，...，x
t
，...，xm)向量化，得到m个包含上下文语义信息的向量x＝(x1，...，x
t
，...，xn)，其中n为语句中单词的数量，xi为语句中包含第i个词上下文信息的词向量；
33.进一步的，所述(一)中步骤3)的具体操作是：以语句向量x作为bilstm模型的输入，提取每个词向量xi的上下文特征得到输入的语句向量对应的隐藏层状态向量h＝(h1，...，h
t
，...，hn)，其中h
t
为t时刻隐藏层状态向量；
34.进一步的，所述(一)中步骤4)的具体操作是：以语句向量x与各时刻隐藏层状态向量h作为attention层输入，针对目标词与其余词汇之间的相似度对每个词汇给予不同的程度的注意，最终得到融合各词汇信息和隐层状态的向量z＝(z1，...，z
t
，...，zn)；具体的，首先利用余弦距离计算目标词x
t
与文献中其余词汇xj之间的相似度分数score(x
t
,xj)，再使用softmax函数计算得出相应的归一化注意力权重参数α
t，j
，根据该权重参数与隐层向量h
计算得出文献级别的全局向量g＝(g1，...，g
t
，...，gn)，其中g
t
表示目标词x
t
对应的全局向量，最后再将全局向量g与对应的隐层向量h拼接，使用tanh函数计算得出attention层的输出向量z；
35.进一步的，所述(一)中步骤5)的具体操作是：以向量z作为tanh层的输入计算得出各词向量可能的标签对应的分数，将神经网络输出的分数记为p，p的维度为k x n，k为每个词可能的标签数量，n为词的个数；最后的crf层综合考虑本层转移矩阵t和矩阵p，计算出最优的标签序列作为标注结果；
36.进一步，所述(三)中步骤1)的分子、靶点编码：在(一)获取文献中的化合物的活性信息时，可同时获得分子的标准格式文件；分子文件由开源工具转换成smiles文件，其中包含了结构信息，靶点文件可开源获得pdb或者其他格式的构象文件，使用autoencoder的前馈神经网络将二者编码；
37.进一步，所述(三)中步骤2)的具体操作是：
38.使用相同的分类模型，而通过对训练集中的包括实例、属性和类标的变化来形成不同的训练输入，以最终形成不同的分类器；当给出原始的训练集合后，通过包括删除、添加、抽取在内的方法来操作训练集合中的部分实例来形成新的训练集合；使用boosting方法来实现产生基分类器。
39.进一步的，所述(三)中步骤3)的具体操作是：
40.首先解决如何决定测试数据所属的特征空间；得到包括原子类型、分子的结构、化学键在内的数据属性，利用属性值的统计信息，用数据属性对评价数据集进行划分；
41.然后，在这些划分上测定每个基分类器的分类正确率，正确率最高的被选出来对测试数据进行分类；
42.进一步的，所述(三)中步骤4)的具体操作是：
43.在确定采用哪种组合方法之前，首先要确定基分类器的分类结果的表示形式，对每一个类标y
x
，首先将所有基分类器给出的预测类标可能是它的概率相加并求平均值，然后这些值最大的类标作为最终的分类结果；分类公式如下：
[0044][0045]
本发明与现有技术相比的有益效果：
[0046]
命名实体识别获取文献中活性数据集的方法在文本向量化时采用生物医学nlp领域的特定预训练模型pubmedbert，可涵盖大部分生物医学术语，避免将一些专有名词分解为无意义子词，可实现更加准确的语义编码，采用bilstm可利用两个方向的句子信息，提高标签标注的准确率，加入attention层关注目标词汇与全文其它词汇的相似度减少标签不一致的情况，最后利用crf层对得到的标签序列进行修正，得到更加合适的标签序列。
[0047]
构建能量参数训练网络，训练神经网络以输出权重参数的方法对包括但不限于以上能量项，能够对全原子进行能量参数优化。生物构象所包含的原子之间相互作用的能量项可以分为范德华力、静电能、溶剂、氢键、二硫键合等多项，根据原子间的相互作用，建立明确的建模特征，训练出结合构象中主要起作用的项并优化权重，使得建立更加有针对性的准确度更高的能量函数成为可能。通过对能量函数进行参数优化的迭代训练，对输入样本进行能量项参数预测计算，与目前多能量函数对接模块相比，这种模型具有更高的准确性和针对性。
附图说明
[0048]
图1命名实体识别建模方法架构图；
[0049]
图2权重参数优化模型流程图；
[0050]
图3分类系统流程图；
[0051]
图4 autoencoder网络架构图；
[0052]
图5 boosting生成分类器结构图；
[0053]
图6本发明方法的整体框架图；
[0054]
图7 hiv药物分子结构式图；
[0055]
图8 bace1抑制剂分子结构图；
[0056]
图9为凝血酶抑制剂分子结构图；
[0057]
图10为几丁质分子结构式图；
[0058]
图11为糖苷键分子图。
具体实施方式
[0059]
下面通过实施例结合附图来对本发明的技术方案做进一步解释，但本发明的保护范围不受实施例任何形式上的限制。
[0060]
实施例1
[0061]
(一)命名实体识别获取文献中活性数据集
[0062]
本发明采用的技术方案是首先从文献数据库中以命名实体识别的方法获取已报导的活性数据，这些数据来源广泛、规模巨大，且通过实验具有准确的实测活性数据，是一个高质量的多源异构数据集，以此作为机器学习模型的数据集；针对将现有通用领域的命名实体识别模型应用于化合物实体识别时的不足，本实施例采用以下方法进行对化合物进行实体识别，以生物医学领域文献化合物实体识别为例，包括以下步骤：
[0063]
步骤1：从现有论文数据库中获取活性化合物、靶点的文献信息；将成对的句子作为pubmedbert的输入，判断2个句子的前后顺序。从结构上来说，pubmedbert是将多个transformer编码器堆叠在一起进行特征提取，每个transformer编码器由self-attention层和前馈神经网络层组成。self-attention是transformer的核心机制。
[0064][0065]
使用self-attention机制的意义在于，它不仅将单词对于全文的重要性进行了编码，还摒弃了传统的循环神经网络结构，在解决了传统模型长期依赖问题的同时，大大加快了模型的并行计算能力。
[0066]
步骤2：使用pubmedbert预训练模型将输入文献信息向量化；采用基于pubmedbert的预训练模型将包含m个语句的输入文献d＝(x1，...，x
t
，...，xm)向量化，得到m个包含上下文语义信息的向量x＝(x1，...，x
t
，...，xn)，其中n为语句中单词的数量，xi为语句中包含第i个词上下文信息的词向量。
[0067]
获取pubmedbert层的输出后，将其与pos embedding，chunking进行拼接，加入到前向和后向lstm中进行编码，并根据上下文，使用softmax函数给出当前单词对应标签的概率。lstm层的主要结构：
[0068]it
＝σ(x
t
·
wi+h
t-1
·
wi′
+bi)
[0069]ft
＝σ(x
t
·
wf+h
t-1
·
wf′
+bf)
[0070]ot
＝σ(x
t
·
wo+h
t-1
·
wo′
+bo)
[0071][0072][0073][0074]
这里，σ是sigmoid函数，i、f、o和c分别表示输入门、遗忘门、输出门和记忆单元；是点积运算，ω和b代表输入门、遗忘门和输出门的权重矩阵和偏置向量；x
t
指的是t时刻的网络输入，同时对应于词表示层的输出。
[0075]
步骤3：使用bilstm获取包含上下文信息的状态向量；
[0076]
以语句向量x作为bilstm模型的输入，提取每个词向量xi的上下文特征得到输入的语句向量对应的隐藏层状态向量h＝(h1，...，h
t
，...，hn)，其中h
t
为t时刻隐藏层状态向量。
[0077]
bilstm、网络层的输出作为attention层的输入。在bilstm层之上使用一个新的关注层在文档级捕获类似的单词注意力。在关注层，我们引入关注度矩阵a来计算当前目标词与文档中所有词之间的相似度。关注权重值，即关注度矩阵，是通过比较当前目标词表示x
t
和文档中的j个词表示xj而得出的
[0078][0079]
这里，分数被称为对齐函数，对于该函数使用了余弦距离。
[0080][0081]
步骤4：利用attention机制实现对不同词汇的不同关注；
[0082]
以语句向量x与各时刻隐藏层状态向量h作为attention层输入，针对目标词与其余词汇之间的相似度对每个词汇给予不同的程度的注意，最终得到融合各词汇信息和隐层状态的向量z＝(z1，...，z
t
，...，zn)；具体的，首先利用余弦距离计算目标词x
t
与文献中其余词汇xj之间的相似度分数score(x
t
,xj)，再使用softmax函数计算得出相应的归一化注意力权重参数α
t，j
，根据该权重参数与隐层向量h计算得出文献级别的全局向量g＝(g1，...，g
t
，...，gn)，其中g
t
表示目标词x
t
对应的全局向量，最后再将全局向量g与对应的隐层向量h拼接，使用tanh函数计算得出attention层的输出向量z；具体如下：
[0083]
attention层的输出作为crf层的输入。给定输入x，输出预测结果y的得分计算公式：
[0084][0085]
这里，转移矩阵元素表示标签从yi转移到y
i+1
的概率，yi为y中元素。表示第i个词语标记为yi的概率。在给定输入x情况下，输出预测结果y的概率：
[0086]
[0087]
其中，y
x
表示所有可能的标签组合，表示真实标签。模型的目标是最大化p(y|x)，这一过程通过对数似然来实现，训练过程中的似然函数：
[0088][0089]
最终预测时，输出得分最高的结果。
[0090][0091]
步骤5：使用crf层获取最优标签序列；以向量z作为tanh层的输入计算得出各词向量可能的标签对应的分数，将神经网络输出的分数记为p，p的维度为k x n，k为每个词可能的标签数量，n为词的个数；最后的crf层综合考虑本层转移矩阵t和矩阵p，计算出最优的标签序列作为标注结果。
[0092]
步骤6：根据序列标注结果，将标注为化合物、靶点及其详细信息的实体抽取保存。
[0093]
(二)构建能量参数生成模型，优化模型输出权重参数，形成多个准确的能量函数模型。
[0094]
经过命名实体识别获得活性信息后，可以获得足够数量标记活性的数据集。然后随机选择完整数据集中的部分数据作为小样本集当作机器学习模型的输入数据。通过加入了参数生成器组件的nmt进行训练，生成一组权重系数。这个权重系数描述了同一能量项在不同药-靶结合过程中能量计算的偏置情况。经过参数生成神经网络反复迭代训练此参数生成神经网络，输出最终的权重系数组合。应用此神经网络输出的权重系数组合与rosetta原有的能量项计算公式组成新的优化后的能量函数模型进行活性预测。
[0095]
能量函数近似于生物分子构象的能量，这个量称为
△etotal
，根据能量项的线性组合计算得出ei，ei为几何自由度θi、化学恒等式aai的函数，并按每个项ω的权重缩放，如以下公式所示：
[0096]
δe
total
＝∑ω
iei
(θi，aai)
[0097]
对接能量函数描述了对原子堆积、静电和溶剂化极为重要的非键合原子对之间相互作用的能量，同样具有模拟氢键和二硫键作用的潜力。能量函数同时解释了用于描述蛋白质中骨架和侧链扭转偏好的统计电位。能量函数中包含对体系结构特征概括很重要的能量术语。蛋白质的能量函数术语示例参考如表1：
[0098]
表1蛋白质的ref15能量参数术语
[0099]
[0100][0101]
本实施例针对包括但不限于以上能量项，能够对全原子进行能量参数优化。不难发现，各类生物分子能量函数所包含的原子之间相互作用的能量项可以分为范德华力、静电能、溶剂、氢键、二硫键合等多项，根据原子间的相互作用，建立明确的建模特征，训练出结合构象中主要起作用的项并优化权重，使得建立更加有针对性的准确度更高的能量函数成为可能。
[0102]
本方法采用现有神经机器翻译模型(nmt),不需要对其模型体系结构进行改变，只需要加入参数生成器组件，即可生成神经网络中的权重，采用decoupled即编码器参数根据源语言生成，解码器参数根据目标语言生成的方式，表示为：
[0103]
θ
(enc)
＝g
(enc)
(ls)，θ
(dec)
＝g
(dec)
(l
t
)
[0104]
在这种情况下编码阶段和解码阶段分离，编码器在编码时并不知道所解码的语言，编码器的中间表示是通用的，可以将其翻译为任意目标语言。
[0105]
将参数生成网络设计为简单的线性变换模式：
[0106]g(enc)
(ls)＝w
(enc)
ls[0107]g(dec)
(l
t
)＝w
(dec)
l
t
[0108]
其中，可以将其理解为对参数的低秩约束(low-rank constraint)。
[0109]
使用贝叶斯优化模型对参数生成网络进行调优，不仅能够获得较好的性能，同时比随机搜索更省时，在使用贝叶斯优化时，首先通过样本点对高斯过程进行估计与更新，然后通过选择函数来确定新的采样点。因此贝叶斯寻优的重点在于高斯过程与选择函数上。
[0110]
1.高斯过程
[0111]
一个完整的高斯过程仅由均值函数m(x)与协方差函数k(x,x
′
)确定，其中均值函数为一个向量，协方差函数为矩阵。从而高斯过程gp就可以表示为：
[0112]
f～gp(m,k)
[0113]
现假设有一组样本点d＝{(x1:t,y1:t)},其协方差矩阵为：
[0114][0115]
一个新样本xt+1加入会更新上述协方差矩阵k,假设有k＝[k(x
t+1
,x1),k(x
t+1
,x2),
…
,k(x
t+1
,x
t
)]那么更新后的协方差可表示为：
[0116][0117]
有了更新后协方差矩阵就可以通过前t个样本估计出ft+1的后验概率分布：
[0118]
p(f
t+1
|d
1:t
,x
t+1
)}～n(u,σ2)
[0119]
u＝k
t
k-1f1:t
[0120]
σ2＝k(x
t+1
,x
t+1
)-k
t
k-1k[0121]
根据新加入的样本点对先验中高斯过程进行更新，使其能更好拟合现实情况。
[0122]
2.选择函数
[0123]
确定了先验概率分布后就需要通过选择函数确定用于更新先验的采样点。这是决定贝叶斯优化能否成功的重要因素。通过该采样点获得后验分布，使该分布更贴切模拟现实情况。
[0124]
本方法使用ei准则作为选择函数，ei准则a(
·
)如下：
[0125][0126]
其中f
best
为数据集d上的最大值；e
y～f(x|d)
；ey～f(x|d)为期望函数；与φ(
·
)分别是高斯分布累计概率函数以及概率密度函数。ei准则最大优点就是能在explore与exploit两种策略间保持平衡，当explore时选择均值最大的点，exploit时选择方差大的点。
[0127]
通过对参数生成网络进行参数优化的迭代训练，对输入样本进行能量项参数预测计算，与目前多能量函数对接模块相比，这种模型具有更高的准确性和针对性。参数生成模型输入输出如下：
[0128]
输入：完整样本集p，小样本集p，迭代次数为n。
[0129]
输出：超参数空间ω。
[0130]
全原子能量函数使用的基准测试方法为对接。
[0131]
权重参数优化模型流程图如图2，具体步骤如下：
[0132]
步骤1：将小样本集使用参数生成模型进行训练，通过贝叶斯参数优化模型进行优化，输出权重系数组合，匹配生物分子构象能量计算公式，形成初步优化后的能量函数模
型。
[0133]
1)对小样本集p进行初始化操作，将输出集合置为空集。
[0134]
待优化的模型是一个加入了参数生成器的神经机器翻译模型，能够生成近似每个生物分子构象相关的能量权重系数。将小样本集p中的样本输入该模型进行训练，能够输出不同能量权重系数组合至参数空间ω。
[0135]
2)对p中每一个样本输出的能量权重系数组合通过能量函数模型r(初次迭代时该模型即为rosetta能量函数)进行评价。
[0136]
3)通过贝叶斯参数寻优模型对参数生成网络进行调优，设置迭代次数为1.形成初步优化后的能量函数模型r(1)。
[0137]
4)对小样本集p中每一个样本使用r(1)进行评价，输出预测活性值存储至输出集合a。
[0138]
步骤2：将完整数据集作为输入，重复步骤1中的工作，将集合a中的输出结果与报导的实际活性值进行比较，设置阈值，对优化后的能量权重系数生成模型是否准确给出判断。将完整样本集中验证为准确的部分作为输入，将判定为活性的样本取出，判定为活性的样本按照8：1：1的比例划分为训练集、测试集和验证集，形成新的样本集，再次对参数生成神经网络进行迭代训练。迭代训练后的能量函数模型使用验证集验证模型准确度，由于输入数据质量更高，此次训练后会输出一个准确性更高的能量函数模型。
[0139]
①
将输出集合a中通过全原子能量函数计算所得的预测活性值与实际活性值进行比较，设置阈值，超出阈值范围的样本被认为计算不准确。
[0140]
②
将完整样本集p中被判定为准确的样本放回，判定为不准确的样本从p中剔除，完成后获得新样本集p1。
[0141]
步骤3：重复步骤2，进行多次迭代，使用参数优化模型不断优化权重参数生成网络，最终形成较为准确的能量函数模型。
[0142]
1)将样本集p1输入经过多次优化后的参数生成网络，输出能量权重系数组合组成能量函数r(n)，输出活性预测值存储至集合a，参数存储至输出集合b，通过与文献报导活性值进行比较，将计算准确的样本再次输入至参数生成网络。
[0143]
2)计算准确的样本生成新一轮集合p(t)。
[0144]
3)while(未满足迭代次数为n)
[0145]
{
[0146]
进行下一次迭代，令t＝t+1，对p(t-1)重复执行步骤3中1)的操作，生成新一轮集合p(t)；
[0147]
对p(t)生成的预测活性值空间a进行评价，将参数集合b清空，存储新的参数组合；
[0148]
利用新的参数空间更新能量函数为r(t+1)。
[0149]
}
[0150]
4)输出集合b中参数的最优解。
[0151]
对于步骤1和步骤3中的原子能量函数r(t)做出解释，以范德华相互作用力为例，在函数的最小值处(d
i，j
＝σ
i，j
)分成两个可以分别加权的分量：吸引力和排斥力，通过这种方式分解函数，从而可以改变组件权重。
[0152][0153]
(三)样本聚类以及构建多分类器
[0154]
分类器组合方法能显著地降低分类器的错误率。分类的错误率的减少可以通过减少方差和偏差来实现，而分类器组合能显著的减小方差，多分类器组合技术是将多个不同的单分类器组合成一个分类器，组合的目的是利用多个分类器的差异来改善最终分类器的分类性能。
[0155]
分类模型建立的过程就是使用训练数据进行学习的过程，使用模型分类过程就是对类标签未知的数据进行分类的过程。分类系统流程图如图3所示。
[0156]
我们使用训练m个类别的数据所得到的多分类器来对训练集外数据进行判断得到结果。多分类器是组合多个分类器对实例进行分类的系统，其中每个分类器被称作基分类器。在分类阶段，每个基分类器都参与对测试实例的分类，然后运用某种组合方法，综合所有基分类器的分类结果以形成一个最终分类结果。
[0157]
本实施例样本聚类和多分类器采用的技术方案是：
[0158]
经过步骤(二)之后，可以构建一个基于能量函数的活性预测模型，此模型对属于某分子族或具有某种特性或者结构的化合物具有良好的预测准确度。相对来说，这类分子在活性表现上最接近于最后一次用于迭代训练能量函数模型的样本集。
[0159]
步骤1：在本实施例中，从完整样本集中随机选取适量小样本集时随机抽取多次，使用(二)中的方法训练出对应次数数量的能量函数模型。同理，这些模型对其最后一次迭代训练所输入的数据集具有相对较好的活性预测准确度。
[0160]
步骤2：假设迭代训练了m个活性预测模型，则对应m个子样本集，每个样本集在特征上具有一定的相似性，此过程将完整的样本集简单聚类成m个分子族此分类器的作用是将新的化合物分类成m个分子族之一。
[0161]
1)分子、靶点编码
[0162]
在步骤(一)获取文献中的活性信息时，可同时获得分子的标准格式文件。分子文件可由rdkit等开源工具转换成smiles文件，其中包含了结构信息。在得到分子的唯一序列表示后，我们利用自编码器(autoencoder)对分子信息进行编码，通过压缩数据提高数据分类的效率。autoencoder是一种能够通过无监督学习，学到输入数据高效表示的人工神经网络。
[0163]
autoencoder由输入层，隐藏层和输出层构成，如图4所示：
[0164]
它利用反向传播算法尝试学习一个h
w,b
(x)≈x的函数让目标值等于输入值。autoencoder尝试逼近一个恒等函数，使得输出接近于输入x。为了使这个函数有意义，需要加入一些限制条件(比如说限制隐藏神经元的数目)，这里我们引入稀疏性限制来限制隐藏层。对于使用sigmoid作为神经元的激活函数的情况下，若神经元的输出为11表示该神经元被激活，否则称为未被激活，则稀疏性的含义是指在大多数情况下神经元都是未被激活的。可以使用神经元的输出作为神经元的激活度，即对于隐藏层的第j个神经元，其激活度为：
[0165]
[0166]
则对于m个样本，其平均激活度为：
[0167][0168]
假设令其中ρ是一个常数，表示的是稀疏性参数，通常可以取一个接近于0的常数，如取ρ＝0.05。为了使得模型比较稀疏，我们希望平均激活度能够尽可能接近稀疏性常数，通常可以取相对熵来度量平均激活度与稀疏性参数之间的差异程度。熵的公式如下：
[0169][0170]
对于上述的自编码器模型，其隐含层的第j个神经元的平均激活度为：
[0171][0172]
稀疏性常数为：ρ，则对于隐含层的相对熵为：
[0173][0174]
其中，s2表示的是隐含层节点的个数，相对熵又称为kl散度，即：
[0175][0176]
相对熵是一种用来度量两个分布之间的差异的方法。对于上述的相对熵，若时达到最小值，最小值为0，否则差距越大，相对熵的值越大。对于稀疏自编码器的损失函数，其与神经网络的损失函数一致，可以表示为：j(w,b)
[0177]
则对于稀疏自编码器，其损失函数即在神经网络的损失函数的基础上加上稀疏性的约束即可，即为：
[0178][0179]
其中，β控制稀疏性的权重。在更新的过程中，原本在神经网络中，更新公式为：
[0180][0181]
对于稀疏自编码器，其更新公式为：
[0182][0183]
2)生成基分类器
[0184]
由于分类器是由分类模型在训练集上训练得到的。通过使用不同的分类模型，或者不同的训练集都可以产生不同的分类器。相应地，当前存在的基分类器的生成的方法也基本都属于这两类，我们采取的是第二类方法，使用相同的分类模型，而通过对训练集中的实例，属性，类标等的变化来形成不同的训练输入，以最终形成不同的分类器。当给出原始的训练集合后，可以通过删除，添加，抽取等方法来操作训练集合中的部分实例来形成新的训练集合。
[0185]
我们使用boosting方法来实现产生基分类器。具体的，我们赋给每一个类别一个权重，根据权重的大小决定其被抽取的概率，且初始的样本权重是相等的，在每次迭代中权重都会发生变化。运行过程是：我们首先赋予我们得到的结构相似的m个类别组成的数据集中的每一个类别相同的权值α。
[0186][0187]
然后应用弱学习算法，在这些加权的训练集合上学习出一个分类器，再根据这个分类器的错误率ε对每一个类别的数据重新赋权。赋权的原则是增加上次被分类器错误分类的数据的权值，相应的减小被分类器正确分类的实例的权值。调整后分类正确的样本权重会降低，分类错误的样本权重会上升，更新完权重再进行迭代训练。
[0188]
若令d表示分布，某个样本被正确分类，则权重更新为：
[0189]d(m+1，i)
＝d
(m，i)
*exp((-α)/sum(d))
[0190]
某个样本被错误分类，则权重更新为：
[0191]d(m+1，i)
＝d
(m，i)
*exp((α)/sum(d))
[0192]
boosting生成分类器方法对应图5如示
[0193]
3)选取基分类器
[0194]
当生成多个基分类器后，需要对分类器进行筛选。使用经过某种方法挑选出来的基分类器子集形成的组合，跟使用全部的基分类器相比，可能具有更好的分类性能。我们使用的是分类器的动态选取技术。
[0195]
具体地，第一步，首先解决的是如何决定测试数据所属的特征空间。我们已经得到了原子类型、分子的结构、化学键等属性，那么就可以利用属性值的统计信息，用属性对评价数据集进行划分。
[0196]
具体的，第二步，我们在这些划分上测定每个基分类器的分类正确率。在分类阶段，根据测试数据的属性值，确定对应的划分，并计算每个基分类器在这些划分上的平均分类正确率，正确率最高的被选出来对测试数据进行分类。
[0197]
4)组合基分类器
[0198]
在确定了多分类器组合的基分类器后，如何组合这些基分类器对测试数据进行分类是多分类器组合中的一个基本问题。在确定采用哪种组合方法之前，首先要确定基分类器的分类结果的表示形式。不同的组合方法适用于不同的表示形式。当前，分类结果的表示形式主要有：
[0199]
(a)基分类器的分类结果仅仅给出一个类标。
[0200]
(b)基分类器的分类结果是一个所有类标按可能性大小的一个排序。
[0201]
(c)基分类器的分类结果是一个向量，向量的每一个分量都给出了可能是每一个类标的程度大小，一般用概率来表示。
[0202]
针对这三种输出结果形成，我们选择采用均值法。具体地，对每一个类标y
x
，首先将所有基分类器给出的预测类标可能是它的概率相加并求平均值，然后这些值最大的类标作为最终的分类结果。分类公式如下：
[0203]
[0204]
(四)整体系统框架图如图6所示。
[0205]
应用实例
[0206]
1、从文献数据库中获得相关文献，然后输入ner模块，获取输出的样本活性信息，可以构建完整样本集。
[0207]
2、随机选取小样本集迭代训练出多个能量函数模型，具体过程如下：随机选取样本组成样本集用于训练；选取能量计算项，使用参数生成模型为能量项生成权重参数，使用贝叶斯优化方法对权重参数进行调整，输出一组较优权重系数组合；权重系数组合分别乘以能量计算项形成新的能量函数活性预测模型；每次把完整样本集输入预测模型，将高于活性阈值的样本再次输出入到参数生成模型中，用以优化权重系数组合，从而优化能量函数模型，此过程迭代进行，直至本轮优化模型后判断活性样本数量与输入的上轮判断活性样本数量达到95％以上为止；上述所有过程执行多次，形成对应数量的能量函数预测模型。
[0208]
3、把2中每个能量函数模型最后一次迭代用以训练的样本当成多个特征的样本集，编码并训练多分类器。新化合物应用本方法预测活性之前先分类成多个分子族之一，再选用对应的能量函数模型进行预测，最终输出结果。
[0209]
4、选择已验证药物测试本虚拟筛选方法的准确性。
[0210]
(1)在研发针对hiv/aids治疗药物的过程中，许多方案都是对反转录酶进行抑制，分为核苷类和非核苷类。核苷类反转录酶抑制剂在dna组装过程中通过转录酶的作用被组装到不断增长的dna链中，由于在核糖部分c-3’位置缺乏oh基团，使其不能与磷酸二酯键结合，此类药物的组装干扰了病毒链的增长。如图7所示，应用中选用齐多夫定、拉米夫定两种药物两种核苷类抑制剂。非核苷类反转酶抑制剂不是结合在酶的活性位点，而是在酶活性位点附近的疏水口袋，应用中选用利匹韦林和依曲韦林两种非核苷类化合物进行测试。
[0211]
(2)如图8所示，用于治疗阿尔兹海默症的抑制剂(a)具有一定的酶和细胞水平活性，在p’侧链引人一个羟基有助于形成氢键作用，由此得到的抑制剂(b)、(c)具有更高的bace1抑制活性。抑制剂(a)相应的甲醚衍生物(d)的抑制活性同样下降。此外，去掉一级侧链的羰基后，化合物的bace1抑制活性丧失。
[0212]
(3)如图9所示，用于直接凝血酶抑制剂治疗血栓的化合物(a)的活性相对偏低，需进一步的提高。经研究发现，通过脂肪链、烷氧基链和烷基胺链桥连得到的化合物(b-d)，其活性分别都提高了将近一个数量级。其中，胺基衍生物(d)在本系列化合物中活性最高。
[0213]
(4)大分子药物以糖类为例，其分子式均可以写成c
x
(h2o)
x
。单糖可以通过各种方式互相连接在一起形成多糖(或寡糖，又称低聚糖)，是一种长链分子。许多糖类药物可以含有一个或多个基团被其他基团取代(可以是蛋白、脂类)或移除。例如图10所示，几丁质是一种被重复的n-乙酰氨基葡萄糖(一种含氮原子的葡萄糖)片段所组成的糖类大分子。
[0214]
与其他分子相比，糖类分子特有的糖苷键分子内能量使得糖类在进行活性计算时具有特定的灵活性。如图11所示，与其他小分子相比，多糖分子由于是由多个单糖相互连接在一起，分子之间具有较高的结构相似性，单糖之间存在结构模仿，其结构的特殊性很容易被本发明多分类器识别出与其他多肽、脂类以及小分子药物的区别，从而在迭代训练能量函数模型时能够有效构建针对此类分子的活性预测模型，在多分类器实现新化合物匹配对应能量函数模型时能够被正确分类。
[0215]
另选一糖类大分子——当归多糖铁复合物作为测试药物之一。当归多糖铁复合物
由当归多糖和三价铁组成的大分子复合物，其结构为以三价铁通过氧桥和羟基桥聚合而成ferrihydrite聚合铁核为结构中心，聚合铁核上稳定地螯合一层当归多糖链(asp)32,形成铁核分子，铁核分子外部再包绕一层亲水性的鞘状当归多糖链(asp)12,分子式为{[(fe2o3
·
2.2h2o)1043(asp)32](asp)12},分子量为270 000da。
[0216]
5、选择以上已验证几种药物测试本虚拟筛选方法的准确性，将以上分子作为待预测的新化合物输入模型进行活性预测。首先将其分子式转化成smiles格式，将其输入多分类器。分类结果如下：用于治疗阿兹海默症的bace1和用于治疗hiv的核苷类药物齐多夫定、拉米夫定被分成一类，核苷类药物利匹韦林和依曲韦林被分成一类，凝血酶抑制剂被分成一类，两种糖类大分子被分成一类。
[0217]
6、选用对应的能量函数预测模型进行活性预测，输出活性预测结果，与实际活性值进行比对，验证那本虚拟筛选方法的有效性。将上述化合物使用对应的能量函数模型进行能量计算后，输出结果为数值e，单位：(kcal/mol)，此数值的绝对值越大，活性越高。将此结果与晶体数据集中真实活性值ea的绝对值比较，误差比越小，代表计算准确度越高。设置误差比10％为活性阈值，小于10％则认为活性预测准确。其误差比f计算如下：
[0218][0219]
例如聚糖类分子对接分数为669.728kcal/mol，真实晶体实测值为733.140kcal/mol，则误差比为8.65％。经检验，上述所有分子预测误差比均不超过10％，该能量函数模型可以有效实现活性预测。
[0220]
结果分析：在迭代训练模型的过程中，不断输入判断为化合物，有指向性的使模型向具备某些特征、属性的分子族进行活性预测的拟合，是模型针对特定种类分子进行活性预测时具有较高的精度。这些化合物可能具备相似的骨架(如品字形环状化合物)或者相同分子基团、片段，其相似结构的表现出的某一或某些化学活性特性被能量函数模型捕捉。这时多分类器会将具备某些相同特点的分子划分为同一个已有的类别之一，再用该类别分子对应的能量函数进行预测，这类能量函数模型可能对带有特殊基团或者相似骨架的分子较为敏感，最后输出较为理想的活性结果。
技术特征：1.一种结合机器学习和构象计算提高药-靶活性预测精度的方法，其特征在于，所述方法包括命名实体识别获取文献中活性数据集，训练神经网络以输出权重参数，样本聚类以及构建多分类器；具体步骤如下：(一)命名实体识别获取文献中活性数据以获得活性数据集，包括以下步骤：1)从现有论文数据库中获取活性化合物和靶点的文献信息；2)使用pubmedbert预训练模型将输入文献信息向量化；3)使用bilstm获取包含上下文信息的状态向量；4)利用attention机制实现对不同词汇的不同关注；5)使用crf层获取最优标签序列；6)根据序列标注结果，将标注为化合物、靶点及其详细信息的实体抽取保存；(二)构建能量参数生成模型，优化模型输出权重参数，形成多个准确的生物分子能量函数模型包括以下步骤：1)从(一)中获得的活性数据集抽取一部分作为小样本集，将小样本集使用加入参数生成器的神经机器翻译模型进行训练，生成一组权重系数，通过反复迭代训练此参数生成神经网络，输出最终的权重系数组合，所述的参数生成神经网络是指所述加入参数生成器的神经机器翻译模型，应用此参数生成神经网络输出的权重系数组合与rosetta原有的能量项计算公式组成新的优化后的能量函数模型进行活性预测，公式如下：δe
total
＝∑ω
i
e
i
(θ
i
，aa
i
)能量函数近似于生物分子和靶蛋白耦合构象的能量，这个能量称为δe
total
，根据能量项的线性组合计算得出ei，ei为几何自由度θ
i
、化学恒等式aa
i
的函数，并按每个项ω的权重缩放；2)将(一)中获得的完整的活性数据集输入步骤1)所得的能量函数模型，将输出结果与文献报导的实际活性值进行比较，设置活性阈值，在能量函数模型中验证，将判定为活性的样本取出，把高于活性阈值而被判定为活性的样本按照8∶1∶1的比例划分为训练集、测试集和验证集，形成新的样本集，使用训练集再次对参数生成神经网络进行迭代训练，输出新的权重系数组合，从而形成新的能量函数模型；经过迭代训练后的新能量函数模型使用测试集用以测试能量函数模型的泛化误差，验证集用以评估模型准确度以及设定相关超参数；3)重复步骤2)，进行多次迭代，不断优化参数生成神经网络，验证准确度时设置阈值，达到阈值时终止步骤2)，最终形成一个相对较为准确的能量函数模型；重复步骤1)、2)、3)，从完整样本集中随机选取小样本集随时抽取多次，获得多个对具有某些特征有偏好识别能力、更敏感的能量函数模型，所获得的能量函数模型数量与随机抽取样本集的次数一致；这些模型对迭代训练过程中其最后一次训练所输入的数据集具有相对较好的活性预测准确度，在此依据活性表现将最后一次训练所输入的数据完成样本聚类；(三)构建多分类器构建多分类器以区分新化合物分子，以匹配(二)中获得的相对应的能量函数模型预测活性，具体构建方法如下：1)使用自动编码器进行分子、靶点编码
2)使用boosting方法和softmax方法生成基分类器3)选择、组合分类器。2.根据权利要求1所述的一种结合机器学习和构象计算提高药-靶活性预测精度的方法，其特征在于，所述方法所述(一)中步骤2)的具体操作是：采用基于pubmedbert的预训练模型将包含m个语句的输入文献d＝(x1，...，x
t
，...，x
m
)向量化，得到m个包含上下文语义信息的向量x＝(x1，...，x
t
，...，x
n
)，其中n为语句中单词的数量，x
i
为语句中包含第i个词上下文信息的词向量。3.根据权利要求1所述的一种结合机器学习和构象计算提高药-靶活性预测精度的方法，其特征在于，所述(一)中步骤3)的具体操作是：以语句向量x作为bilstm模型的输入，提取每个词向量x
i
的上下文特征得到输入的语句向量对应的隐藏层状态向量h＝(h1，...，h
t
，...，h
n
)，其中h
t
为t时刻隐藏层状态向量。4.根据权利要求1所述的一种结合机器学习和构象计算提高药-靶活性预测精度的方法，其特征在于，所述(一)中步骤4)的具体操作是：以语句向量x与各时刻隐藏层状态向量h作为attention层输入，针对目标词与其余词汇之间的相似度对每个词汇给予不同的程度的注意，最终得到融合各词汇信息和隐层状态的向量z＝(z1，...，z
t
，...，z
n
)；具体的，首先利用余弦距离计算目标词x
t
与文献中其余词汇x
j
之间的相似度分数score(x
t
，x
j
)，再使用softmax函数计算得出相应的归一化注意力权重参数α
t，j
根据该权重参数与隐层向量h计算得出文献级别的全局向量g＝(g1，...，g
t
，...，g
n
)，其中g
t
表示目标词x
t
对应的全局向量，最后再将全局向量g与对应的隐层向量h拼接，使用tanh函数计算得出attention层的输出向量z。5.根据权利要求1所述的一种结合机器学习和构象计算提高药-靶活性预测精度的方法，其特征在于，所述(一)中步骤5)的具体操作是：以向量z作为tanh层的输入计算得出各词向量可能的标签对应的分数，将神经网络输出的分数记为p，p的维度为k x n，k为每个词可能的标签数量，n为词的个数；最后的crf层综合考虑本层转移矩阵t和矩阵p，计算出最优的标签序列作为标注结果。6.根据权利要求1所述的一种结合机器学习和构象计算提高药-靶活性预测精度的方法，其特征在于，所述(三)中步骤1)的分子、靶点编码：在(一)获取文献中的化合物的活性信息时，可同时获得分子的标准格式文件；分子文件由开源工具转换成smiles文件，其中包含了结构信息，靶点文件可开源获得pdb或者其他格式的构象文件，使用autoencoder的前馈神经网络将二者编码。7.根据权利要求1所述的一种结合机器学习和构象计算提高药-靶活性预测精度的方法，其特征在于，所述(三)中步骤2)的具体操作是：使用相同的分类模型，而通过对训练集中的包括实例、属性和类标的变化来形成不同的训练输入，以最终形成不同的分类器；当给出原始的训练集合后，通过包括删除、添加、抽取在内的方法来操作训练集合中的部分实例来形成新的训练集合；使用boosting方法来实现产生基分类器。8.根据权利要求1所述的一种结合机器学习和构象计算提高药-靶活性预测精度的方法，其特征在于，所述(三)中步骤3)的具体操作是：首先解决如何决定测试数据所属的特征空间；得到包括原子类型、分子的结构、化学键
在内的数据属性，利用属性值的统计信息，用数据属性对评价数据集进行划分；然后，在这些划分上测定每个基分类器的分类正确率，正确率最高的被选出来对测试数据进行分类。9.根据权利要求1所述的一种结合机器学习和构象计算提高药-靶活性预测精度的方法，其特征在于，所述(三)中步骤4)的具体操作是：在确定采用哪种组合方法之前，首先要确定基分类器的分类结果的表示形式，对每一个类标y
x
，首先将所有基分类器给出的预测类标可能是它的概率相加并求平均值，然后这些值最大的类标作为最终的分类结果；分类公式如下：
技术总结本发明涉及一种结合机器学习和构象计算提高药-靶活性预测精度的方法，属于药物筛选技术领域，所述方法包括命名实体识别获取文献中分子数据集，训练神经网络以输出权重参数，样本聚类以及构建多分类器；本发明方法结合了机器学习算法与晶体构象能量计算的方法提高虚拟筛选的精度，使得药物筛选更高效、成本更低，且结果的准确度和可靠性提高。且结果的准确度和可靠性提高。且结果的准确度和可靠性提高。
技术研发人员：刘昊 周源东 陈淼 王晓薇 夏祎敏 刘其琛
受保护的技术使用者：中国海洋大学
技术研发日：2022.04.22
技术公布日：2022/7/5
