1. 引言
乳腺癌是目前世界上最常见,致死率较高的癌症之一 [1] 。乳腺癌的发展与雌激素受体密切相关。有研究发现 [2] ,雌激素受体α亚型(Estrogen receptors alpha, ERα)在不超过10%的正常乳腺上皮细胞中表达,但大约在50%~80%的乳腺肿瘤细胞中表达;而对ERα基因缺失小鼠的实验结果表明,ERα在乳腺发育过程中扮演了十分重要的角色 [3] 。目前,抗激素治疗常用于ERα表达的乳腺癌患者,其通过调节雌激素受体活性来控制体内雌激素水平。ERα被认为是治疗乳腺癌的重要靶标,能够拮抗ERα活性的化合物可能是治疗乳腺癌的候选药物 [4] 。比如,临床治疗乳腺癌的经典药物他莫昔芬和雷诺昔芬就是ERα拮抗剂。而IC50和pIC50作为影响ERα活性的重要指标,准确预测IC50和pIC50的值对药物的选取具有重要指导意义。因此,构建化合物对Erα生物活性的定量预测模型,对科学选取抗乳腺癌药物提供了新思路。
2. 变量对生物活性影响的重要性分析
2.1. 模型选取
本文针对1974个化合物的729个分子的数据进行预先处理,通过数据清洗发现存在与目标值无相关性的全零列,判断出来共有225列全零列。去除全零列可以降低模型的复杂度,提高模型的稳健性。
随机森林算法是最常用也是最强大的监督学习算法之一,它兼顾了解决回归问题和分类问题的能力。随机森林是通过集成学习的思想,将多棵决策树进行集成的算法,能够提供更好的预测性能 [5] [6] 。许多研究表明 [7] ,RF算法对噪声、过拟合和异常值具有较高的预测精度。此外,与SVR等其他算法相比,该算法也具有许多优点,如处理高维度的复杂数据和具有更少的参数。给定一个数据集:
(1)
式中:n为数据集样本的数量;m为每个样本的特征个数。将数据集划分为训练集和测试集,其中训练集为:
(2)
式中:L表示训练集的数量。
2.2. RF参数优化
随机森林算法把数据分为训练集(Trading Set)、验证集(Validation Set)、测试集(Test Set),每次随机选出n组数据,用训练集训练出n个模型,测试集对n个模型进行评价,选出最终模型 [7] 。
在确定决策树算法的划分依据(criterion)时,可选择基尼系数(Gini)、信息熵(Entropy)、信息增益(Gain)等参数。为使模型获得更好的参数,获得最佳泛化性能,对模型进行预剪枝处理,信息熵计算公式表示为:
(3)
式中,
表示信息熵,k为类别总数,D为样本总数,
为属于某个类别的样本数。
为使评估模型更加准确可信,在模型中使用网格搜索和交叉验证(GS-CV)的方式对参数进行优化 [8] ,GS-CV首先遍历给定的参数组合来优化模型表现。进而将训练集继续拆分成训练集和测试集,对模型进行交叉验证。最终确定以entropy为划分依据,并限制决策数的最大深度max_depth = 6,决策树数量为n_estimators = 40,以防止过拟合现象。
2.3. 模型训练
采取随机有放回抽样,抽取70%的数据作为训练集,30%的数据作为测试集。设置随机数种子为0以复现算法。
1) 利用Bootstrap方法对训练集随机采样 [9] ,得到各子训练样本集合和测试样本集合。
2) 利用步骤1)得到的新的训练集分别建立多个基模型。在每棵树的节点处,先从m个特征中随机选取t个特征,使用这些特征中的最好的分裂方式对每一个节点进行分裂,使用的分裂方法是CART [10] 。
3) 将测试集带入训练好的决策树模型中,其预测值为:
。
4) 将所有决策树的预测值的众数作为随机森林模型的预测结果。
2.4. 重要性分析
通过上述四个步骤,建立模型;根据变量对生物活性影响的重要性进行排序,根据模型计算的重要性结果排序得到前20个对生物活性最具有显著影响的分子描述符(即变量)分别别是;LipoaffinityIndex、MDEC-23……MlogP,对生物活性最具有显著影响的前20个分子的重要度排序如图1所示。
![](//html.hanspub.org/file/10-2570826x13_hanspub.png?20230509100453479)
Figure 1. The top 20 descriptions of molecules that have the most significant effects on biological activity
图1. 前20个对生物活性最具有显著影响的分子描述
3. IC50值及pIC50值预测
3.1. 岭回归算法及优化算法
岭回归是一种改良的最小二乘估计法,通过放弃最小二乘法的无偏性,以损失部分信息、降低精度为代价获得回归系数,它是更为符合实际、更可靠的回归方法,对存在离群点的数据的拟合要强于最小二乘法 [11] 。不同与线性回归的无偏估计,岭回归的优势在于它的无偏估计,更趋向于将部分系数向0收缩。因此,它可以缓解多重共线问题,以及过拟合问题。
在标准线性回归中,通过最小化真实值
和预测值
的平方误差来训练模型,这个平方误差值也被成为残差平方和(RSS, Residual Sum of Squares):
(4)
最小二乘法即最小化残差平方和,为
(5)
将其转化为矩阵形式:
(6)
求解为:
(7)
将不适定问题转化为适定问题,在矩阵
的对角线元素上加入一个小的常数值
,然后取其逆得系数:
(8)
是单位矩阵,对角线全是1,类似于“山岭”,
是岭函数,改变其数值可以改变单位矩阵对角线的值。
随后,代价函数
在RSS的基础上加入了对系数值的惩罚,矩阵形式为:
(9)
损失函数:
(10)
式中,
为第i个训练样本特征值组合预测函数,
为第i个训练样本的真实值。
通过梯度下降的方法 [12] ,以损失函数限制,优化模型精度,进行迭代求解,需要手动指定超参数,迭代求得最优参数的过程如图2所示。
(11)
(12)
为学习速率,
旁边的整体表示方向沿着这个函数下降的方向线,最后就能找到山谷的最低点,然后更新w值使用,面对训练数据规模十分庞大的任务,能够找到较好的结果。
回归性能评估:
(13)
式中,
为预测值,
为真实值。
3.2. 模型训练
使用第二节提取的20个分子特征和及分子pIC50 (相比于IC50_nM,pIC50的数量级相差较小,更适合于回归模型的建立)为数据集,对岭回归模型进行训练。
岭回归算法中正则化力度越大,惩罚项占据主导地位,使得每个自变量的权重系数趋近于零;正则化力度越小,惩罚项的影响也越来越小,导致每个自变量的权重系数震荡的幅度变大,如图3所示。
![](//html.hanspub.org/file/10-2570826x38_hanspub.png?20230509100453479)
Figure 3. Relationship between regularization force and weight coefficient
图3. 正则化力度与权重系数关系
另外,对迭代次数、回归偏置设置(模型结果更具稳健性)、是否进行标准化处理(进一步减小因噪声造成的误差)、优化器选择(在数据量增大的时候可以选择随机平均梯度法,进一步对收敛速度进行优化)设置见表1。
![](Images/Table_Tmp.jpg)
Table 1. Predictor parameter settings
表1. 预估器参数设置
3.3. 模型结果
将化合物的分子描述符与pIC50指标建立映射关系,以80%的数据作为训练集,20%的数据作为测试集。得到岭回归输出模型,保留小数点后三位有效数字后的回归模型如式所示。20个分子描述(根据第二节所得)的权重系数如图4所示,(仅标注大于等于1%的权重系数),偏置为4.924,计算模型的均方误差仅为0.808。
(14)
![](//html.hanspub.org/file/10-2570826x40_hanspub.png?20230509100453479)
Figure 4. Molecular description weight coefficient
图4. 分子描述权重系数
3.4. 模型验证
通过训练好的模型,对已有50个化合物进行pIC50值和对应的IC50值进行计算(表2)。从图5中可以看到,预测值可随样本变化迅速并准确地靠近真实值,并且从图6和图7中可以看到,两个图的图像走向是对应的,这也印证了所建立的模型的准确性。
根据训练好的模型,对已有50个化合物进行IC50值得预测,根据得到得pIC50得值,按数学关系对IC50值进行求解:
(15)
![](//html.hanspub.org/file/10-2570826x42_hanspub.png?20230509100453479)
Figure 5. The comparison of the predicted value and the true value of pIC50
图5. pIC50预测值与真实值对比
![](//html.hanspub.org/file/10-2570826x43_hanspub.png?20230509100453479)
Figure 6. The predicted value of pIC50
图6. pIC50值预测值
![](//html.hanspub.org/file/10-2570826x44_hanspub.png?20230509100453479)
Figure 7. The calculated value of IC50_nM
图7. IC50_nM计算值
![](Images/Table_Tmp.jpg)
Table 2. Predicted values of IC50_nM and pIC50
表2. IC50_nM和pIC50预测值
4. 结论
本文采用了随机森林算法,以信息理论为基础,将化合物的分子描述符对雌激素受体α亚型的活性影响进行特征重要性排序,得到可用于算法判断的20个高效分子描述符,使用岭回归算法实现对ERα生物活性的定量预测,选择线性度较好的pIC50指标作为算法训练的目标值,相较于IC50,可获得更稳定的求解结果,所用建模方法在具有目标值的数据集上可以获得较好的表现,为科学选择抗乳腺癌药物提供了新思路。