1. 引言
乳腺癌目前发病率在2~3/万之间,高居全球第一,而且还在呈上升趋势,年龄也越来越年轻化。在研究ERα基因缺失小鼠的实验结果中,发现ERα是治疗乳腺癌的重要靶标,因此能够拮抗ERα活性的化合物可能是治疗乳腺癌的候选药物 [1] 。但是想要成为候选药物,除了需要具备良好的生物活性外,还需要在人体内具备良好的药代动力学性质和安全性,合称为ADMET (Absorption吸收、Distribution分布、Metabolism代谢、Excretion排泄、Toxicity毒性)性质 [2] 。但一个化合物的活性再好,如果其ADMET性质不佳,比如很难被人体吸收,或者体内代谢速度太快,或者具有某种毒性,那么其仍然难以成为药物,因而还需要进行ADMET性质优化。
传统药物研发渠道的平均成本为26亿美元,大概耗时12年,因此如何在降低成本和时间的同时确保药物的有效性成为药物公司的重大难题,基于机器学习、深度学习辅助药物各个阶段的研发越来越成为各大公司的首选。基于图注意力网络,构造分子图作为分子结构特征的药物ADMET分类预测模型进行药物研发的虚拟筛选,据有良好的精准性 [3] 。采用Chemoffice 2004中的MOPAC-PM3算法筛选量化吡喃酮类化合物的量子化学结构,利用人工神经网络中的径向基网络建立分子结构描述符与生物活性间的相关模型,有效的提高了对吡喃酮类化合物结构的预测精度 [4] 。基于RegNet-1d模型和积分梯度法的ERα拮抗剂的生物活性预测方法,通过搭建RegNet-1d深度学习模型,并以积分梯度法为理论基础进行数据结构优化,变量对生物活性影响的相关性分布,以此筛选合适的分子描述符变量,时优化后的模型预测准确率略有下降但所需测量的数据量大大减少,节约了药物研发的时间和成本 [5] 。采用分子描述、支持向量机、遗传算法三种机器学习建立ADMET的QSAR预测模型,验证结果得出可推广应用至药物代谢、毒性评估等方面 [6] 。
本文采用“华为杯”第十八届中国研究生数学建模竞赛D题中的数据,包括1974个化合物样本,每个样本都有729个分子描述符变量,5个ADMET性质数据。将从分子描述符出发构建预测模型,基于融合遗传算法的随机森林算法来预测化合物的IC50值和对应的pIC50值预测值。同时再借助SVM与Adaboost二分类模型,对化合物Caco-2、CYP3A4、hERG、HOB、MN的5中成分进行分别预测,建立ADMET分类预测模型,从而能找到既能满足较高的化合物活性,也能拥有较好的ADMET性质,助于抗乳腺癌药物的研发。
2. 构建ERα生物活性的定量预测模型
2.1. 特征选择
在不破坏原始数据可解释性的前提条件下,依据对分子描述符的显著性影响排名,进行特征选择。本文对各变量相关性进行分析,在最大程度上保留原始数据信息的同时,将分子变量描述符数据的维度从729维降至20维。在降维和筛选变量之前首先需要对全部729个变量数据进行整定。第一步:为减少计算量,删除冗余分子描述符,即需要过滤掉方差为0的特征 [7] 。第二步:对变量数据进行归一化处理使各个特征的尺度控制在相同的范围内,这样可以便于在计算分子描述符之间相关性。第三步:对以上预处理后的数据进行特征选择,将分子变量描述符数据的维度从729维降至20维。
第一步的整定算法流程,利用Python的sklearn包中Variance Threshold方法,他是一个简单的特征选择基准方法,该方法就是去除所有没有达到指定阈值的特征。默认是去除所有零方差的数据。第二步计算方法为:Min-Max标准化(Min-Max Normalization) (线性函数归一化)也称为离差标准化,是对原始数据的线性变换,使得结果映射到0~1之间。利用第一步的剔除后变量数据导入第二步进行变量归一化处理,根据上述数据整定算法可以实现对原始数据效果更佳的筛选 [8] ,采用Python语言编程实现,得到数据整定结果,整定前后样本数皆为1974,而操作变量数由729降为504。第三步利用过滤式中的互信息法 [9] 以及lasso回归算法进行特征选择 [10] ,设计了综合筛选模型,并调用Pandas工具包对各变量之间的相关度分析,去除相关度低的部分变量,最终得到符合条件的主要影响变量。采用Python语言对上述两种特征选择方法实现,得出前20个对生物活性最具有显著影响的分子描述符(即变量),maxHsOH','BCUTc-1h','minHBa','minwHBa','SaaCH','MLFER_A','maxHBa','MAXDN2','BCUTp-1h','gmin','maxHBd','maxHCsats','minssO','hmin','minHBint10','MDEO-12','minHBint6','ATSc4','MDEC-22','C2SP2',构建这些化合物对ERα生物活性的定量预测模型。
2.2. 基于改进的随机森林算法对ERα的活性预测
随机森林算法属于bagging算法的一种,也属于bagging算法的一种加强算法,是将多棵CART回归树集成的一种有监督学习算法,其样本的数据集输入为式x为SMILES,y为输出的活性值,具体公式:
(2-1)
迭代次数为t次,即是对训练集进行
次分别采样,得到最终的集合E,所得集合的算术平均值就是最后的模型输出。
在对模型进行训练时采取7:3的训练验证集比例,在训练过程中融入遗传算法加快模型收敛速度,寻找更优解,遗传算法是通过选择、交叉以及变异等机制,模拟出一个人工种群的进化过程,借鉴生物界自然选择和自然遗传机制的随机优化搜索算法,为避免决策树陷入过拟合,本文使用遗传算法来优化参数本文使用的决策树共有五个超参数:n_estimators (随机森林包含的弱分类器数量)、max_depth (树的最大深度)、min_samples_leaf (叶子节点包含的样本数)、min_samples_split (分枝所包含的最少样本个数)、max_features (选的最大特征数)。将初始种群的是数量设置为100,维度设为5,随机初始化种群,适应度函数为MSE。经过50次迭代后n_estimators = 15,max_depth = 15,min_samples_leaf = 2,min_samples_split = 3,max_features = 7为最优解。当使用原始随机森林算法时得到的准确度(R2)和均方误差(MSE)分别为0.61,0.76,利用改进后的随机森林算法获得的结果为0.73,0.52其中准确度提高了0.12,均方差减小了0.24,由实验结果可得当使用改进后的随机森林算法拟合原数据能力更强。通过改进的随机森林算法预测模型对特征1974个化合物中30%数据的IC50值和对应的pIC50值验证效果如下图1所示。
Figure 1. Validation of IC50 values and corresponding pIC50 values
图1. IC50值和对应的pIC50值验证效果
运用此模型,再对50个化合物的生物活性值进行IC50值和对应的pIC50值预测,具体预测结果如下表1所示。
Table 1. Prediction results of IC50 values and corresponding pIC50 values
表1. IC50值和对应的pIC50值预测结果
3. 构建ADMET性质预测模型
对于所提供的1974个化合物的ADMET是数据,分别构建化合物Caco-2 (小肠上皮细胞渗透性)、CYP3A4 (能否够被CYP3A4代谢)、hERG (是否具有心脏毒性)、HOB (口服生物利用度)、MN(是否具有遗传毒性)的分类预测模型,然后使用所构建的5个分类预测模型,鉴于化合物标签均为2分类模型,因此我们使用广泛且经典SVM分类模型和Adaboost分类模型进行分析比较,并通过多个评价指标来选择该化合物相应的模型。
3.1. SVM分类模型
由于数据样本是离散的,因此采用机器学习中的支持向量机(SVM)算法进行网络训练,空心圆点和黑心圆点代表两类不同类别的样本; H为分类线,H1,H2分别为平行于分类线的直线,它们经过离分类线最近的那些少量的样本点,两者间距离称为分类间隔。H线将两个不同的类正确隔离,同时使分类间隔最大化。设样本集为
,并满足:
该分类的间隔等于
,其间隔最大等价于求
的值最小。满足上式且
最小的分类平面叫做最优分类面,H1,H2
两条平行直线上的那些训练样本点称为支持向量,利用拉格朗日方法可以把求解最优分类面的原始问题
转化为其对偶问题,即在条件
下对ai求解以下最大的函数值如下:
(3-1)
ai为原问题中与(1)式对应的拉格朗日乘子,该问题就是求解二次凸规划的优化问题,存在唯一最优解,可以证明,有些乘子ai不为零,它也就是支持向量。求解该问题得到最优平面的w*和b*,此时最优分类函数
(3-2)
求和只对支持向量进行;b*是偏移量。根据泛函分析中的度量空间理论,倘若有一种核函数
满足Mercer条件的核函数替换线性算法的内积就可以找到原输入空间中对应的非线性算法。如果用特征空间的
代替x,则上式转化为下式
(3-3)
而相应的分类函数变为
(3-4)
SVM分类算法中,如果定义不同的内积函数,就能实现多项式逼近、贝叶斯分类器、径向基函数(RBF)方法等选用不同的核函数就可以构造不同的SVM [11] 。在本模型中我们选择径向基函数(RBF)为模型的核核函数构造SVM分类器。
3.2. Adaboost分类模型
在集成学习原理中,集成学习按照个体学习器之间是否存在依赖关系可以分为两类,第一个是个体学习器之间存在强依赖关系,另一类是个体学习器之间不存在强依赖关系。前者的代表算法就是是boosting系列算法。在boosting系列算法中,Adaboost是最著名的算法之一。Adaboost既可以用作分类,也可以用作回归,本次二分类任务中使用了其分类功能 [12] 。AdaBoost二元分类问题算法流程如下。
输入为样本集
,输出为{−1, +1},弱分类器算法,弱分类器迭代次数K。输出为最终的强分类器
。
1) 初始化样本集权重函数如下
(3-5)
2) 对于
:
a) 使用具有权重Dk的样本集来训练数据,得到弱分类器Gk(x)
b) 计算
的分类误差率
(3-6)
c) 计算弱分类器的系数
(3-7)
d) 更新本集的权重分布
(3-8)
这里Zk是规范化因子
(3-9)
3) 构建最终分类器为式:
(3-10)
对于Adaboost多元分类算法,其实原理和二元分类类似,最主要区别在弱分类器的系数上。比如Adaboost SAMME算法,它的弱分类器的系数如下所示
(3-11)
3.3. 测试结果及分析
化合物Caco-2、CYP3A4、hERG、HOB、MN的分类预测模型通过建立的SVM模型和Adaboost模型,使用7:3的比列划分训练集和测试集,如图2、图3所示。
如上所示图混淆矩阵主对角线为测试集中预测正确的样本,副对角线为测试集中预测错误的样本,通过实验对比各模型在化合物Caco-2、CYP3A4、hERG、HOB、MN的分类结果,以准确率,精度,召回率,F1值5种评价指标已经药物的实际作用综合选择模型。实验结果如下表2、表3所示。
所以在选择模型时不能只看某一方面数据,要做考虑到实际药物摄入与人体的友好性相关,所以在选择检测化合物Caco-2模型时,为缩小样本空间,确保药物的精度选择Adaboost模型,其精准率较SVM高0.02%。在选择化合物CYP3A4时,考虑到该药物的代谢情况,优先选择精准率较高的SVM分类器。在选择化合物hERG时,考虑到该化合物对心脏有毒性,在精度相同时,选择F1值较大的,应为F1综
Table 2. Evaluation of the Adaboost model’s prediction results
表2. Adaboost模型预测结果评估
Table 3. Evaluation of the SVM model’s prediction results
表3. SVM模型预测结果评估
合考量了正确率与召回率。及在测试该化合物时选取SVM模型。在选择化合物化合物HOB时,考虑到该化合物口服生物利用度,选择精度较高的Adaboost模型。在选择化合物MN时,考虑到该化合物的遗传毒性,选择精度较高的Adaboost模型。
4. 秃鹰搜索(BES)与化合物筛选模型的建立
为了能从操作变量的取值区间内找到对抑制ERα具有更好的生物活性,同时具有更好的ADMET性质的化合物,将使用一种智能搜索算法,将模型1作为种群的适应度函数,通过迭代计算使种群向最佳适应度位置(药物活性最高时的20操作变量的值)移动,直至完成求解过程。其次,从智能搜索算法的历史最佳搜索位置中依次选取全局最优位置,将此位置输入模型2至模型6完成ADMET性质预测,直至筛选出ADMET性质满足要求的药物。
秃鹰搜索算法:本文模型一使用改进的随机森林建模,其内部包含多颗决策树作为弱分类器,由于决策树工作时呈现非线性,使得模型一总体上呈现较强的非线性 [13] 。在模型一中搜索最优解是一种复杂的非线性过程,而在这类问题上具有显著的竞争优势,因此选择这一算法来进行建模。
BES寻优过程可分为三个阶段。在第一阶段(选择空间),鹰选择猎物数量最多的空间,从中心点移动到选定的搜索区域,具有较强全局寻优能力;在第二阶段(空间搜索),鹰在选定的空间内移动寻找猎物,具有较强全局寻优能力;在第三阶段(俯冲),鹰从第二阶段确定的最佳位置摆动,并确定最佳狩猎,即全局最优位置。计算如下:
1) 选择阶段:
秃鹰在选择的搜索空间中识别并选择最佳的区域(根据食物的数量),在那里它们可以捕食。公式(4-1)用数学方法表示了这种行为:
(4-1)
其中,这里是控制位置变化的参数,取值在[1.5, 2]之间。是[0, 1]之间的一个随机数。Pbest是当前最优位置,Pmean表示利用前面所有点的信息。Pi为第只秃鹰的位置。
2) 搜索阶段:
秃鹰在选定的搜索空间内搜索猎物,并在螺旋形空间内以不同的方向移动,以加快搜索速度。俯冲的最佳位置用以下数学公式表示
(4-2)
(4-3)
(4-4)
(4-5)
其中:(i)与(i)分别为螺旋方程的极角与极径,a与R是控制螺旋轨迹的参数,变化范围分别为(5, 10)、(0.5, 2)。rand为(0, 1)内随机数,x(i)与y(i)表示极坐标中秃鹰位置,取值均为(−1, 1)。秃鹰位置跟新如以下公式所示:
(4-6)
3) 俯冲阶段:
秃鹰从搜索空间的最佳位置摇摆到它们的目标猎物。所有的点也都朝着最好的点移动。下面三个公式从数学上说明了这种行为:式中:增加了秃鹰向最佳点和中心点的移动强度.
(4-7)
(4-8)
(4-9)
俯冲过程,位置根据以下式格式跟新如下,其中c1与c2表示秃鹰向最佳与中心位置的运动强度,取值均为(1, 2)。
(4-10)
(4-11)
药物筛选策略:药物筛选策略分为两部分,第一部分利用黑鹰搜索算法从模型一寻找全局最优位置,即最高pIC50值所对应的20个操作变量的值。第二部分用模型2~6对第一部分中得到的全局最优位解进行ADMET性质友好性筛选,挑选符合条件的解。
① 黑鹰搜索寻找模型一全局最优位置
设置初始变量的维度为20,正好与选取的20个操作变量一一对应,使用模型1的预测结果作为种群适应度计算。种群大小设置为100,迭代次数50,算法中的超参数使用文献 [14] 提供的值。算法首先从0到1之间随机选取操作变量初始值,接着算法进入迭代求解,每次迭代都经历选择空间、空间搜索、俯冲三个阶段,直到达到最大迭代次数,或满足求解精度,黑鹰搜索将每次迭代的结果存放至历史最优解列表,它代表者黑鹰搜索寻优的轨迹,常用收敛曲线表示。随着迭代次数的增加,得到的最优解逐渐接近真实最优解,并且最终收敛在真实最优解附近。
② ADMET性质友好性筛选
ADMET性质友好性主要是指:满足化合物的小肠上皮细胞渗透性较好、能够被CYP3A4代谢、不具有心脏毒性、口服生物利用度较好、不具有遗传毒性,五者中的三个及以上指标 [15] 。也就是模型2~6预测所预测的结果中至少有满足以下条件中的至少三条:Caco-2为‘1’、CYP3A4为‘1’、hERG为‘0’、HOB为‘1’、MN为‘0’。
从黑鹰搜索的历史最优解中从后往前依此取出一个最优解,并将最优解分别输入至模型2~6,五个二分类模型根据输入值预测Caco-2、CYP3A4、hERG、HOB、MN。
并且使用以下函数作为评分器,进行打分。其中Modeli(x)表示模型i对化合物x的预测结果,target = {1, 1, 0, 1, 0},评分器表示输入的化合物满足ADMET性质友好性中的数量。当评分器输出结果大于或等于3时,认为该化合物DMET性质友好,将其作为抗胰腺癌候选药物输出。如果该值低于3,则放弃放弃该药物,并从黑鹰搜索的历史最优解中取出下一种药物(诺该化合物与前面值相同,则自动往后取出下一种化合物)如下式,判断器ADMET性质对人体友好性。
(4-12)
5. 结果分析
运行上述程序次得到如下表4所示结果如下,其中pIC50为9.113,活性较高。其ADMET性质中Caco-2为1,CYP3A4为0,hERG为0,即该化合物对小肠上皮细胞渗透性较好、能够被CYP3A4代谢、不具有心脏毒性。虽然该化合物MN为1、HOB为0,但总体上表现为ADMET友好性。
Table 4. Properties of candidate compounds
表4. 候选化合物性质
表5是黑鹰搜索的候选化合物原始坐标,由于原数据经过为归一化处理,因此需要将获得的候选物坐标还原成真实值。如表6所示,其值满足所选20个操作变量的取值范围,连续性等,故该化和物能够被选为真实候选药物。因此当选取的20个操作变量的值位于该化合物附近时,可认为其具有与该化合物相似的性质。
Table 5. Real coordinates of candidate compounds
表5. 候选化合物真实坐标
Table 6. Results of Black Hawk search
表6. 黑鹰搜索结果
根据小白鼠试验发现ERα活性的化合物是治疗乳腺癌的候选药物,但要想成为药物,就需要兼顾ADMET性质,因此本文首先采取基于改进的随机森林算法对ERα的活性预测建立模型一,再通过SVM分类模型和Adaboost分类模型构建化合物Caco-2 (小肠上皮细胞渗透性)、CYP3A4 (能否够被CYP3A4代谢)、hERG (是否具有心脏毒性)、HOB (口服生物利用度)、MN (是否具有遗传毒性)的分类预测模型建立模型二,以上两种模型均使用7:3的比列划分训练集和测试集。最后通过秃鹰搜索算法将前两个模型结合建立最终的药物筛选模型,从而能得到最优解,实验结果表明利用此方法能对备筛化合物库进行预筛,去除ADMET、分子稳定性、溶解性等性质较差的化合物,通过减少高通量筛选的压力,从而提高筛选效率;同时数据库质量的提升还可以降低药物开发后期的失败率,以避免化合物药代动力学引起的开发失败,降低开支。
NOTES
*通讯作者。