基于随机森林回归的乳腺癌抗药生物活性预测

Prediction of Antibiotic Activity of Breast Cancer Drug Resistance Based on Random Forest Regression

摘要: 在药物研发中，雌激素受体α亚型(Estrogen receptors alpha, ERα)被认为是治疗乳腺癌的重要靶标，能拮抗ERα活性的化合物可能是治疗乳腺癌的候选药物。因此本文旨在以生物活性值pIC50作为因变量，作用于ERα靶标化合物的分子结构描述符作为自变量，构建关于ERα靶标化合物的生物活性预测模型，进而挑选出有效的抗癌候选药物。首先采用方差过滤法、随机森林、XGBoost以及灰色关联分析对自变量进行筛选，得到MDEC-23等16个与pIC50相关性强，且变量间相关性弱的分子结构描述符。其次建立随机森林回归生物活性预测模型，将预测结果与支持向量回归、梯度提升回归树、XGBoost模型和MLP回归模型预测结果进行对比分析，结果表明随机森林回归模型能更好地拟合数据，在R2、MAE、MSE上优于其它模型，更适应于对生物活性pIC50值的预测，同时也表明筛选出的分子结构描述符在一定程度上能治疗乳腺癌。

Abstract: In drug development, the Estrogen receptors α subtype (ERα) is considered to be an important target for the treatment of breast cancer, and compounds that antagonize ERα activity may be candidates for the treatment of breast cancer. Therefore, this paper aims to use the biological activity value pIC50 as the dependent variable and the molecular structure descriptor acting on the ERα target compound as the independent variable to construct a prediction model for the biological activity of the ERα target compound, and then select effective anti-cancer drug candidates. Firstly, the independent variables were screened by variance filtering, random forest, XGBoost and gray correlation analysis, and 16 molecular structure descriptors such as MDEC-23 with strong correlation with pIC50 and weak correlation between variables were obtained. Secondly, a random forest regression biological activity prediction model is established, and the prediction results are compared and analyzed with the prediction results of support vector regression, gradient boosting regression tree, XGBoost model and MLP regression model, and the results show that the random forest regression model can better fit the data, is better than other models in R2, MAE and MSE, and is more suitable for predicting the bioactive pIC50 value, and also shows that the screened molecular structure descriptors can treat breast cancer to a certain extent.

1. 引言

乳腺癌是最常见的恶性肿瘤类型，也是全球女性癌症死亡的第二大原因 [1] 。2020年全球新发乳腺癌226.14万例，占总体癌症发病的11.7%，死亡68.50万例，占总体癌症死亡的6.9% [2] ，图1是2020年全球位于前十位癌症的数据统计图，由该图可知乳腺癌已成为全球范围内发病率最高的癌症，且死亡率已居所有恶性肿瘤的第五位 [3] 。据“我国进展期乳腺癌共识指南2020 (CABC3)”报道，国内乳腺癌发病率逐年升高，每年10万人中约有545.29人患乳腺癌 [3] [4] 。因此研发治疗乳腺癌的药物是十分必要的。

Figure 1. Statistics chart of the top 10 cancer in the world in 2020

图1. 2020年全球前十位癌症数据统计图

在药物研发中，雌激素受体α亚型(Estrogen receptors alpha, ERα)被认为是治疗乳腺癌的重要靶标 [5] [6] [7] ，因此能够拮抗ERα活性的化合物可能是治疗乳腺癌的候选药物。其次，为了节约时间与成本通常采用建立“化合物活性预测模型”的方法来筛选潜在活性化合物 [8] 。具体做法是：针对与疾病相关的某个靶标(此处为ERα)，收集一系列作用于该靶标的化合物及其生物活性数据，然后以一系列分子结构描述符作为自变量，化合物的生物活性值作为因变量，构建化合物的“定量结构–活性关系模型”，然后使用该模型预测具有更好生物活性的新化合物分子，或者指导已有活性化合物的结构优化。传统机器学习方法尤其是随机森林(rand forest, RF) [9] 、支持向量机和人工神经网络在药物活性方面能够达到较高的预测精度 [10] 。

综上所述，本文以作用于ERα靶标化合物的分子结构描述符作为自变量，生物活性pIC50作为因变量，利用数据挖掘技术、机器学习等方法，先对变量进行筛选，变量合理验证等过程后构建关于ERα靶标化合物的生物活性预测模型，挑选出有效的抗癌候选药物，进而对生物活性pIC50进行预测，为乳腺癌的药物研发上提供理论支撑。

2. 数据描述与自变量筛选

2.1. 数据描述

共有1975 ´ 730个数据，如表1所示，其中第一列为化合物的结构，最后一列为生物活性pIC50值，其余列均为化合物的分子结构描述值。

Table 1. The first 10 rows of data display table

表1. 前10行数据展示表

Figure 2. The standard deviation of the independent variable

图2. 自变量的标准差图

数据描述性分析

计算出729个自变量的方差、标准差、均值、四分位点值，部分值如表2所示，由这些统计量可知，自变量之间的取值范围较大，建模前应该对数据进行标准化处理，再结合图2，存在大量数据波动范围小的自变量，其自变量对建模效果产生的影响较小，可利用方差过滤法去除。

Table 2. Table of statistical values for the top 10 independent variables

表2. 前10个自变量的统计量值表

2.2. 自变量筛选

在做预测模型时常常需要对变量进行筛选以达到降低变量个数、提升模型预测效果以及更便于实际应用的目的。本文中自变量共有729个，对于变量的筛选是很有必要的。

Figure 3. Variable filtering flowchart

图3. 变量筛选流程图

2.2.1. 筛选流程

变量筛选流程如图3所示，通过简单观察数据，发现存在大量取值为0的列，以及存在较多数据波动不大的列(见图2)，所以首先选择低方差过滤法进行变量选择，共过滤掉225个单一值特征变量，取得504个候选变量；其次利用随机森林、极致梯度提升(eXtreme Gradient Boosting, XGBoost) [11] [12] 对504个候选变量分别进行筛选，再取交集，共得到35个候选变量，其中随机森林和XGBoost的变量累积重要度的取值为80%；最后使用灰色关联分析得到20个候选变量(35个候选变量的灰色关联度值见表3)，其中灰色关联算法的分辨系数值选择0.5；最终利用最大信息数(Maximal Information Coefficient, MIC)、距离相关系数、以及皮尔逊相关系数确定建模变量，共得到16个变量，其中MIC基于信息的非参数性探索，用于衡量两个变量X和Y之间线性或非线性的关联强度，常用于机器学习的特征选择性，皮尔逊相关系数删除与生物活性pIC50存在较大线性相关的变量，距离相关系数删除与pIC50存在较大非线性相关的变量。

根据图3了解到在变量筛选上：

1) 当累积重要度选择为0.8时，XGBoost筛选的效果比随机森林的好，首先在筛选变量个数上少于随机森林筛选的变量个数，其次在49个变量中就有35个变量与随机森林筛选结果相同，保留了随机森林20%的结果；

2) 利用机器学习筛选出的变量还存在较高的关系，如自变量之间、自变量与因变量之间存在线性相关或着非线性相关，必须对挑选的变量进行合理性验证。

Table 3. Grey correlation values for 35 candidate variables

表3. 35个候选变量的灰色关联度值

2.2.2. 变量合理性验证

20个变量与pIC50的MIC值、皮尔逊相关系数值、距离相关系数值如表4所示，由该表可知：

1) 此20个特征变量与生物活性pIC50之间的距离相关系数和最大信息数MIC都较大，体现了其对建模目标影响较大，验证了本题变量选择的合理性；

2) 所有变量的Pearson相关系数绝对值都未超过0.6，可知这些特征变量与生物活性pIC50的线性相关性比较低，则在预测模型的选择上，非线性的回归模型比传统线性回归模型会有更好的表现；

3) 根据距离相关系数判别，可知还存在与生物活性pIC50中相关的特征变量，如nHBAcc、C1SP2。

Table 4. Verify the information table properly

表4. 合理验证信息表

注：距离相关系数判别：中相关0.4~0.6，弱相关0.2~0.4，极弱或无相关0~0.2。

为了避免自变量之间存在较高的相关性而对建模结果产生的影响，计算出这20个变量之间的距离相关系数，绘制得到图4，由该图可知20个自变量之间存在距离相关系数大于0.8的现象，存在高度相关关系，如vpc-5和vpc-4，maxHBint10和SHBint10，maxHBint10和minHBint10等等，因此还需要对变量进一步筛选。

Figure 4. Graph of distance correlation coefficient between variables

图4. 变量间的距离相关系数图

2.2.3. 变量最终确认

在选择出的对生物活性影响最大的20个自变量(见表4)重，由于变量之间含有部分高度相关，若直接选用这20个变量参与建模，则并不能最大程度的得到最大量的有用信息，造成信息的浪费，因此还需要重新挑选出一个不超过20个变量的特征子集。根据变量之间距离相关系数图(见图4)，对于变量之间相关性大于0.8的变量，考虑将其中一个自变量删除，最终剩余16个变量如表5所示。

Table 5. 16 molecular descriptors

表5. 16个分子描述符

3. 基于随机森林回归生物活性预测模型

3.1. 随机森林回归的基本思想

首先随机森林回归模型作为非线性生物活性预测模型。随机森林 [13] [14] 的基本思想是利用bootstrap重抽样方法从原始样本中抽取多个样本，对每个bootstrap样本构建决策树，然后将所有决策树预测平均值作为最终预测结果。随机森林回归可以看成是由很多弱预测器(决策树)集成的强预测器。

本文实现的RF是将多个二叉决策树打包组合而成的，训练RF便是训练多个二叉决策树。在训练二叉决策树模型的时候需要考虑怎样选择切分变量、切分点以及怎样衡量一个切分变量、切分点的好坏。针对于切分变量和切分点的选择，采用穷举法，即遍历每个特征和每个特征的所有取值，最后从中找出最好的切分变量和切分点；针对于切分变量和切分点的好坏，一般以切分后节点的不纯度来衡量，即各个子节点不纯度的加权和，其计算公式如下：

$G\left(x_i,v_{ij}\right)=\frac{n_{left}}{N_s}H\left(X_{left}\right)+\frac{n_{right}}{N_s}H\left(X_{right}\right)$ (1)

其中，x_i为某一个切分变量，v_ij为切分变量的一个切分值，n_left、n_right、N_s分别为切分后左子节点的训练样本个数、右子节点的训练样本个数以及当前节点所有训练样本个数，X_left、X_right分别为左右子节点的训练样本集合，H(X)为衡量节点不纯度的函数，在本题中选用MSE作为模型的不纯度函数。

3.2. 随机森林回归模型的建立

3.2.1. 建模过程

1) 论最优特征子集的选取

根据最终变量的确定，选择了16个变量参与构建生物活性预测模型，具体变量见表5。

2) 数据标准化

数据标准化的公式如下所示：

$x^{\ast }=\frac{x-\bar{x}}{\sigma }$ (2)

进行数据标准化的原因有：① 将不同量级的数据统一转化为同一个量级，保证数据之间的可比性；② 将数据拉回成均值为0，标准差为1的数据有利于回归模型的收敛。

3) 数据划分

按8:2的比例将1974行数据划分成训练集和测试集，样本比为1580:394。用训练集训练模型，再用训练好的模型在测试集上验证效果。

4) 建立模型

5) 模型预测效果评估

本文采用 $R^2$ 、MSE、MAE评价模型预测效果，公式如下：

$R^2=1-\frac{{\displaystyle \underset{i}{\sum }{\left({\hat{y}}_l-y_i\right)}^2}}{{\displaystyle \underset{i}{\sum }{\left({\bar{y}}_l-y_i\right)}^2}}$ (3)

$MSE=\frac{1}{m}{\displaystyle \underset{i=1}{\overset{m}{\sum }}{\left(y_i-{\hat{y}}_l\right)}^2}$ (4)

$MAE=\frac{1}{m}{\displaystyle \underset{i=1}{\overset{m}{\sum }}\left|y_i-{\hat{y}}_l\right|}$ (5)

其中， $y_i$ ， ${\hat{y}}_l$ 分别是测试集上的真实值和预测值。

3.2.2. 建立随机森林回归模型

对随机森林参数进行设置，设置树的个数为500棵，最大深度为10，得到随机森林回归模型，将其在测试集上的预测效果显示如图5所示，从中可以看出预测值与真实值走势大致相同，因此认为建立的模型有效。

4. 模型预测效果对比

为了说明所建立的随机森林模型有效性与优越性，按照同样的步骤对同一个训练集建立支持向量回归(Support Vector Regression, SVR)、梯度提升回归树、XGBoost回归以及多层感知机(MLP)回归四个模型。将五个模型的真实值与预测值对比图绘制如图6所示，由该图也可以看到随机森林回归很大程度上能更好地拟合真实值，更适应于ERα生物活性的预测。

Figure 5. True value vs predicted value

图5. 真实值vs预测值

Figure 6. Model comparison diagram

图6. 模型对比图

Table 6. Indicator evaluation table

表6. 指标评价表

为了避免偶然性，每个模型训练50次，计算五个模型的 $R^2$ 值、MSE值、MAE值(见表6)，发现随机森林模型的 $R^2$ 值最高，相对于支持向量回归、梯度提升回归树、XGBoost回归以及MLP回归依次提高了8.557%、6.126%、1.265%、20.518%，同时与文献 [15] 建立的GWO-KELM模型和支持向量回归相比，本文的 $R^2$ 分别提高了1.722%、13.22%。从MSE值、MAE值上，随机森林回归模型的值也是最小的，与文献 [15] 建立的GWO-KELM模型相比分别降低了4.583%、3.586%。说明了本文所选的分子结构描述符更能准确地表达pIC50值，以及随机森林回归在生物活性pIC50上的预测效果更好。

5. 总结

基于机器学习方法在预测药物活性方面能够达到较高的预测精度的优点，本文利用随机森林等方法实现抗癌候选物的研发，与传统方法相比极大地节约了时间和成本，以及减少人工带来的误差。为了保证模型的可用性以及提高模型的预测效果，第一步选择对变量进行特征选取。在变量筛选上，本文层层递进，先挑选了对靶标ERα影响较大的20个化合物，再利用MIC、距离相关系数、皮尔逊相关系数对变量进行合理性选择，避免多重共线性问题，最终挑选出16个分子结构描述符，其中约1/2的变量与何毅 [9] 和秦雅琴等人 [8] (见表5标红色的分子结构描述符)所筛选的变量一致，保留了大部分能治疗乳腺癌的分子结构。第二步建立随机森林回归模型，同时为了说明该模型的有效性以及优越性，建立了支持向量回归 [16] 、梯度提升回归树、XGBoost回归以及MLP回归模型。最后结合评价指标值、预测效果图的结果以及参考文献可知随机森林回归模型的预测效果在这五个模型当中是最优的，表明随机森林更适应于生物活性值的预测，该结论与叶丹等人 [17] 所得的结论相同，所以在对ERα的生物活性pIC50值的预测时，本文建议可选择随机森林回归。