1. 引言
乳腺癌是目前世界上最常见，致死率较高的癌症之一。乳腺癌的发展与雌激素受体密切相关，有研究发现，雌激素受体α亚型(Estrogen Receptors alpha, ERα)在不超过10%的正常乳腺上皮细胞中表达，但大约在50%~80%的乳腺肿瘤细胞中表达；而对ERα基因缺失小鼠的实验结果表明，ERα确实在乳腺发育过程中扮演了十分重要的角色。目前，抗激素治疗常用于ERα表达的乳腺癌患者，其通过调节雌激素受体活性来控制体内雌激素水平。因此，ERα被认为是治疗乳腺癌的重要靶标，能够拮抗ERα活性的化合物可能是治疗乳腺癌的候选药物。比如，临床治疗乳腺癌的经典药物他莫昔芬和雷诺昔芬就是ERα拮抗剂。目前，在药物研发中，为了节约时间和成本，通常采用建立化合物活性预测模型的方法来筛选潜在活性化合物。具体做法是：针对与疾病相关的某个靶标(此处为ERα)，收集一系列作用于该靶标的化合物及其生物活性数据，然后以一系列分子结构描述符作为自变量，化合物的生物活性值作为因变量，构建化合物的定量结构–活性关系(Quantitative Structure-Activity Relationship, QSAR)模型，然后使用该模型预测具有更好生物活性的新化合物分子，或者指导已有活性化合物的结构优化。
2. 数据来源
本文使用的数据集来自“中国研究生创新实践系列大赛”平台举办的“华为杯”第十八届中国研究生数学建模竞赛D题，数据包括1974个化合物的生物活性值所对应的729个分子结构描述符和化合物对ERα生物活性影响的情况，如表1、表2所示。
Table 1. 729 molecular structure descriptors for 1974 compounds
表1. 1974个化合物对应的729个分子结构描述符
Table 2. Effect of compounds on the biological activity of ERα
表2. 化合物对ERα生物活性的影响情况
3. 数据预处理
对原始数据和样本数据进行分析得：由于样本数量较大，取得途径并未详细告知，原始数据可能在采集过程中或传输过程中因种种原因出现数据缺失、数据异常等问题，导致数据质量下降。在这种情况下，如果直接使用这些不良数据进行后续得抗乳腺癌候选药物得优选建模，得到得结果不具有说服性。因此，考虑在进行用模型对分子描述符筛选之前，对样本原始数据进行预处理。
对729个分字描述符的预处理如下：首先，删除数值没有变化的变量进行数据清洗；其次，通过3σ准则将原始数据中的异常值进行处理以提高数据的质量；然后，计算每列分子描述符的变异系数，对变异系数小于设定值的分子描述符进行删除；最后利用相关系数删除和因变量不相关的变量。按照要求对经过上述方法筛选出来的变量用：多元线性回归模型、LASSO回归模型、随机森林模型进行分子描述符的筛选。最后结合熵权法和Topsis方法对所有的指标进行评价，对指标进行排序，选取前50个对生物活性最具有显著影响的分子描述符 [1] 。思路流程图如图1所示。
Figure 1. Flow chart of data preprocessing ideas
图1. 数据预处理思路流程图
3.1. 重复数据的处理
对于1974个化合物的各个分子描述符的数据进行分析，发现数据存在数值没有整列变化的情况。因此利用Python的“drop”函数进行数据的处理，删除了225个重复的分子描述符，之后利用Python筛选出重复的行有94行，并将其删除 [2] 。
3.2. 异常值的处理
利用Python中的“outrange”函数，用循环的方式，实现所有变量的异常值筛选，然后用该分子描述符的平均值替换异常值。通过重复数据处理后的504个分子描述符，经过3准则筛选后，筛查出了122个分子描述符后，然后用平均值替换异常值。
计算进行了异常值处理的504个分子描述符的变异系数，取其绝对值，通过分析，设变异系数的阈值为0.3，变异系数小于0.3的删除。经过3 准则筛选后，删除了71个分子描述符后，剩余433个分子描述符 [3] 。
3.3. 利用相关系数处理数据
构造包含x和y的数据框，利用Python计算相关系数，通过分析，设相关系数的阈值为0.025，则删除相关系数小于0.025的变量。通过变异系数处理后的433个分子描述符，经过相关系数法筛选后，删除了324个分子描述符后，剩余109个分子描述符。我们将用109个分子描述符进行之后的数学建模。
3.4. 利用模型筛选变量
3.4.1. 线性回归模型
1) 数据标准化
在数据分析之前，我们通常需要先将数据标准化(normalization)，利用标准化后的数据进行数据分析。数据标准化也就是统计数据的指数化。数据标准化处理主要包括数据同趋化处理和无量纲化处理两个方面。利用Python进行上述标准化处理，原始数据均转换为无量纲化指标测评值，即各指标值都处于同一个数量级别上就可以进行综合测评分析。
因此在利用模型筛选分子描述符之前我们需要先进行z-score标准化，及利用分子描述符原始数据的均值和标准差进行数据的标准差，其公式为：
2) 线性回归模型
线性回归(Linear Regression)是利用数理统计中的回归分析，来确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法，运用十分广泛。利用通过上述方法筛选出来的109个分子描述符建立线性回归模型
对变量系数
的重要性利用python计算，得出的结果如表3所示。
Table 3. Table of coefficient importance after linear regression
表3. 线性回归后的系数重要性表
3.4.2. LASSO回归模型
LASSO回归模型(Least Absolute Shrinkage and Selection Operator)是一种压缩估计，它通过构造一个惩罚函数得到一个较为精炼的模型，使得它压缩一些回归系数，即强制系数绝对值之和小于某个固定值；同时设定一些回归系数为零。因此保留了子集压缩的优点，是一种处理具有复共线性数据的有偏估计。LASSO回归的特点是在拟合广义线性模型的同时进行变量筛选和复杂度调整，因此，不论目标因变量是连续的，还是二元或者多元离散的，都可以用LASSO回归建模进行预测。这里的变量筛选是指不把所有的变量都放入模型中进行拟合，而是有选择的把变量放入模型从而得到更好的性能参数。
将上述通过了数据预处理后的109个分子描述符进行构建LASSO回归模型，利用Python计算得到结果如表4所示。
Table 4. Table of values after LASSO regression
表4. LASSO回归之后的数值表
3.4.3. 随机森林模型
随机森林(Random Forest)是一种比较新的机器学习模型集成学习方法，通过反复二分数据进行分类或回归，使得计算量大大降低。随机森林通过自助法重采样技术，从原始训练样本集N中有放回地重复随机抽取K个样本生成新的训练样本集，然后根据自助样本集生成n个分类树组成随机森林，其实质是对决策树算法的一种改进，将多个决策树合并在一起，每棵树的建立依赖于一个独立抽取的样本集。并且，随机森林算法自带特征筛选机制，即随机森林能够评估各个特征在相应问题上的重要性 [4] 。
将上述筛选出的109个分子描述符带入构建的随机森林算法，利用Python计算其的结果如表5所示。
Table 5. Random forest results table
表5. 随机森林结果表
3.4.4. 利用Topsis法筛选变量
Topsis (Technique for order preference by similarity to an ideal solution)是根据有限个评价对象与理性化目标的接近程度进行排序的方法，是在现有的对象中进行相对优劣的评价。Topsis法是多目标决策分析中一种常用的有效方法，又称为优劣解距离法。其基本原理是通过检测评价对象与最优解和最劣解的距离来进行排序，若评价对象最靠近最优解又最远离最劣解，则为最好，否则不为最优。其中最优解的各项指标值都达到各评价指标的最优值，最劣解的各指标值都达到各评价指标的最差值。将经过线性回归、LASSO回归，随机森林计算出的数据整合到一个表上，利用Python实现Topsis方法运行并进行分子描述符重要性排序，排序结果如表6所示。
3.5. 结果结论
通过对数据预处理和模型筛选，筛选出对生物活性影响有重要影响的分子描述符，并通过Topsis法根据分子描述符对生物活性影响的重要性进行了排序。提取前50个对生物活性最具有显著影响的分子描述符，如表7所示。
Table 6. Topsis method score table
表6. Topsis法得分表
Table 7. The top 50 molecular descriptors with the most significant impact on biological activity
表7. 前50个对生物活性最具显著影响的分子描述符
本文首先删除数值没有变化的变量进行数据清洗；其次，通过3准则将原始数据中的异常值进行处理以提高数据的质量；然后，计算每列分子描述符的变异系数，变异系数特别小的说明整列数据无太大变化，则影响力较小，故对变异系数小于设定值的分子描述符进行删除；最后利用相关系数删除和因变量不相关的变量。按照要求对经过上述方法筛选出来的变量用：多元线性回归模型、LASSO回归模型、随机森林模型进行分子描述符的筛选。最后结合熵权法和Topsis方法对所有的指标进行评价，对指标进行排序。极大地保证了分子描述符筛选过程和结果的合理性和完善性。
4. 大脑情感网络模型和PSO-AGA算法
与其他模型相比，BEL模型具有收敛速度快、运行时间短、精度高等优点。PSO-AGA算法能够在全局搜索和局部搜索能力之间实现平衡，提高了检测精度和稳定性。与其他优化算法相比，PSO-AGA算法具有更好的性能和预测精度，BEL模型和PSO-AGA具有广泛的应用价值 [4] 。
4.1. 改进的大脑情感模型
大脑情感学习(Brain Emotional Learning, BEL)模型在模拟人类大脑中的杏仁体(amygdala)组织和眶额皮质(orbitofrontal cortex)组织之间信息传递方式的基础上，建立的是一种新型机器学习模型。情绪神经科学研究表明，情绪的脑机制主要是前额叶皮层和边缘系统。边缘系统的核心部分由杏仁核和海马体组成，用于情感学习和记忆。感觉刺激到达丘脑后，可以通过下通道直接发送到杏仁核，也可以先通过上通道发送到感觉皮层，在那里仔细处理感觉刺激，然后将信号发送到杏仁核 [5] 。
本文使用Chen Dan等人提出的一种改进的BEL模型，该模型包含三个隐藏层，即感觉皮层、眼窝前额皮层和杏仁核。该模型能够提高数据的预测精度 [6] 。
4.2. PSO-AGA
该算法在搜索能力上有较大的提升空间，但收敛速度较慢。PSO算法具有一定的记忆功能，但容易陷入局部最优 [4] 。将遗传算法与粒子群算法相结合，不仅提高了收敛速度，而且提高了全局搜索能力。可以在更大的空间内搜索粒子，而不局限于之前的最优位置，这将提高算法的灵活性，增强优化能力。
在PSO-AGA算法中，首先初始化遗传参数，如最大最小交叉概率、最大最小突变概率、迭代次数等；其次，对种群和染色体进行初始化，然后使用PSO算法更新种群，使用AGA算法进一步更新种群。然后将最优染色体保存到下一代，直到达到最大迭代次数。最后，将最优染色体，即最优权重和阈值代入BEL模型，输出最终预测结果 [7] 。
4.3. 基于大脑情感网络模型和PSO-AGA算法的预测效果检验
数据首先进入丘脑，其中的最大值会经过下通道直接传入杏仁核，然后其他数据的会先在感觉皮层处经过一个初步的处理，处理过后的数据会传入杏仁核和眶额叶皮层，模型的输出值即为预测值。杏仁核的输出值减去眶额叶皮层的输出值，其中使用PSO-AGA算法去优化其中的参数和权值，使最终输出结果更好，预测效果更优 [8] 。得到基于大脑情感网络模型和PSO-AGA算法的预测效果检验的结果为，拟合系数为0.80931，均方误差根为0.16188，平均绝对误差为0.32860。
4.4. 与其他模型的预测效果进行对比
使用多种方法构建化合物对ERα生物活性的定量预测模型，找出ERα生物活性值与选出的50个分子描述符之间的逻辑关系。首先，考虑到筛选出的分子描述符之间存在潜在的化学联系机理，同时操作变量高度非线性、相互关联，此时传统的基于线性关系的回归预测模型可能不再适用，基于上述考虑，文本计划使用基于大脑情感网络模型和PSO-AGA算法构建预测模型，并与传统的多元线性回归模型、支持向量机、随机森林模型和梯度提升树回归模型四种模型进行比较，得到基于大脑情感网络模型和PSO-AGA算法构建的预测模型更有先进性和有效性。使用基于大脑情感网络模型和PSO-AGA算法构建的预测模型对文件“ERα _activity.xlsx”的test表中的50个化合物进行IC50值和对应的pIC50值预测。
4.4.1. 基于多元线性回归对ERα生物活性的预测模型求解并检验
利用Python求解基于多元线性回归对ERα生物活性的预测模型，首先将筛选出的50个分子描述符(maxdssC、nF8Ring、…、MDEO-22、minHBint5)设为自变量(
)建立多元线性回归方程，计算各个分子描述符对应的回归系数，利用Python将“ERα _activity.xlsx”文件中中的IC50_nM列及对应的pIC50对多元线性回归模型进行训练和检验。得到线性回归模型得拟合系数为0.54750，均方误差根为0.93934，平均绝对误差为0.75577。
4.4.2. 基于随机森林对ERα生物活性的预测模型求解并检验
利用Python构建基于随机森林对ERα生物活性的预测模型，并用Python将“ERα_activity.xlsx”文件中中的IC50_nM列及对应的pIC50对随机森林模型进行训练和检验，得知随机森林模型得拟合系数为0.71770，均方误差根为0.74193，平均绝对误差为0.54589。
4.4.3. 基于支持向量机对ERα生物活性的预测模型求解并检验
利用Python构建基于支持向量对ERα生物活性的预测模型，并用Python将“ERα_activity.xlsx”文件中中的IC50_nM列及对应的pIC50对支持向量机模型进行训练和检验，可知支持向量机模型得拟合系数为0.54978，均方误差根为0.93697，平均绝对误差为0.73247。
4.4.4. 基于梯度提升树回归模型对ERα生物活性的预测模型求解并检验
利用Python构建基于梯度提升树对ERα生物活性的预测模型，并用Python将“ERα_activity.xlsx”文件中中的IC50_nM列及对应的pIC50对梯度提升树回归模型进行训练和检验，可知梯度提升树回归模型的拟合系数为0.66972，均方误差根为0.80252，平均绝对误差为0.56077。
4.5. 模型验证与结果分析
将基于大脑情感网络模型和PSO-AGA算法构建的预测模型与4.4中各个预测模型的检验整合成表，其中R2为模型的拟合系数，RMSE为均方误差根，MAE为平均绝对误差，模型预测效果对比如表8所示。
Table 8. Summary table of model prediction results tests
表8. 模型预测结果检验汇总表
通过表中的数据发现基于大脑情感网络模型和PSO-AGA算法构建的预测模型的拟合程度最好，拟合系数远高出其他模型，均方根误差和平均绝对误差也是各模型最小的。综上所述，体现了基于大脑情感网络模型和PSO-AGA算法构建的预测模型在当前的很多数据集上，相对其他算法有着很大的优势。
通过Python将构建的预测模型进行求解，将对文件“ERα_activity.xlsx”的test表中的50个化合物进行IC50值和对应的pIC50值预测结果填入表格，如表9所示。
5. 总结
通过探究能够拮抗ERα活性的化合物，找到可能影响乳腺癌治疗药物的化合物活性的分子符对于研制乳腺癌的治疗候选药物，减少乳腺癌的致死率具有较大的实际临床意义。基于大脑情感模型的PSO-AGA算法构建影响拮抗ERα活性模型能够科学的建ERα生物活性值与分子描述符之间的逻辑关系，对优化ERα拮抗剂的生活物性具有较强的指导意义 [9] 。
本文又建立了基于多元线性回归的对ERα生物活性预测模型、基于支持向量机对ERα生物活性预测模型、基于随机森林对ERα物活性预测模型、基于梯度提升树回归模型对ERα生物活性预测模型。并且对不同的模型进行综合评价，极大地保证了基于大脑情感模型的PSO-AGA 算法构建影响拮抗ERα活性模型预测结果的合理性和完善性。
基金项目
重庆理工大学研究生教育高质量发展行动计划资助成果。
参考文献
NOTES
*通讯作者。