技术领域
[0001] 本发明涉及一种抗乳腺癌候选药物的效果预测模型构造方法，属于预测模型技术领域。
相关背景技术
[0002] 乳腺癌是全球范围内女性最常见的恶性病症，严重威胁女性的生命与健康，是多种致癌因子作用在乳腺上皮细胞上发生增殖失控的现象。近年来乳腺癌治疗手段日趋成熟，靶向治疗、手术治疗、化学治疗等治疗方法在治疗乳腺癌方面非常适用，但因乳腺癌易变异易转移且治疗方法副作用较大，所以因乳腺癌死亡的人数仍然得不到有效控制。因此，挑选有效的治疗药物以及候选药物对降低死亡率具有重要意义。
[0003] 大约75％的人受雌激素的影响可以接受靶向雌激素(此处为抑制ERα)的治疗，ERα生物活性值越低代表拮抗ERα生物活性越大，抑制ERα活性作用越大；但由于乳腺癌细胞的一系列耐药性，只有分子靶点阳性的患者才可以使用靶向药物治疗，而且由于细胞的耐药性，在一段时间后可能会有至少一半的患者出现耐药性，并非对所有患者有效。所以一个化合物想成为乳腺癌候选药物除了需要考虑拮抗乳腺癌雌激素活性，还要考虑药物在人体内的药代动力学性质和安全性，即ADMRT性质，说明该化合物要在小肠上皮细胞渗透性好、在人体中较易代谢、更适合口服、不具有毒性以及遗传毒性。
[0004] 乳腺癌疾病的种类繁多,病因复杂,在医学研究和开发对抗乳腺癌毒品问题已经不容忽视。药物研究与发展需要选择潜在活性化合物作为候选药物。复合成为候选药物的标准不仅需要良好的生物活性并且需要在人体内具有良好的药代动力学性质和安全性。本文根据提供的ERα受体拮抗剂信息来构建生物活性的定量分类预测模型和ADMET性质预测模型,筛选出具有抑制ERα的尽可能好的生物活性分子描述符和具有良好的ADMET性质和范围的化合物。最终不断优化模型预测能力以至实现乳腺癌药物治疗的效果达到最好。
具体实施方式
[0052] 为使本发明的目的、技术方案和优点更加清楚明了，下面通过附图及实施例，对本发明进行进一步详细说明。但是应该理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限制本发明的范围。
[0053] 除非另有定义，本文所使用的所有的技术术语和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同，本文中在本发明的说明书中所使用的术语只是为了描述具体的实施例的目的，不是旨在于限制本发明。
[0054] 根据研究目标，挑选出能够使化合物对抑制ERα更好的生物活性且具有更好的A、D、M性质的分子描述符。为了结果的精准性，变量需要进行排序并筛选出对生物活性影响最具显著影响的20个分子描述符，进而通过筛选后的分子描述符建立预测模型。一方面，根据筛选后的分子描述符构建化合物对ERα生物活性的定量预测模型，由于变量之间存在高度非线性关系，故传统的预测方法不再适用。BP神经网络具有较强的非线性映射能力和较高的计算精度。因此，我们对数据进行标准化处理之后从1974个化合物中(附件数据集)抽取90％的数据用于神经网络的训练，同时采用剩下的10％的数据验证模型的精度评价模型的合理性并通过均方误差说明模型的拟合效果。通过构建的模型根据10个化合物筛选出的20个分子描述符对其IC50和PIC50进行预测；另一方面，根据原始的分子描述符和化合物的ADMET数据构建化合物Caco‑2、CYP3A4、hERG、HOB、MN的分类预测模型，对其中10个化合物进行相应的预测。最后要得到分子描述符的取值以同时满足最佳的抑制ERα的生物活性和尽可能好的ADMET性质。
[0055] 根据ERα拮抗性信息，通过数据挖掘、数据分析方法，对化合物的生物活性以及ADMET性质分别进行定量预测和分类预测，同时能够对优化过程进行积极响应，提供相应的预测技术。本文技术路线图如图1所示：
[0056] 本申请模型实现依据随机森林集成模型，随机森林是由统计学家Leo Breiman在20世纪90年代最早提出的一种基于Bagging思想的集成模型。它是一个包含多个决策树的分类器，其输出的类别由单棵决策树输出类别的众数而定。
[0057] 随机森林在使用时一般通过基尼系数作为切分点的标准，通过类权重来决定类标签。两者主要由随机森林加权权重决定，加权权重相当于给予特征贡献小的类较大的惩罚系数，同时给予特征贡献大的类较小的惩罚系数。
[0058] 用于切分点的不纯度值计算如下：
[0059]
[0060] ΔG＝G(NL)‑G(NR)(2)
[0061] 决定终节点的类标签计算表达式如下：
[0062]
[0063] 其中N为未分离的节点，NL和NR分别为分离后的左节点和右节点，Wi为类样本权重，ni表示各节点各类样本的数量，ΔG衡量不纯度的减少量，其值越大表明分离效果越好。
[0064] 对于变量的贡献度评价用V表示，GINI值用GI表示。对于m个特征值X1,X2,…Xm，计算每个特征值的基尼系数：
[0065]
[0066] 其中，K表示有K个类别(本文设定K＝100),ρmk表示节点m中第k个类别所占的比例。
[0067] 特征xJ在节点m分枝前后基尼系数变化量为：
[0068]
[0069] 其中GIl和GIr表示分枝后两个新节点的基尼系数。
[0070] 如果特征Xj在决策树i中出现的节点在集合M中，那么Xj在第i棵树的重要性为：
[0071]
[0072] 从而随机森林中所有树的基尼系数总和为：
[0073]
[0074] 为了使所有特征贡献度总和为1，故需要对每一个特征重要性进行归一化处理，表达式如下：
[0075]
[0076] Pr{Pw＝0}＜α＜1 (9)
[0077] 通过随机森林筛选后，考虑到自变量之间具有高度的非线性耦合性，耦合性越高代表独立性越差，降低变量间的耦合程度能大大减少模块间的影响，因此有必要对数据进行再次过滤。从相互高度相关的变量中保留贡献值较大的变量，以此代表高耦合变量组(即高度相关的变量组成的变量组)的全部信息。
[0078] 传统的Pearson相关是英国统计学家皮尔逊提出，只适用于两个变量是线性关系的情况。但本题中各变量之间使复杂的非线性关系，故Pearson相关分析不适用。因此通过Spearman秩相关分析寻找变量间的相关性指标。
[0079] Spearman秩相关系数用来估计两个变量Y1,Y2的相关性，是衡量两个变量依赖性的非参数指标，变量之间的Spearman相关系数满足以下公式：
[0080]
[0081] 相关系数越大，两者相关性越高。在统计学中，对相关性强弱有如下(见表1)说明；
[0082]
[0083] 本文采用一种迭代更新的算法，首先找寻自变量和因变量之间高度相关的分子描述符以保证预测模型的合理性。然后找寻自变量和因变量之间相关程度尽可能低的分子描述符以减少变量之间的非线性耦合性影响。
[0084] 本发明还采用BP神经网络，BP神经网络是一种多层神经网络的反馈，具有任意复杂的模式分类和优异的非线性函数映射能力。在结构上，BP神经网络具有输入层、隐藏层和输出层。在不同情况下，层数及各层神经元个数可任意设定，随着结构的变化其性能也会有所不同。在转移的过程中,输入信号从输入层进入，逐层通过隐藏层进行处理,直至输出层。每一层的神经元状态将直接影响下一层神经元的状态。如果输出层不能满足实际的预期数据则转换为反向传播,根据预测误差改变权重和阈值,以便接近预期的输出结果。
[0085] BP神经网络中的每一个输入量经过隐含层后，相应的权值和阈值都有一定的变化，在隐藏层加权后到达输出层进行输出。根据BP神经网络原理对于各层有：
[0086] 输出层：
[0087] hk＝f(netk),k＝1,2,…,m (11)
[0088]
[0089] 隐含层：
[0090] yj＝f(netj),j＝1,2,…,m (13)
[0091]
[0092] 其中:
[0093]
[0094] 设计神经网络结构示意图如图2所示。
[0095] 本发明还采用遗传算法进行结合，遗传算法是一种基于自然选择和遗传学原理进化出的计算模型，利用优胜劣汰的自然选择和生物界繁衍的基因重组、突变的遗传机制的全局自适应概率搜索算法，通过数学的方式进行计算机仿真运算。遗传算法问题中潜在解集的一个种群开始，按照适者生存的原理，逐步演化产生出几乎最优的解。在进化过程中，由于个体适应度大小选择个体，借助遗传算子进行选择、交叉、变异，从而产生出新的种群，因此进化出的后代种群会比前代更加能适应环境，最后的结果也是问题中的最优解。
[0096] 在种群中选择优质个体、淘汰劣质个体的操作叫选择，目的是为了把优化个体遗传给下一代。同样生物遗传基因的重组也起了核心作用，即把两个父代个体的部分结构加以替换重组产生新的个体。变异算子作为交叉算子的辅助算子，结合交叉算子，使得交叉和变异相互配合和竞争，使其具备兼顾全局和局部均衡搜索的能力。
[0097] 计算时将实际问题中的变量进行编码形成染色体，随机产生种群，计算出每个个体的适应度，然后通过终止条件判断是否是最优解，若是则结束计算，若不是则重新进行遗传操作产生新的种群，此过程循环进行，直至选出满足优化准则的最优解。
[0098] 本实施例中，建立了随机森林筛选模型，对给定的1974个化合物的729个分子描述符进行初步变量选择，并弱化了自变量间的强耦合性，筛选出了对生物活性具有显著影响的分子描述符。具体来说，为了提高模型筛选效率并保证筛选的合理性，先对原始数据进行预处理，剔除了其中存在的大量零值，进一步对保留的有效数据进行归一化处理，而后通过随机森林对归一化处理后的数据进行初步筛选，针对1974个化合物的729个变量分别计算自变量贡献度和累计贡献度，对排名前70的自变量进行可视化分析，如图3所示。
[0099] 上图中变量贡献度在第70个变量后均小于0.0023，所以选取前70个自变量做下一步分析。根据自变量和因变量的相关系数进行第一步筛选，最终选出与因变量之间高度相关的38个自变量如表2所示：
[0100]
[0101] 初次筛选后，保留下来的38个自变量能够充分解读。但是由于自变量之间的耦合性导致自变量之间表达的信息重合度高，对资源造成浪费。故需要通过Spearman等级相关分析对变量进行二次筛选以实现信息高度表达的前提下使得变量维数尽可能小。初次筛选后自变量之间的Spearman相关系数如图4所示，方形框中为相关性较高的变量组。
[0102] 经过不断更新迭代筛选，最终得到20个对生物活性具有显著性影响的分子描述符，其相关性热力图如图5所示。
[0103] 根据以上步骤，对1974个化合物的729个分子描述符进行筛选后，根据变量对生物活性影响的重要性(贡献度)进行排序，得到对生物活性最具显著影响的20个分子描述符以及贡献度如表3所示：
[0104]
[0105] 考虑到BP神经网络具有较高的计算精度，本文通过筛选出的20个分子描述符，建立了基于BP神经网络的IC50值和对应的PIC50值的预测模型。从1974个原始样本中随机挑取80％用于模型训练，得到了训练后生物活性定量模型的重要参数。再利用剩余20％的样本进行预测，借助均方误差分析了模型的拟合效果并验证了模型的合理性。
[0106] 本文采用含有一个隐藏层的三层多输入结构构建用于化合物对ERa生物活性的定量预测模型。
[0107] (1)数据整理
[0108] 根据第一题中的筛选方法对本题中50个化合物对应的729个分子描述符进行筛选，对筛选出的20个分子描述符进行归一化处理。将每一行作为一组输入训练集，将数据中1974组数据按80％和20％的比例分为训练数据(80％)和测试集数据(20％)。
[0109] (2)设计网络层结构
[0110] 输入输出层：通过筛选得到20个分子描述符作为输入，PIC50值作为输出。故输入层神经元个数n＝20，输出层神经元个数m＝1。
[0111] 隐藏层：隐藏神经元个数的确定十分重要。隐藏神经元个数的设定过多，会加大网络计算量并且容易导致系统过拟合。隐藏层神经元个数过少，则会影响网络的性能，达不到预期的效果。网络中隐藏层神经元的数目与实际问题的复杂程度、输入和输出层的神经元数以及对期望误差的设定有着直接的联系。目前，没有明确的公式可以确定隐藏层中神经元的数目。本文出于对精度的考虑，通过多次调试，最终确定隐藏层神经元为15个。
[0112] (3)激活函数的选取
[0113] 在多层神经网络中，上层节点的输出和下层节点的输入之间具有一个函数关系，这个函数称为激活函数。激活函数为BP神经网络提供了非线性化泛化能力，使得神经网络可以将每一层的输入和输出建立非线性联系。如果没有激活函数，则每一层节点的输入都是上层输出的线性函数。本文中采用一种双曲正切Sigmoid函数(tansing)，在tansing函数中，输出将被限制在(‑1，1)区间内。使用计算公式为：
[0114]
[0115] (4)均方误差
[0116] 模型的预测问题本质上属于一种回归问题，与分类问题不同，解决回归问题的神经网络一般只有一个输出，故模型结果的评价标准通常用均方误差(MSE)来衡量。其定义为：
[0117]
[0118] 其中yi为数据的真实值，yi'为神经网络给出的预测值。
[0119] 首先采用绍的随机森林的方法进行预测输出，得到预测结果的均方误差为3.0759，其预测输出和实际输出的对比如图6所示。
[0120] 从图中可以看出随机森林的预测效果并不理想，所以通过BP神经网络进行预测。采用三层结构的BP神经网络进行训练，将筛选出的20个分子描述符作为输入变量，设置最大迭代次数为100次，期望误差为0.00004，学习速率设置为0.1，经过1974个化合物的不断训练(其中80％)和测试(其中20％)，得到预测结果的均方误差为0.514，该值较量纲选取而言较小。其预测输出和实际输出的对比如图7所示。
[0121] 从图中可以看出预测输出的值和实际输出的值重合部分较多，说明预测结果较为理想，建立的预测模型较为合理。为了更清晰的描述误差的取值情况，做出每个样本点的误差分布及其变化情况，从而得到神经网络的预测误差图如图8所示。
[0122] 最终构建得到输出为PIC50值的预测模型，根据模型对50个化合物进行IC50值和PIC50值的预测。
[0123] 构建分类预测模型：
[0124] 本文基于BP神经网络模型，利用筛选方法得到5种ADMET性质对生物活性最具显著影响的20个分子描述符，根据三种模型发给发分别建立有关化合物Caco‑2、CYP3A4、hERG、HOB、MN的分类预测模型。
[0125] 为了构建化合物5种ADMET性质的分类预测模型，首先根据上文的筛选方法分别找到5种ADMET性质对生物活性最具显著影响的20个分子描述符，其贡献度如表5所示：
[0126]
[0127] 对50个化合物进行相应的预测。本文采取BP神经网络模型构建化合物Caco‑2、CYP3A4、hERG、HOB、MN的分类预测模型。根据附件所提供的1974个化合物对应729个分子描述符，随机抽取80％用于模型的训练，20％用于模型的测试，得到BP神经网络模型预测的准确率如表6所示。
[0128]分子描述符 Caco‑2 CYP3A4 hERG HOB MN
准确率 0.867816 0.913793 0.896552 0.804598 0.913793
[0129] 根据预测效果，最终选取BP神经网络针对5个化合物分别构建5个分类预测模型，对50个化合物进行相应的预测，得到部分预测结果如表7所示。
[0130]
[0131] 为了寻找同时满足抑制ERa具有更好的活性以及具有更好的ADMET性质的分子描述符并确定其取值，本文首先采用问题1中的化合物筛选方法初次筛选出20个分子描述符，然后借助Logistic检验的方法得到7个显著变量，按显著能力排序分别为MDEC‑23、LipoaffinityIndex、C1SP2、TopoPSA、VC‑5、ETA_Shape_Y、SHBint10并对其进行归一化。根据归一化后的7个分子描述符对1974个化合物样本中的80％进行神经网络训练，输出PIC50和性质A、D、E的取值，通过剩余20％的样本进行网络测试，经过隐藏层节点调试，最后构建出输入层为7个自变量，隐藏层为12个节点，输出层为4个因变量的神经网络模型。如图9所示。
[0132] 最后根据遗传算法随机生成7个自变量的取值进行归一化作为所构建神经网络模型的输入，将模型输出的PIC50、A、D、E的值求和作为目标函数，通过最大化目标函数优化从而得到满足题目要求的自变量取值并将其反归一化，具体结果如表8所示：
[0133]
[0134] 由表中数据，寻找的分子描述符分别为MDEC‑23、LipoaffinityIndex、C1SP2、TopoPSA、VC‑5、ETA_Shape_Y、SHBint10，其对应取值为38.4、23、20、1059.1、1.5、0.6、0时，能够使化合物对抑制ERα具有更好的生物活性，同时具有更好的A、D、E性质。
[0135] 以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换或改进等，均应包含在本发明的保护范围之内。
