Modeling and Simulation 建模与仿真, 2023, 12(2), 1583-1592 Published Online March 2023 in Hans. /journal/mos  /10.12677/mos.2023.122147
基于RFR 模型的抗乳腺癌候选药物优化
宛翔天，杨家麒，刘兴宇，李江涛
上海理工大学，机械工程学院，上海
收稿日期：2023年2月21日；录用日期：2023年3月22日；发布日期：2023年3月29日
摘
要 在抗乳腺癌药物研发中，为了节省时间成本，建立化合物预测模型来筛选活性化合物是一种有效的方法。本文根据提供的乳腺癌靶标雌激素受体α亚型(Estrogen receptors alpha, ER α)拮抗剂信息，利用Lasso 回归与随机森林相结合的方法，对数据降维并筛选出影响生物活性的主要变量，明确了模型建立的方向；在此基础上，建立了关于生物活性的随机森林回归模型并进行了预测。
关键词
抗乳腺癌药物研发，预测模型，Lasso 回归，随机森林回归
Optimization of Candidate Drugs for  Anti-Breast Cancer Based on RFR Model
Xiangtian Wan, Jiaqi Yang, Xingyu Liu, Jiangtao Li
School of Mechanical Engineering, University of Shanghai for Science and Technology, Shanghai  Received: Feb. 21st , 2023; accepted: Mar. 22nd , 2023; published: Mar. 29th , 2023
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. Accord-ing to the information of Estrogen receptors alpha (ER α) antagonist, which is the target of breast cancer trea
tment, 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 biolog-ical activity was established and predicted.
宛翔天等
Keywords
Anti-Breast Cancer Drug Development, Prediction Model, Lasso Regression, Random Forest
Regression Array Copyright © 2023 by author(s) and Hans Publishers Inc.
This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).
/licenses/by/4.0/
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为样本个数，1,2,,729
n= ；x m为Lasso初步降维筛选所得自变量，1,2,,157
r= 。
m= ；x r为随机森林回归二次降维筛选所得自变量，1,2,,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]。
普通的最小二乘法线性回归对应的线性模型如下：
()1122d d f x w x w x w x b =++⋅⋅⋅++                            (1)
用向量表示为，当用空间里的一个直观的模型取逼近这些数据点时，往往无法完全经过这些点，即模型理论值与数据值存在一定的偏差，其中b 为偏差值[5]。此时我们需要这些误差的总和最小，而由于实际值分布较为分散，往往在理论模型的上下两侧，所以若是将这些误差直接相加，则会出现正负抵消的现象，从而无法正确地估计误差的大小。因此，需要将每一项的误差项平方，即将误差项正向化，根据这
个方法，得出最小二乘法的误差项公式：
()()()()()222T T T T T 111N N N i i i i i i i i i L w f x y w x y w x y w X Y wX Y ====−=−=−=−−∑∑∑        (2)
将上式展开得，
()()()T T T T T T T T 2L w w X Y wX Y w X w X Y Y Y =−−=−+                  (3)
为了得出使误差项最小得w ，对其求偏导，得出
() 1T T w X X
X
Y −=                                (4)
宛翔天 等
但是在这种情况下，为了满足矩阵可逆的条件下，只对于目标值y 的个数大于变量值x 的个数，一旦变量值x 的个数过多，就会引起多重共线性和过拟合的现象[6]。
Lasso 回归为了解决这样得问题引入了扰动项构建正则化框架，使得矩阵变为可逆矩阵 ()()arg min w
L w p w λ+                                    (5)
其中，()P w w =，此时Lasso 得回归模型为
()2
21ˆlasso arg min w w y Xw w λ=−+                            (6) 其中11P
j j w w ==∑为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)给出：
()()()L R i ij L R s s
n n G x v H X H X N N +=+                          (7) 其中x i 是分割变量之一；v ij 是分割变量的一个值；n L 、n R 、N s 分别表示分割后左右子节点的训练样本数和当前节点的所有训练样本数；X L 和X R 分别为左右子节点的训练样本集；H (X )是杂质函数。分类和回归任务的杂质函数H (X )不同，本文采用均方误差(MSE)作为杂质函数。因此，G (x , v )的表达式由式(8)给出：
()()()221
,i L j R i L j R y X y X s G x v y y y y N ∈∈  =−+−
∑∑                      (8) 而决策树中节点的训练过程等价于寻分割变量和G (x , y )最小的分割点。
在RFR 的基础上进行特征提取，计算特征的重要性。节点k 的特征重要性由式(9)给出：
宛翔天 等
k k k L L R R n w G w G w G =×−×−×                            (9)
其中w k 、w L 、w R 分别为节点k 及其左右子节点的训练样本数与总训练样本数的比值；G k 、G L 、G R 分别为节点k 及其左右子节点的杂质。在得到每个节点的重要性后，我们可以通过式(10)计算特征的重要性：
j
j
i m
m n f n =∑∑                                    (10)
其中f i 为特征i 的重要性，j 为特征i 上划分的节点，m 为所有节点。由式(11)给出归一化后的特征重要性公式：
, i
norm i j j all features f f f ∈=∑                              (11)
将f norm 由大到小排序后，可以得到norm
f ∗和累积重要性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值。

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。