抗乳腺癌候选药物的生物活性预测建模

Predictive Modeling of Bioactivity of Anti-Breast Cancer Drug Candidates

摘要: 本文构建化合物对Erα生物活性的定量预测模型，结合集成学习方法与逻辑回归方法，使用自动的参数优化方法，使各算法达到最优泛化性能。首先使用随机森林算法，以信息理论为基础，将化合物的分子描述符对雌激素受体α亚型的活性影响进行特征重要性排序，得到可用于算法判断的20个高效变量；再根据这20个高效分子描述符，利用岭回归算法实现对ERα生物活性的定量预测。结果表明，该模型可以准确预测Erα生物活性，为科学选择抗乳腺癌药物提供了新思路。

Abstract: In this paper, a quantitative prediction model for ERα biological activity of compounds was con-structed, combined with integrated learning method and logistic regression method, and automatic parameter optimization method was used to achieve the optimal generalization performance of each algorithm. First, based on information theory, random forest algorithm was used to rank the characteristic importance of the effects of molecular descriptors of compounds on the activity of es-trogen receptor α subtypes, and 20 efficient variables were obtained. Based on these 20 molecular descriptors, ridge regression algorithm was used to quantitatively predict the biological activity of ERα. The results show that this model can accurately predict the biological activity of Erα, which provides a new idea for scientific selection of anti-breast cancer drugs.

1. 引言

乳腺癌是目前世界上最常见，致死率较高的癌症之一 [1] 。乳腺癌的发展与雌激素受体密切相关。有研究发现 [2] ，雌激素受体α亚型(Estrogen receptors alpha, ERα)在不超过10%的正常乳腺上皮细胞中表达，但大约在50%~80%的乳腺肿瘤细胞中表达；而对ERα基因缺失小鼠的实验结果表明，ERα在乳腺发育过程中扮演了十分重要的角色 [3] 。目前，抗激素治疗常用于ERα表达的乳腺癌患者，其通过调节雌激素受体活性来控制体内雌激素水平。ERα被认为是治疗乳腺癌的重要靶标，能够拮抗ERα活性的化合物可能是治疗乳腺癌的候选药物 [4] 。比如，临床治疗乳腺癌的经典药物他莫昔芬和雷诺昔芬就是ERα拮抗剂。而IC₅₀和pIC₅₀作为影响ERα活性的重要指标，准确预测IC₅₀和pIC₅₀的值对药物的选取具有重要指导意义。因此，构建化合物对Erα生物活性的定量预测模型，对科学选取抗乳腺癌药物提供了新思路。

2. 变量对生物活性影响的重要性分析

2.1. 模型选取

本文针对1974个化合物的729个分子的数据进行预先处理，通过数据清洗发现存在与目标值无相关性的全零列，判断出来共有225列全零列。去除全零列可以降低模型的复杂度，提高模型的稳健性。

随机森林算法是最常用也是最强大的监督学习算法之一，它兼顾了解决回归问题和分类问题的能力。随机森林是通过集成学习的思想，将多棵决策树进行集成的算法，能够提供更好的预测性能 [5] [6] 。许多研究表明 [7] ，RF算法对噪声、过拟合和异常值具有较高的预测精度。此外，与SVR等其他算法相比，该算法也具有许多优点，如处理高维度的复杂数据和具有更少的参数。给定一个数据集：

$D=\left\{\left(X_i,y_i\right)\right\}\left(\left|D\right|\right)=n,X_i\in R^m,y_i$ (1)

式中：n为数据集样本的数量；m为每个样本的特征个数。将数据集划分为训练集和测试集，其中训练集为：

$S_L=\left\{\left(X_1,y_1\right),\left(X_2,y_2\right),\cdots ,\left(X_L,y_L\right)\right\}$ (2)

式中：L表示训练集的数量。

2.2. RF参数优化

随机森林算法把数据分为训练集(Trading Set)、验证集(Validation Set)、测试集(Test Set)，每次随机选出n组数据，用训练集训练出n个模型，测试集对n个模型进行评价，选出最终模型 [7] 。

在确定决策树算法的划分依据(criterion)时，可选择基尼系数(Gini)、信息熵(Entropy)、信息增益(Gain)等参数。为使模型获得更好的参数，获得最佳泛化性能，对模型进行预剪枝处理，信息熵计算公式表示为：

$\text{Entropy}\left(D\right)=-{\displaystyle \underset{k=1}{\overset{K}{\sum }}\frac{\left|C_k\right|}{\left|D\right|}}{\log }_2\frac{\left|C_k\right|}{\left|D\right|}$ (3)

式中， $\text{Entropy}\left(D\right)$ 表示信息熵，k为类别总数，D为样本总数， $C_k$ 为属于某个类别的样本数。

为使评估模型更加准确可信，在模型中使用网格搜索和交叉验证(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) 将测试集带入训练好的决策树模型中，其预测值为： $\text{M1}\left(\text{X}\right),\text{M2}\left(\text{X}\right),\text{M3}\left(\text{X}\right),\cdots ,\text{Mk}\left(\text{X}\right)$ 。

4) 将所有决策树的预测值的众数作为随机森林模型的预测结果。

2.4. 重要性分析

通过上述四个步骤，建立模型；根据变量对生物活性影响的重要性进行排序，根据模型计算的重要性结果排序得到前20个对生物活性最具有显著影响的分子描述符(即变量)分别别是；LipoaffinityIndex、MDEC-23……MlogP，对生物活性最具有显著影响的前20个分子的重要度排序如图1所示。

Figure 1. The top 20 descriptions of molecules that have the most significant effects on biological activity

图1. 前20个对生物活性最具有显著影响的分子描述

3. IC₅₀值及pIC₅₀值预测

3.1. 岭回归算法及优化算法

岭回归是一种改良的最小二乘估计法，通过放弃最小二乘法的无偏性，以损失部分信息、降低精度为代价获得回归系数，它是更为符合实际、更可靠的回归方法，对存在离群点的数据的拟合要强于最小二乘法 [11] 。不同与线性回归的无偏估计，岭回归的优势在于它的无偏估计，更趋向于将部分系数向0收缩。因此，它可以缓解多重共线问题，以及过拟合问题。

在标准线性回归中，通过最小化真实值 $y_i$ 和预测值 ${\hat{y}}_i$ 的平方误差来训练模型，这个平方误差值也被成为残差平方和(RSS, Residual Sum of Squares)：

$\text{RSS}={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(y_i-{\hat{y}}_i\right)}^2}$ (4)

最小二乘法即最小化残差平方和，为

$J_{\beta }\left(\beta \right)=\arg \min {\displaystyle \underset{i=1}{\overset{p}{\sum }}{\left(y_i-x_i{\beta }_i-{\beta }_0\right)}^2}$ (5)

将其转化为矩阵形式：

$J_{\beta }\left(\beta \right)=\arg \min {\left(y-X\beta \right)}^{\text{T}}\left(y-X\beta \right)$ (6)

求解为：

$\beta =\left(X^{\text{T}}X\right)X^{\text{T}}y$ (7)

将不适定问题转化为适定问题，在矩阵 $\left(X^{\text{T}}X\right)$ 的对角线元素上加入一个小的常数值 $\lambda$ ，然后取其逆得系数：

${\hat{\beta }}_{ridge}={\left(X^{\text{T}}X+\lambda I_n\right)}^{-1}\left(X^{\text{T}}Y\right)$ (8)

$I_n$ 是单位矩阵，对角线全是1，类似于“山岭”， $\lambda$ 是岭函数，改变其数值可以改变单位矩阵对角线的值。

随后，代价函数 $J_{\beta }\left(\beta \right)$ 在RSS的基础上加入了对系数值的惩罚，矩阵形式为：

$J_{\beta }\left(\beta \right)={\displaystyle \underset{i=1}{\overset{p}{\sum }}{\left(y_i-X\beta \right)}^2}+\lambda {\displaystyle \underset{j=0}{\overset{n}{\sum }}{\beta }_i^2}={\displaystyle \underset{i=1}{\overset{p}{\sum }}{\left(y_i-X\beta \right)}^2}+\lambda {\Vert \beta \Vert }^2$ (9)

损失函数：

$J\left(a\right)={\left(h\left(x_1\right)-y_1\right)}^2+{\left(h\left(x_2\right)-y_2\right)}^2+{\left(h\left(x_3\right)-y_3\right)}^2+\cdots +{\left(h\left(x_m\right)-y_m\right)}^2$ (10)

式中， $h\left(x_1\right)$ 为第i个训练样本特征值组合预测函数， $y_i$ 为第i个训练样本的真实值。

通过梯度下降的方法 [12] ，以损失函数限制，优化模型精度，进行迭代求解，需要手动指定超参数，迭代求得最优参数的过程如图2所示。

$w_1=w_1-\alpha \frac{\vartheta \cos t\left(w_0+w_1{x}_1\right)}{\vartheta x_1}$ (11)

$w_0=w_0-\alpha \frac{\vartheta cost\left(w_0+w_1{x}_1\right)}{\vartheta x_1}$ (12)

$\alpha$ 为学习速率， $\alpha$ 旁边的整体表示方向沿着这个函数下降的方向线，最后就能找到山谷的最低点，然后更新w值使用，面对训练数据规模十分庞大的任务，能够找到较好的结果。

Figure 2. Gradient descent method

图2. 梯度下降法

回归性能评估：

$\text{MSE}=\frac{1}{m}{\displaystyle \underset{i=1}{\overset{m}{\sum }}{\left(y^i-\bar{y}\right)}^2}$ (13)

式中， $y^i$ 为预测值， $\bar{y}$ 为真实值。

3.2. 模型训练

使用第二节提取的20个分子特征和及分子pIC₅₀ (相比于IC₅₀_nM，pIC₅₀的数量级相差较小，更适合于回归模型的建立)为数据集，对岭回归模型进行训练。

岭回归算法中正则化力度越大，惩罚项占据主导地位，使得每个自变量的权重系数趋近于零；正则化力度越小，惩罚项的影响也越来越小，导致每个自变量的权重系数震荡的幅度变大，如图3所示。

Figure 3. Relationship between regularization force and weight coefficient

图3. 正则化力度与权重系数关系

另外，对迭代次数、回归偏置设置(模型结果更具稳健性)、是否进行标准化处理(进一步减小因噪声造成的误差)、优化器选择(在数据量增大的时候可以选择随机平均梯度法，进一步对收敛速度进行优化)设置见表1。

Table 1. Predictor parameter settings

表1. 预估器参数设置

3.3. 模型结果

将化合物的分子描述符与pIC₅₀指标建立映射关系，以80%的数据作为训练集，20%的数据作为测试集。得到岭回归输出模型，保留小数点后三位有效数字后的回归模型如式所示。20个分子描述(根据第二节所得)的权重系数如图4所示，(仅标注大于等于1%的权重系数)，偏置为4.924，计算模型的均方误差仅为0.808。

$\begin{array}{c}h\left(x\right)=w^{\text{T}}x+b\\ =-1.367x_1+1.238x_2-3.566x_3-1.298x_4-4.362x_5-2.997x_6+1.361x_7+1.909x_8\\ +2.439x_9+4.353x_{10}3.054x_{11}+1.664x_{12}-1.904x_{13}+2.392x_{14}-4.864x_{15}+3.559x_{16}\\ -3.891x_{17}+1.123x_{18}-8.451x_{19}+2.145x_{20}+4.924\end{array}$ (14)

Figure 4. Molecular description weight coefficient

图4. 分子描述权重系数

3.4. 模型验证

通过训练好的模型，对已有50个化合物进行pIC₅₀值和对应的IC₅₀值进行计算(表2)。从图5中可以看到，预测值可随样本变化迅速并准确地靠近真实值，并且从图6和图7中可以看到，两个图的图像走向是对应的，这也印证了所建立的模型的准确性。

根据训练好的模型，对已有50个化合物进行IC₅₀值得预测，根据得到得pIC₅₀得值，按数学关系对IC₅₀值进行求解：

${\text{IC}}_{50}={10}^{\left(9-{\text{pIC}}_{50}\right)}$ (15)

Figure 5. The comparison of the predicted value and the true value of pIC₅₀

图5. pIC₅₀预测值与真实值对比

Figure 6. The predicted value of pIC₅₀

图6. pIC₅₀值预测值

Figure 7. The calculated value of IC₅₀_nM

图7. IC₅₀_nM计算值

Table 2. Predicted values of IC50_nM and pIC50

表2. IC₅₀_nM和pIC₅₀预测值

4. 结论

本文采用了随机森林算法，以信息理论为基础，将化合物的分子描述符对雌激素受体α亚型的活性影响进行特征重要性排序，得到可用于算法判断的20个高效分子描述符，使用岭回归算法实现对ERα生物活性的定量预测，选择线性度较好的pIC₅₀指标作为算法训练的目标值，相较于IC₅₀，可获得更稳定的求解结果，所用建模方法在具有目标值的数据集上可以获得较好的表现，为科学选择抗乳腺癌药物提供了新思路。