基于神经网络抗乳腺癌候选药物优化建模

Optimization Modeling of Anti Breast Cancer Drug Candidates Based on Neural Network

摘要: ERα在乳腺发育过程中具有十分重要的作用，被认为是治疗乳腺癌的重要靶标。在化合物候选过程中，其不仅需要具备良好的生物活性，还需要具备良好的药代动力学性质和安全性。本文采用数据挖掘技术，根据提供的ERα拮抗剂信息，构建化合物生物活性的定量预测模型和ADMET性质的分类预测模型，从而为同时优化ERα拮抗剂的生物活性和ADMET性质提供预测服务，具有一定的现实意义。

Abstract: ERα plays a very important role in breast development and is considered as an important target for the treatment of breast cancer. In the candidate process of compounds, they need not only good biological activity, but also good pharmacokinetic properties and safety. In this paper, data mining technology is adopted according to the provided ERα Antagonist information, to build a quantitative prediction model of compound biological activity and a classification prediction model of ADMET properties, so as to simultaneously optimize ERα. The biological activity and ADMET properties of antagonists provide prediction services, which has certain practical significance.

1. 引言

据肿瘤数据库GLOBOCAN近期发布的全球癌症统计报告数据显示，乳腺癌已超肺癌，成为全球癌症新发病例数首位，同时也是中国人数癌症发病的第四大原因。乳腺癌的发展与雌激素受体密切相关，经典的雌激素受体，包括ERα、ER 这两种亚型，并且两者的结构相似，它们位于细胞核内，通过调节特异性靶基因的转录发挥“基因型”调节效应，雌激素受体如图1所示。

Figure 1. Estrogen receptor

图1. 雌激素受体

有研究发现，雌激素受体α亚型(Estrogen receptors alpha, ERα)在不超过10%的正常乳腺上皮细胞中表达，但大约在50%~80%的乳腺肿瘤细胞中表达；而对ERα基因缺失小鼠的实验结果表明，ERα [1] 确实在乳腺发育过程中扮演了十分重要的角色。目前，抗激素治疗常用于ERα表达的乳腺癌患者，其通过调节雌激素受体活性来控制体内雌激素水平。因此，ERα被认为是治疗乳腺癌的重要靶标，能够拮抗ERα活性的化合物可能是治疗乳腺癌的候选药物，具体做法如图2所示。

一个化合物想要成为候选药物，除了需要具备良好的生物活性(此处指抗乳腺癌活性)外，还需要在人体内具备良好的药代动力学性质和安全性，合称为ADMET (Absorption吸收、Distribution分布、Metabolism代谢、Excretion排泄、Toxicity毒性)性质。其中，ADME主要指化合物的药代动力学性质，描述了化合物在生物体内的浓度随时间变化的规律，T主要指化合物可能在人体内产生的毒副作用。一个化合物的活性再好，如果其ADMET性质不佳，比如很难被人体吸收，或者体内代谢速度太快，或者具有某种毒性，那么其仍然难以成为药物，因而还需要进行ADMET性质优化。一般需着重考虑化合物的5种ADMET性质，即小肠上皮细胞渗透性、细胞色素P450酶、化合物心脏安全性评价、人体口服生物利用度、微核试验。

Figure 2. Selection process of candidate drugs

图2. 候选药物的选取过程

2. 寻找影响生物活性物的主要因素

2.1. 数据预处理

为保证变量降维的可靠性，需预先对变量数据进行处理，检查是否存在空值的位点，如存在，则需使用当前变量数据均值替代缺失值；检查是否存在变量异常值，经查阅文献，若当前变量数据65%以上为恒定值，则可将该变量判定为无意义变量，故将该变量当作异常变量进行剔除。

1) 缺失值处理

本文给出1974个化合物对应729个分子描述变量的原始数据，将其数据表导入SPSS软件中，进行遍历搜索，若该位点为空，则判定为缺失值，最终搜索的缺失值数目为0。

2) 异常值处理

缺失值处理结束之后，为了数据的可靠性，再进行异常值处理，由于各分子描述符变量所对应不同的化合物，故各分子描述符变量下的数据并没有离群这一说法，因为并不是同一样本下测得的数据。经查阅相关文献，当某个分子描述符变量下存在65%以上的数据都相同，则表明，该分子描述符变量对各个化合物的生物活性没有规律性影响，故可将该分字符变量当作异常变量，并进行剔除操作。

3) 数据标准化

在变量降维之前，需要对各变量数据进行标准化操作，主要是因为各个变量的性质、量纲、可用性等特性可能会存在差异，导致降维无法得出具有代表性的变量。标准化能使各变量具有均等权重，同时，下文需要通过随机森林回归算法对预处理后的247个变量进行降维，故数据标准化具有重要作用。

目前，数据标准化的方法有很多，一般有极值法、标准差法以及三折线法等等。然而不同的标准化方法，可能对系统变量的评价结果也会有所不同。本文采用标准化方法——标准差法(Z-Score标准化)。通过标准化法 [2] 对这247个变量进行标准化处理，对应的计算公式为

$\overline{x_i}=\frac{x_i-\mu }{\delta }$ (1)

式中， $\overline{x_i}$ 为标准化后的第i个化合物对应的数据； $x_i$ 为第i个化合物； $\mu$ 为该变量的所有化合物下的数据均值； $\delta$ 为该变量的所有化合物下的数据标准差。

2.2. 变量降维

随机森林 [3] 进行回归分析时，可以在样本分类的基础上进行计算，由此获得各个重要性特征变量对目标值的影响程度即贡献度。贡献度越大，说明该特征变量对目标值的影响程度越大；反之则越小。关于某节点k的贡献度计算公式如下式所示。

$n_k={\omega }_k\ast G_k-{\omega }_{\text{left}}\ast G_{\text{left}}-{\omega }_{\text{right}}\ast G_{\text{right}}$ (2)

式中， ${\omega }_k$ ：节点k中训练样本与总训练样本数目的比例。

${\omega }_{\text{left}}$ ：节点k左节点中训练样本与总训练样本数目的比例。

${\omega }_{\text{right}}$ ：节点k右节点中训练样本与总训练样本数目的比例。

$G_k$ ：节点k处的不纯度。

$G_{\text{left}}$ ：节点k左节点的不纯度。

$G_{\text{right}}$ ：节点k右节点的不纯度。

由上述公式，可得出任一特征变量对目标值的贡献度公式为：

(3)

借助Matlab工具对随机森林回归算法进行应用，对预处理过后的247个变量进行贡献度计算，贡献度分布情况如图5所示。并将这247个变量针对贡献度进行由高到低排列，并设定贡献度阈值为0.2%，由此初步筛选得到贡献度排名前60个变量。

在不确定变量数据是否满足正态分布的前提下，不可直接使用Person相关性分析，本文采用Spearman相关系数 [4] 表示两变量间的关联程度，由此通过Matlab程序进行高耦合变量过滤。在统计学中，Spearman相关系数通常用来衡量两个变量间依赖性的非参数指标，用希腊字母 $\rho$ 表示。对与样本容量为n的样本，相关系数 $\rho$ 为：

$\rho =\frac{{{\displaystyle \sum }}_i\left(x_i-\bar{x}\right)\left(y_i-y\right)}{\sqrt{{{\displaystyle \sum }}_i{\left(x_i-\bar{x}\right)}^2{{\displaystyle \sum }}_i{\left(y_i-\bar{y}\right)}^2}}$ (4)

式中， $x_i,y_i$ 代表秩次， $\rho$ 表示相关系数， $\left|\rho \right|$ 越接近1，表明两变量相关程度越大。

通过Matlab得出上文60个变量间的相关系数，然后进行两两筛选，对具有高耦合性的变量组剔除低贡献度的变量，最后综合得出20个主要变量，以确保其独立性，主要变量间的相关度分布图如图3所示。

Figure 3. Distribution of correlation between main variables

图3. 主要变量间的相关度分布图

3. 建立化合物对ERα生物活性的定量预测模型

3.1. BPERα生物活性预测模型建立与验证

1) 数据整理

将上文中选取的20个主要变量对应的1974个化合物样本数据进行整理，然后将数据导入SPSS软件中进行标准化处理，处理以后的数据再导入Matlab中，将其中70%化合物样本数据作为训练集，将30%化合物样本数据作为测试集。

2) 网络结构建立

① 设置输入输出层：将问题1中选取的20个主要变量所对应的1974个化合物样本数据，经标准化之后导入Matlab中作为输入，输入层神经元数目设置为20个；将pIC₅₀值导入Matlab作为输出，输出层神经元数目设置为1个。

② 设置隐含层：隐含层的确定对整个网络训练过程起十分重要的作用。隐含层设置过少会导致网络训练性能达不到要求；而隐含层设置过多又会导致过拟合现象。经反复训练，最终确定隐含层神经元数目为20个。

根据BP (Back Propagation)网络 [5] 建立ERα生物活性预测模型后，得出20个主要变量对应的化合物样本测试集平均误差百分比在4%左右，个别样本所到最大误差百分比在12%左右，由此说明建立的预测模型比较良好。化合物样本总体误差比如图4，预测误差如图5，预测值与真实值如图6。

Figure 4. Overall error ratio of compound sample

图4. 化合物样本总体误差比

Figure 5. Compound sample error

图5. 化合物样本误差

Figure 6. Comparison between predicted value and real value

图6. 预测值与真实值的比较情况

3.2. ERα生物活性结果预测

根据上文Spearman相关系数对pIC₅₀值进行预测，得出pIC₅₀预测结果，可知pIC₅₀值是IC₅₀值的负对数，因此这两目标变量之间存在很强的相关性。再通过非线性拟合 [6] 得到pIC₅₀与IC₅₀之间的数学表达式，进而求得IC₅₀列中的值。

此外，均方误差(Mean Squared Error, MSE)也是评价预测模型好坏的重要指标，是预测值与真实值差值平方的平均值，本文模型的均方误差为0.254，误差较小，进一步说明本文建立的预测模型良好。

通过非线性拟合 [6] 得到pIC₅₀与IC₅₀之间的数学表达式，假设pIC₅₀与IC₅₀之间的逻辑关系式为：

$Y\left({\text{pIC}}_{50}\right)=a\ast \ln \left(X\left({\text{IC}}_{\text{50}}\right)\right)+b$ (5)

将文中pIC₅₀值与IC₅₀值导入Matlab中，进行拟合，拟合函数曲线如图7所示。由曲线拟合情况和基本信息可以看出，函数拟合优秀。

Figure 7. Fitting function curve

图7. 拟合函数曲线

4. 建立各ADMET性质的预测模型

Logistic回归为概率型非线性回归模型 [7] ，是研究二分类观察结果与一些影响因素之间关系的一种多变量分析方法。

考虑具有n个独立变量的向量 $x=\left(x_1,x_2,\cdots ,x_n\right)$，设条件概率 $P\left(y=1|x\right)=p$ 为根据观测量相对于某事件x发生的概率。Logistic回归模型可以表示为：

$P\left(y=1|x\right)=\pi \left(x\right)=\frac{1}{1+{\text{e}}^{-g\left(x\right)}}$ (6)

这里的 $f\left(x\right)=\frac{1}{1+{\text{e}}^{-x}}$ 称为Logistic函数。其中 $g\left(x\right)=w_0+w_1{x}_1+\cdots +w_n{x}_n$

在x条件下y不发生的概率为：

$P\left(y=0|x\right)=1-P\left(y=1|x\right)=1-\frac{1}{1+{\text{e}}^{-g\left(x\right)}}=\frac{{\text{e}}^{g\left(x\right)}}{1+{\text{e}}^{g\left(x\right)}}$ (7)

故事件发生的概率与不发生的概率之比如下式，也称作发生比。

$\frac{P\left(y=1|x\right)}{P\left(y=0|x\right)}=\frac{p}{1-p}={\text{e}}^{g\left(x\right)}$ (8)

对发生比取对数可得到：

$\ln \left(\frac{p}{1-p}\right)=g\left(x\right)=w_0+w_1{x}_1+\cdots +w_n{x}_n$ (9)

根据Logistic回归对得出的各组20个主要变量进行计算得出各组回归模型的 $w_0,w_1,\cdots ,w_{20}$，如下表1所示。并对1974个化合物的ADMET数据，分组进行训练，并用测试集数据进行预测。各组ADMET性质数据均分布在0~1之间，各组数据分布情况如下图8所示。

Table 1. wi Values of regression models in each group

表1. 各组回归模型的w_i值

(a) Caco-2 (b) CYP3A4 (c) hERG (d) HOB (e) MN

Figure 8. Distribution of ADMET property data in each group

图8. 各组ADMET性质数据分布情况

由图8可知，1974个化合物样本的Logistic回归值分布在0~1之间，以0.5为分界线，大于0.5的数据归1，小于0.5的数据归0。通过数据处理得出：Caco-2组测试数据的正确性为89.9%；CYP3A4测试数据的正确性为87.6%；hERG测试数据的正确性为88.2%；HOB测试数据的正确性为91.7%；MN测试数据的正确性为86.7%，说明本文回归模型较好。

5. ERα和ADMET性质同性影响的较优变量

首先对ADMET中的Caco-2、CYP3A4、hERG、HOB、MN进行正向化以保证各个指标值为1时性质为优，值为0时性质为差。进行多指标综合评价 [8] 后，对数据预处理后的247个操作符变量根据各个指标(共6个)求得相应的spearman系数，找到对ERα和ADMET中三个以上具有同相关性的变量，思路流程图如图所示。为确定各个变量的取值范围，本文通过题目给出的数据样本，选取90%样本分布区间作为变量的取值范围。然后将上述找到的变量与ADMET各性质间建立回归预测模型，进而得出ADMET性质的综合评价方程，再通过PSO粒子群算法寻ADMET性质优值及最优解，最后通过此时的变量最优值去验证ERα生物活性的优劣性。

1) 数据整理与ADMET评价模型建立

将8.2中选取的10个主要变量分别与Caco-2、CYP3A4、hERG、HOB、MN建立Logistic回归预测模型，为建立ADMET的评价模型，令：

$F=F_1+F_2+F_3+F_4+F_5$ (10)

F——ADMET的评价指标；

F₁——Caco-2的评价指标；

F₂——CYP3A4的评价指标；

F₃——hERG的评价指标；

F₄——HOB的评价指标；

F₅——MN的评价指标；

各指标的Logistic回归方程形式为

$F\left(x\right)=1-\frac{1}{1+{\text{e}}^{-g\left(x\right)}}=\frac{{\text{e}}^{g\left(x\right)}}{1+{\text{e}}^{g\left(x\right)}}$ (11)

其中：

$g\left(x\right)=w_0+w_1{x}_1+\cdots +w_n{x}_n$ (12)

由此得到各评价指标的回归系数，如下表2所示：

Table 2. Regression coefficient of each evaluation index

表2. 各评价指标的回归系数

各回归预测正确百分比为：F₁模型为86.75%；F₂模型为89.58%；F₃模型为85.32%；F₄模型为87.45%；F₅模型为91.53%，说明各性质回归模型较好。

2) PSO算法寻ADMET性质优值及最优解

PSO又叫粒子群算法，源于对鸟类捕食行为的研究，其核心思想是利用群体中的个体对信息的共享使整个群体的运动在问题求解空间中产生从无序到有序的演化过程，从而获得问题的可行解。粒子群算法流程图如下图9所示：

Figure 9. Flow chart of particle swarm optimization

图9. 粒子群算法流程图

PSO算法核心公式如下：

$v_i^d=w{v}_i^{d-1}+c_1{r}_1\left(pbes{t}_i^d-x_i^d\right)+c_2{r}_2\left(gbes{t}^d-x_i^d\right)$ (13)

其中：

$c_1$ ：粒子的个体学习因子；

$c_2$ ：粒子的社会学习因子；

w：速度的惯性权重；

$v_i^d$ ：第d次迭代时，第i个粒子的速度；

$x_i^d$ ：第d次迭代时，第i个粒子所在的位置；

$pbes{t}_i^d$ ：到第d次迭代为止，第i个粒子经过的最好的位置；

$gbes{t}^d$ ：到第d次迭代为止，所有粒子经过的最好位置；

$r_1,r_2$ ：[0, 1]上的随机数。

最终通过粒子群寻优，得到ADMET性质 [9] 的模型值为5，即说明Caco-2、CYP3A4、hERG、HOB、MN各性质都为优，此时得出10个变量的最优数值如下表3所示：

Table 3. Optimal values of 10 variables

表3. 10个变量的最优数值

3) 检验ERα生物活性

将上述变量所得得最优值代入问题2建立的BP神经网络中，计算得出此时的pIC₅₀值为6.8，该值说明ERα生物活性较好，故本题筛选出的这10个变量及求出的最优值符合题目要求。

6. 模型评价与改进

6.1. 模型评价

1) 模型充分考虑到729个分子描述符变量的异常情况，各分子描述符变量与生物活性之间存在高度非线性，以及各变量之间存在高度耦合性，故先对729个变量进行异常值筛选，得出247个对生物活性存在较大影响的变量，之后使用随机森林回归、Spearman相关系数进行二次降维，确保了变量与生物活性间的相关性和变量间的独立性。

2) 由于分子描述符变量与生物活性间存在较强的非线性关系，因此本文利用三层BP神经网络，再根据自学习的网络对pIC₅₀进行预测，最后得出预测误差比相对较小，预测模型精度表现良好。

3) 由于ADMET的5种性质的数值符合“0~1”分布，故采用logistic回归对其建立分类预测模型，得出5类模型的正确预测百分比均大于85%，说明模型精度表现良好。

6.2. 模型改进

由于BP神经网络计算的误差是通过不断地逐层传递来修正各层间地权值与阈值，因此本文基于PSO算法全局搜索BP神经网络的权值与阈值，并将得到的优化初值作为BP网络的初始权值与阈值，由此对文中的生物活性进行预测，最终得出训练样本的拟合优度为90.93%，测试集的拟合优度为87.11%，说明网络训练效果比较理想，拟合优度如下图10所示。

Figure 10. Goodness of fit

图10. 拟合优度

通过PSO优化后的BP网络，对ERα生物活性进行预测，得出化合物样本测试集平均误差比为3.74%，比优化前的平均误差比降低了19.05%，说明优化结果比较理想。化合物样本预测误差比如图11所示。

Figure 11. Optimized compound sample prediction error ratio

图11. 优化后的化合物样本预测误差比

7. 结语

1) 本文充分利用随机森林回归、Spearman相关性分析、BP神经网络、Logistic回归、PSO算法研究了如何建立乳腺癌候选药物的预测模型。

2) 根据本文所用的降维方法，得到ERα生物活性的20个主要变量，严格满足题目要求，既保证变量数据与目标数据的高相关性，且保证变量数据间的独立性，适用于高维度且大样本的变量筛选，对医疗药物成分数据建模具有较高参考意义。

3) 本文采用BP神经网络模型表达非线性系统，采用粒子群算法优化BP神经网络初始权值，改善预测误差，适用于医疗药物的影响预测。

4) 本文采用logistic回归算法用于0~1分类系统的预测，并得到预测准确率为85%以上的5类预测模型。

5) 通过同时考虑分子符描述变量对生物活性ERα与ADMET性质产生的正相关或负相关影响，筛选出满足条件的变量，并通过样本数据集，给出变量的取值范围，并给出获得最佳药物性能时变量的取值范围，若想得出更精确的取值范围，则需对这些变量一一进行实验探究。