基于RFR模型的抗乳腺癌候选药物优化

Optimization of Candidate Drugs for Anti-Breast Cancer Based on RFR Model

摘要: 在抗乳腺癌药物研发中，为了节省时间成本，建立化合物预测模型来筛选活性化合物是一种有效的方法。本文根据提供的乳腺癌治疗靶标雌激素受体α亚型(Estrogen receptors alpha, ERα)拮抗剂信息，利用Lasso回归与随机森林相结合的方法，对数据降维并筛选出影响生物活性的主要变量，明确了模型建立的方向；在此基础上，建立了关于生物活性的随机森林回归模型并进行了预测。

Abstract: In the research and development of anti-breast cancer drugs, it is an effective method to establish a compound prediction model to screen active compounds in order to save time and cost. According to the information of Estrogen receptors alpha (ERα) antagonist, which is the target of breast cancer treatment, this paper uses Lasso regression and random forest to reduce the dimension of the data and screen out the main variables that affect the biological activity, and makes clear the direction of establishing the model. On this basis, a random forest regression model about biological activity was established and predicted.

1. 引言

根据世界卫生组织国际癌症研究署发布的《2020年全球癌症负担报告》显示：全球乳腺癌新发病例高达226万例，已经超过肺癌(221万例)成为全球第一大癌，占全球新发癌症病例的11.7%。乳腺癌是目前世界上最常见、致死率较高的癌症之一 [1]。而我国乳腺癌病症情况更加严峻。据统计，2020年我国乳腺癌新发病例约42万，乳腺癌发病数高居全球第一，死亡人数约11.7万，居女性癌症死亡首位 [2]。乳腺癌发病率之高，给我国女性带来了沉重的疾病负担。

因此，我国对乳腺癌药物的研发进程优化势在必行。研究发现，乳腺癌的发展与雌激素受体密切相关，雌激素受体α亚型(Estrogen receptors alpha, ERα)在小于等于10%的正常乳腺上皮细胞中表达，但大约在50%~80%的乳腺肿瘤细胞中表达；利用ERα基因缺失小鼠进行实验验证，发现ERα确实在乳腺发育过程中扮演了十分重要的角色。目前，抗激素治疗常用于ERα表达的乳腺癌患者，其通过调节雌激素受体活性来控制体内雌激素水平。因此，ERα被认为是治疗乳腺癌的重要靶标，能够拮抗ERα活性的化合物可能是治疗乳腺癌的候选药物。当前，在药物研发中，通过建立化合物活性预测模型的方法来筛选有效活性化合物可以加快研发进展和降低研发成本。此处对乳腺癌的具体做法是：针对相关靶标(ERα)，收集一系列作用于靶标的化合物及其生物活性数据，然后以一系列分子结构描述符作为自变量，化合物的生物活性值作为因变量，构建化合物的定量结构–活性关系模型，然后使用该模型预测具有更好生物活性的新化合物分子，或者指导已有活性化合物的结构优化。

2. 模型建立

本文首先对1974个化合物的729个分子描述符(即变量)进行变量选择，根据变量对生物活性影响的重要性进行排序。通过观察处理后的数据发现，变量数大于测量次数，也就说明当用所有的变量去表示目标值的时候，数据矩阵无法做到列满秩，某些变量可以通过其他变量进行表示，即变量与变量之间存在着多重共线性。Lasso回归在解决多重共线性的问题上具有优势。应用Lasso回归筛选出的变量能够很好地对模型进行表达，且各变量之间相互独立，因此选用Lasso回归方法进行主要变量的筛选。随机森林(RF)是利用bootstrap重抽样方法从原始样本中抽取多个样本，对每个bootstrap样本进行决策树建模，然后组合多棵决策树的预测，通过投票得出最终预测结果。它具有很高的预测准确率，并可以根据变量的重要性进行排序。具体流程如图1所示，其中x_n为样本个数， $n=1,2,\cdots ,729$ ；x_m为Lasso初步降维筛选所得自变量， $m=1,2,\cdots ,157$ ；x_r为随机森林回归二次降维筛选所得自变量， $r=1,2,\cdots ,20$ 。

2.1. Lasso回归初步降维

根据经验规则，如果方差膨胀因子VIF > 10，则认为该回归方程存在严重的多重共线性。经检验，数据材料所提供初始样本VIF = 1.55e+06，存在高度多重非线性问题。如果多个解释变量之间高度相关，则不容易区分它们各自对被解释变量的单独影响力。在应对多重共线性的问题中，Lasso回归能够客观的筛选出有效变量从而解决多重共线性问题。在进行建模之前，需要将原始数据表格中变量全部为0的参数删除。从剩余504个变量中，找出主要的操作变量，同时，这些变量之间应不具有多重共线性，它们相互独立，具有代表性。而考虑到Lasso回归收缩估计量存在偏差，因此需要对这些选择出来的参数再次进行最小二乘法回归，并以此回归系数判定这些主要变量对于目标值的实际权重，根据权重系数去除趋近于0的系数的参数，最终就可以得出对于目标值最具影响力的参数变量。

Figure 1. Flow chart

图1. 思路流程图

Lasso回归是1997年由Tibshirani提出的一种压缩估计方法，它的本质其实是正则化的最小二乘法，通过给最小二乘法添加正则化的惩罚函数 [3]，让回归系数的绝对值之和在小于一个常数的约束条件下，使得回归模型残差平方和最小，产生严格等于零的回归系数，从而判定某些变量对于目标值的有效程度，筛选出有效变量 [4]。

普通的最小二乘法线性回归对应的线性模型如下：

$f\left(x\right)=w_1{x}_1+w_2{x}_2+\cdot \cdot \cdot +w_d{x}_d+b$ (1)

用向量表示为，当用空间里的一个直观的模型取逼近这些数据点时，往往无法完全经过这些点，即模型理论值与数据值存在一定的偏差，其中b为偏差值 [5]。此时我们需要这些误差的总和最小，而由于实际值分布较为分散，往往在理论模型的上下两侧，所以若是将这些误差直接相加，则会出现正负抵消的现象，从而无法正确地估计误差的大小。因此，需要将每一项的误差项平方，即将误差项正向化，根据这个方法，得出最小二乘法的误差项公式：

$L\left(w\right)={\displaystyle \underset{i=1}{\overset{N}{\sum }}{\Vert f\left(x_i\right)-y_i\Vert }^2={\displaystyle \underset{i=1}{\overset{N}{\sum }}{\Vert w^{\text{T}}{x}_i-y_i\Vert }^2}}={\displaystyle \underset{i=1}{\overset{N}{\sum }}{\left(w^{\text{T}}{x}_i-y_i\right)}^2}=\left(w^{\text{T}}{X}^{\text{T}}-Y^{\text{T}}\right)\left(wX-Y\right)$ (2)

将上式展开得，

$L\left(w\right)=\left(w^{\text{T}}{X}^{\text{T}}-Y^{\text{T}}\right)\left(wX-Y\right)=w^{\text{T}}{X}^{\text{T}}-2w^{\text{T}}{X}^{\text{T}}Y+Y^{\text{T}}Y$ (3)

为了得出使误差项最小得w，对其求偏导，得出

$w={\left(X^{\text{T}}X\right)}^{-1}{X}^{\text{T}}Y$ (4)

但是在这种情况下，为了满足矩阵可逆的条件下，只对于目标值y的个数大于变量值x的个数，一旦变量值x的个数过多，就会引起多重共线性和过拟合的现象 [6]。

Lasso回归为了解决这样得问题引入了扰动项构建正则化框架，使得矩阵变为可逆矩阵

$\underset{w}{\arg \min }\left[L\left(w\right)+\lambda p\left(w\right)\right]$ (5)

其中， $P\left(w\right)=\Vert w\Vert$ ，此时Lasso得回归模型为

$\hat{w}\left(\text{lasso}\right)=\underset{w}{\arg \min }{\Vert y-Xw\Vert }_2^2+\lambda {\Vert w\Vert }_1$ (6)

其中 ${\Vert w\Vert }_1={\displaystyle \underset{j=1}{\overset{P}{\sum }}\left|w_j\right|}$ 为w上的1范数罚分，它可以稀释解的个数，λ是一个调整参数，选择合适的λ能够调整最小二乘法拟合，并将一些分量得系数参数收缩为0。通过这样得方法，Lasso回归就将所有的变量筛选出对于目标值影响较大的显著变量，将影响较小的变量系数直接压缩至0。

在对数据进行Lasso回归之前，由于数据变量较多，且量纲都不统一，因此，需要对数据先通过Matlab进行标准化，然后用Stata对标准化后的数据进行Lasso回归。当λ作为函数，在不同的λ约束条件下，回归系数也会发生很大的改变，当λ过大时，惩罚力度也会变大，使得所有的回归系数都归于0。

紧接着，通过使用K折交叉验证的方法来选择最佳参数。K折交叉验证，就是将样本数据随机分为K个等分。将第一个子样本作为“验证集”保留不动，其余的子集作为“训练集”来估计这个模型，再以此预测保留的第一个子样本，并计算这个子样本的均方误差。依次类推，将所有的子样本的均方误差加总，再通过调整参数，使得整个样本的均方误差最小，从而具有最佳的预测能力。

Lasso回归不仅解决了高度多重共线性问题，还进行了样本变量的初步降维，一次降维后，筛选出157个变量。

2.2. 随机森林二次降维与结果分析

在Lasso回归降维结果的基础上进行随机森林二次降维。RFR是集成学习的一种方法，它将许多决策树集成在一起以提高泛化能力。在训练阶段，RFR使用bootstrap抽样从输入的训练数据集中收集几个不同的子训练数据集，依次训练几个不同的决策树。目的是使用几个不同的子模型，以增加最终模型预测结果的稳定性。对于切片变量和切片点的选择，我们遍历每个特征和每个特征的所有值，然后通过测量分割后节点的杂质(即每个子节点G(x, v)的杂质的加权和)来确定最佳的分割变量和分割点。G(x, v)的公式由式(7)给出：

$G\left(x_i+v_{ij}\right)=\frac{n_L}{N_s}H\left(X_L\right)+\frac{n_R}{N_s}H\left(X_R\right)$ (7)

其中x_i是分割变量之一；v_ij是分割变量的一个值；n_L、n_R、N_s分别表示分割后左右子节点的训练样本数和当前节点的所有训练样本数；X_L和X_R分别为左右子节点的训练样本集；H(X)是杂质函数。分类和回归任务的杂质函数H(X)不同，本文采用均方误差(MSE)作为杂质函数。因此，G(x, v)的表达式由式(8)给出：

$G\left(x,v\right)=\frac{1}{N_s}\left({\displaystyle \underset{y_i\in X_L}{\sum }{\left(y_i-{\bar{y}}_L\right)}^2}+{\displaystyle \underset{y_j\in X_R}{\sum }{\left(y_j-{\bar{y}}_R\right)}^2}\right)$ (8)

而决策树中节点的训练过程等价于寻找分割变量和G(x, y)最小的分割点。

在RFR的基础上进行特征提取，计算特征的重要性。节点k的特征重要性由式(9)给出：

$n_k=w_k\times G_k-w_L\times G_L-w_R\times G_R$ (9)

其中w_k、w_L、w_R分别为节点k及其左右子节点的训练样本数与总训练样本数的比值；G_k、G_L、G_R分别为节点k及其左右子节点的杂质。在得到每个节点的重要性后，我们可以通过式(10)计算特征的重要性：

$f_i=\frac{{\displaystyle \underset{j}{\sum }{n}_j}}{{\displaystyle \underset{m}{\sum }{n}_m}}$ (10)

其中f_i为特征i的重要性，j为特征i上划分的节点，m为所有节点。由式(11)给出归一化后的特征重要性公式：

$f_{norm,i}=\frac{f_i}{{\displaystyle {\sum }_{j\in allfeatures}{f}_j}}$ (11)

将f_norm由大到小排序后，可以得到 $f_{norm}^{\ast }$ 和累积重要性C_i。

与传统方法不同，该方法可以对不同输出的重要变量进行筛选。由于在RFR中充分利用了不足的样本，特征提取的结果比传统的灵敏度分析更加可信。那些累积重要性通常大于0.95的特征可以被认为是各种输出的重要特征，构建RFR的总体过程如图2所示。

Figure 2. The overall structure process of RFR

图2. RFR的总体结构过程

这里，在Lasso回归降维基础上进行随机森林回归二次降维，根据变量对生物活性影响的重要性进行排序，并给出前20个对生物活性最具有显著影响的变量，结果如表1和图3所示。

3. 模型求解

本节利用上小节中筛选出来的不超过20个分子描述符变量，构建化合物对ERα生物活性的定量预测模型。然后使用构建的预测模型，对50个化合物进行IC50值和对应的pIC50值预测。本质上就是根据已有的变量数据，通过建立的模型去预测对应的pIC50值，同时用已知的样本数据验证所建立的模型的合理性。因为IC50值和对应的pIC50值之间存在函数关系，且pIC50值与生物活性值具有正相关性，即pIC50值越大表明生物活性越高；实际QSAR建模中，一般采用pIC50来表示生物活性值，且pIC50值幅度较小、精度更高，因此预测模型选择针对pIC50，然后利用函数关系求解IC50值。

Table 1. The top 20 variables with the most significant effects on biological activity

表1. 前20个对生物活性最具有显著影响的变量

Figure 3. Importance ranking diagram of variables

图3. 变量重要性排序图

3.1. 预测方法

在数学建模中，预测数据的模型有着很多不同的预测方法，并且不同的方法各有优劣，此处选择了三种不同的预测模型对数据进行预测，分别为径向基函数(RBF)、响应面法(RSM)和随机森林回归(RFR)，并通过三种方法的MAPE值进行对比寻找出最优的预测模型。

径向基函数(RBF)是一个取值仅仅依赖于离原点距离的实值函数。构造神经网络的基本方法为假设某种过程是属于某种函数空间的函数，然后连接成神经网格，运行一段时间该网络的电势趋于最小达到某种动态的平衡，从而可以求出该函数，而选择径向基函数空间是一个比较简单的容易用神经网络实现的方法。

响应面法(RSM)通过较少的试验在局部范围内比较精确的逼近函数关系，并用简单的代数表达式展现出来，计算简单。同时，通过回归模型的选择，可以拟合复杂的响应关系，具有良好的鲁棒性，并且还能获得显式表达。

随机森林回归(RFR)是集成学习的一种方法，它将许多决策树集成在一起以提高泛化能力。在训练阶段，RFR使用bootstrap抽样从输入的训练数据集中收集几个不同的子训练数据集，依次训练几棵不同的决策树。目的是使用几个不同的子模型，以增加最终模型预测结果的稳健性和稳定性。随机林的优点是：它具有非常高的预测准确率，对异常值和干扰具良好的容忍度，且不易出现过度拟合。

RBF神经网络 [7] 由输入层接收训练样本 $x_i\in R^d$ ，隐含层将输入样本通过径向基函数映射到新的空间，设隐含层节点数为M，若径向基函数为高斯函数，则 $c_i\in R^d$ 表示高斯函数的中心向量，δ_i表示高斯函数的核宽，由式(1)实现输入空间到新空间的映射：

$\phi \left(\Vert x-c_i\Vert \right)=\exp \left(-\frac{{\Vert x-c_i\Vert }^2}{{\delta }_i}\right)$ (12)

输出层在新空间实现线性加权求和，设w_i为隐含层与输出层的连接权值， $\phi$ 为径向基函数， $y\in R$ 为输出结果， $R^d\to R$ 的映射函数为：

$y=f\left(x\right)={\displaystyle \underset{i=1}{\overset{M}{\sum }}{W}_i}\varphi \left(\Vert x-c_i\Vert \right),i=1,2,\cdots ,M$ (13)

由此，RBF神经网络完成 $f:R^d\to R$ 的非线性映射。

响应面法先选定采取的近似数学多项式模型的类型，然后将样本数据输入进该数学多项式模型中进行求解，得出一个初步的近似化模型，再用原始的样本数据带入计算得出的数学模型中求解新的输出，与原输出进行比较，计算其残差平方和，并以残差平方和最小为目标，进行数学模型中最佳项的选择条件。

这里采用的是二阶的响应面多项式方程，因为一阶的响应面多项式方程无法拟合复杂的模型，因此不选取一阶响应面模型。而对于高阶的响应面方程需要较多的数据变量个数来满足，此时也不满足，因此我们选用二阶响应面方程模型来进行模型预测。二阶响应面通式如下表示：

$\hat{y}={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+\cdot \cdot \cdot +{\beta }_M{x}_M+{\beta }_{M+1}{x}_1^2+{\beta }_{M+2}{x}_2^2+\cdot \cdot \cdot +{\beta }_{2M}{x}_{2M}^2+{\displaystyle \underset{i\ne j}{\sum }{\beta }_{ij}{x}_i{x}_j}$ (14)

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

其中 $\hat{y}$ 为响应近似值，y_i为响应实际值，n是构造响应面模型的样本数，x_i为样本参数值。

在数学模型确定后，还需要确定反馈方法，这里采用的是顺序替换法，即从常数项开始拟合，每次增加一个项使回归平方和(SSR)最小，每增加一个项后，检查是否可以去掉或替换已经存在的项，同时使SSR更小。顺序替换法的计算量较小，方便快捷。

随机森林回归原理在上文已经介绍，这里就不在赘述。不同的是，前文随机森林回归计算的是Lasso回归降维之后的157个变量，而这里的随机森林回归基于上节给出的前20个对生物活性值影响最大的变量。对五十个化合物进行生物活性值预测，图4为随机森林树、节点图的部分展示。

Figure 4. Random Forest-Tree and node diagram (part)

图4. 随机森林–树、节点图(部分)

3.2. 精度验证

径向基函数和响应面法，误差与准确率比较接近。相比之下，随机森林回归具有更高的准确率，如表2和图5所示，将所有样本中75%的样本作为训练集，25%的样本作为测试集，误差验证方式为MAPE (8.34%)，准确率为91.66%。

$\text{MAPE}=\frac{100\text{\%}}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left|\frac{{\hat{y}}_i-y_i}{y_i}\right|}$ (16)

Table 2. Comparison of accuracy rates of the three methods

表2. 三种方法的准确率比较

Figure 5. Model accuracy test. (a) Comparison of the accuracy of the three methods; (b) Comparison of the error of the RFR model

图5. 模型精度检验 (a) 三种方法的准确率比较；(b) RFR模型误差对比图

平均绝对百分比误差(MAPE)范围[0, +∞)，MAPE为0%表示完美模型，MAPE大于100%则表示劣质模型。

由图5可知，RFR模型在对化合物的pIC50预测值与真实值有较高的拟合度，其8.34%的平均百分比误差与91.66%的准确率明显优于RBF模型和RSM模型。因此，使用基于RFR构建的预测模型，对50个化合物进行pIC50值预测，根据预测的pIC50值计算相对应的IC50值，其中IC50值是pIC50值的负对数，即IC50 = 10^(9-pIC50)，结果如表3和图6所示。

Table 3. Forecast results

表3. 预测结果

Figure 6. Prediction results of pIC50

图6. 目标化合物pIC50预测结果

4. 结论

本文对1974个化合物的729个分子描述符进行变量选择，利用Lasso回归与随机森林回归对数据进行降维，再选择径向基函数、响应面法和随机森林回归这三种不同的预测模型对数据进行预测，并通过三种方法的MAPE值进行对比寻找出最优的预测模型，其中随机森林回归模型最优，并且其精度最好。在建模过程中，已经证明了筛选出的变量有足够的能力去构建预测模型，并且通过调参提高了模型的精度。从优化的结果不难看出药物研发所用化合物的生物活性存在不小差距，而根据预测结果可以得出抗乳腺癌候选药物的目标化合物pIC50的生物活性，以此来筛选出有效活性化合物进行药物研发，对加快研发进展和降低研发成本具有一定的帮助。