抗癌候选药物ERα抑制剂活性预测

夏玉兰，谢济铭，王雅婧，卢梦媛，王锦锐，秦雅琴

1）昆明理工大学交通工程学院，云南昆明 650504；2）温州医科大学第一临床医学院，浙江温州325006

乳腺癌是目前世界上最常见且致死率较高的癌症之一.根据世界卫生组织（Wold Health Organization，WHO）国际癌症研究机构（International Agency for Research on Cancer，IARC）2020年发布的全球癌症负担数据，乳腺癌新发病例及死亡率均列首位，成为全球“第一癌”［1］.然而，现有的乳腺癌临床治疗方法不能完全阻止乳腺癌的复发及转移扩散.因此，研发新型抗乳腺癌药物是治疗乳腺癌和改善乳腺癌患者预后的有效策略［2］.

在抗乳腺癌药物研发过程中，可通过对药物化合物分子结构特征、性质进行研究，给合数理统计的方法建立定量构效关系模型来实现药物活性预测.KIDERA 等［3］运用主成分分析（principal component analysis，PCA）从氨基酸的188 个特征参数中筛选出10 个最佳主成分变量，并采用多元线性回归对多肽药物进行定量表征与活性预测，取得了较好的结果.王青艳等［4］提出一种双层定量结构-活性关系（quantitative structure-activity relationship，QSAR）的数学模型，发现二层次的预测网络比单层次的预测精度更高.

随着计算机数据挖掘技术的不断发展，机器学习因其对高维数据深层次有较好的表达能力，已成为药物化学领域用来解决复杂化合物分类问题和预测新的活性分子的一个重要方法［5-6］.黄斌［7］建立了基于支持向量机（support vector machine，SVM）的药物透血脑屏障活性分类模型.SARDARI等［8］为了确定对分枝杆菌有效的化合物，首先采用k近邻算法（k-nearest neighbor，kNN）和k-means 聚类方法对化合物的分子描述符进行聚类，再利用人工神经网络（artificial neural network，ANN）建模方法建立了最小抑菌浓度预测模型.DUTT等［9］通过决策树（decision tree，DT）算法和移动平均分析法预测G蛋白偶联受体-40 激动活性.陆家兴等［10］采用极端梯度提升（extreme gradient boosting，XGBoost）算法开发了一种基于LINCS-L1000扰动信号的药物引导下细胞活性预测算法，并与PCA、局部线性嵌入（locally linear embedding，LLE）、核主成分分析（kernel principal component analysis，KPCA）和独立成分分析（independent component analysis，ICA）进行了比较，证明XGBoost 模型预测效果较优.然而，这些方法使用的都是浅层机器学习算法，且需要海量数据进行训练测试以提高模型的预测精度，面对已知样本和计算单元受到限制的情况，其对复杂问题的泛化能力难以满足预测精度的要求，无法学习更有用的特征.为解决机器学习算法中因存在超参数导致的预测时间久，准确率较低的问题，BERGSTRA等［11］指出当超参数规模较大时，随机搜索（random search，RS）在超参数寻优方面更高效.李玉娟［12］提出一种改进粒子群算法的深度学习超参数优化方法，明显提升了原算法的收敛速度和搜索能力.WU等［13］提出一种基于高斯过程的贝叶斯超参数优化算法，只要少量样本就可找到最优值.但是，上述文献并未对小样本量且数据特征不明显的数据集展开研究.

本研究通过对抗癌候选药物化合物的分子描述符进行统计分析，消除变量间的多重共线性，采用随机森林（random forest，RF）方法［15］对所有潜在影响因素进行重要度筛选，确定作用于乳腺癌治疗靶标雌激素受体α 亚型（estrogen receptor alpha subtype，ERα）生物活性的关键影响因素，利用现有研究中机器学习方法具有较高预测精度的特点，同时考虑机器学习算法中超参数的设置对预测精度的影响，采用贝叶斯超参数优化（Bayesian hyperparametric optimization，BHO）方法，使用双向长短期记忆（bi-directional long short-term memory network，Bi-LSTM）神经网络［16-17］，改进传统长短期记忆（long short-term memory，LSTM）网络信息单向传递的不足，借助卷积神经网络（convolutional neural networks，CNN）［18］逐一提取分子描述符显著变量的关键特征，获取丰富的状态信息，建立ERα 生物活性集成机器学习定量预测模型（multivariate-Bayesian hyperparametric optimized bi-directional long short-term memory，Mul-BHO-Bi-LSTM），并与基于集成学习的梯度提升决策树（gradient boosting decision tree，GBDT）［20］预测模型结果进行对比验证.

1 实验方法

针对乳腺癌治疗靶标ERα，首先收集一系列作用于该靶标的化合物生物活性数据，然后以诸多分子结构描述符作为自变量，选取生物活性值作为因变量，构建ERα生物活性定量预测模型.预测模型的构建如图1，具体步骤如下：

图1 ERα生物活性定量预测模型的构建过程Fig.1 Construction of an quantitative prediction model for ERα biological activity.

1）首先对实验数据进行描述性统计分析，剔除实验数据全部为0的分子描述符（自变量），然后诊断自变量间的多重共线性，逐个剔除不显著的自变量，获得模型的输入变量.

2）构建RF 回归模型，测度自变量的重要性，依据重要性评分（variable importance measure，VIM，记为MVIM）对变量重要性进行筛选，获得MVIM＞0.01的变量共计20个.

3）充分考虑多类别分子描述符的理化性质及拓扑结构特征，采用BHO方法调整Bi-LSTM模型的结构参量，建立ERα 生物活性集成机器学习定量预测模型Mul-BHO-Bi-LSTM.同时，构建集成预测模型GBDT作为对比模型.

4）将MVIM＞0.01的20个变量输入到所构建的Mul-BHO-Bi-LSTM 与GBDT 模型中，经过训练和测试，对两模型的预测结果进行分析评价.

2 数据分析及处理

2.1 数据介绍

实验所用数据源自公开数据集“2021年中国研究生数学建模竞赛”.该数据集包含：①生物活性值：1 974 个化合物作用于ERα 的生物活性值PIC50，该值通常用于描述化合物生物活性，能够拮抗ERα 活性的化合物可能是治疗乳腺癌的候选药物.②分子描述符，包括1 974 个化合物的729 个分子描述符信息.分子描述符是一系列用于描述化合物结构和性质特征的参数，包括物理化学性质（如分子量和LogP等）与拓扑结构特征（如氢键供体数量和氢键受体数量等）等.

本研究通过统计发现，该数据集中化合物样本量较小，PIC50分布无明显规律，分子描述符数据多且杂，同时具有异常值较多、有效特征不明显、难以直接准确预测等特点.因此，需先对数据做深层次的处理分析，挖掘最影响乳腺癌治疗靶标ERα的PIC50变量，为构建化合物的定量结构-活性关系模型提供数据基础.

2.2 描述性统计分析

对1 974 个化合物的ERα 的生物活性数据及729 个分子描述符数据进行描述性统计分析，结果如表1.由表1 可见，分子描述符的统计值整体差异较大，数据结构规律不明显，难以直接进行分析预测，且有225个自变量所对应的样本数据全部为0，为减少数据冗余，删除这225个自变量数据.

表1 描述性统计分析结果（N=1 974）Table 1 Results of descriptive statistical analysis（N=1 974）

2.3 多重共线性诊断

为提升模型预测准确率，需对1 974 个化合物的504 个分子描述符（自变量）进行多重共线性诊断［14］，并逐个检验选入的变量.如此，得到线性回归模型输出结果自变量和因变量之间的相关系数R=0.89；拟合优度指标R2=0.79；调整后R2=0.75；R2的变化量ΔR2=0.7；标准误差为0.72；回归方程显著性的统计量F的变化量ΔF=19.25，与ΔF关联的自由度（degree of freedom，DF）的分子Df1=319；与ΔF关联的自由度的分母Df2=1 652.表2 为方差分析结果.由表2 可见，在显著性水平为0.05 的情况下，回归模型的整体线性关系显著成立.因此，本研究最终剔除AMR、apol、naAromAtom 和nAromBond 等185 个变量，余下的nAcid、AlogP、ALogp2 和ATSc2 等319 个变量为模型的影响自变量.

表2 方差分析Table 2 Analysis of variance

2.4 变量筛选

由于共线性诊断后的变量仍具有维度高、信息冗余复杂等特征，考虑到物理化学性质和拓扑结构特征等参数的异质性，为实现生物活性PIC50指标特征的有效提取，保证能够及时高效识别分子描述符参数，避免因为丰富的数据维度与属性造成的“维数灾难”，因此在建立模型过程中，仍需对变量做进一步筛选.

考虑到RF 模型在特征变量选择方面具有模型精度高、效率快的优点［15］，且能反向评估各输入参数对目标值的相对重要性，对于变量间的多元共线性具有一定的包容性，同时对噪声数据和非平衡样本数据的鲁棒性较好，符合本研究的数据类型及模型构建的思想，因此，本研究采用RF 回归方法进行变量重要度筛选.

经描述性统计分析及多重共线性诊断后，将余下的变量数据输入至RF模型中并进行重要度分析，最终评价结果如图2.其中，变量VC-5、maxssO、MDEC-23、minsOH 和BCUTc-1l 的重要度指数MVIM超过0.01，即重要性相对较高；其余变量的MVIM值均小于0.01.记MVIM＜0.01 的变量为对生物活性重要度较低的变量.筛选出MVIM＞0.01 的变量共计20 个，分别为：ATSc3、BCUTc-1l、BCUTc-1h、C1SP2、VC-5、CrippenLogP、ndssC、SHBint10、ShsOH、minHBint5、minHsOH、minsssN、minsOH、maxHsOH、maxssO、nHBAcc、MDEC-23、MDEO-12、MLFER_A和XlogP.

图2 随机森林模型特征变量筛选结果Fig.2 The screening results of the characteristic variables by the random forest model.

3 ERα生物活性定量预测模型

在深度学习模型中，循环神经网络（recurrent neural network，RNN）及其变体具有良好的非线性拟合特性与深层特征挖掘能力，在定量数据建模方面表现优异，可应用于类别不一、拓扑特征复杂的分子描述符特征建模.但随着序列长度的增加，RNN梯度消失问题愈发严重，特征提取效率也会大幅降低.LSTM也是提取序列特征的神经网络［16-17］，能有效解决梯度消失问题，提高序列特征提取效率，但单向的LSTM 网络难以全面捕捉变量信息等特征，因此，本研究采用性能更优的Bi-LSTM网络模型作为基本预测模型，将筛选出的20 个变量输入至预测模型中，借助贝叶斯网络进行超参数调优，建立ERα生物活性集成机器学习定量预测模型.

3.1 基于Mul-BHO-Bi-LSTM的定量预测模型

考虑到分子描述信息具有非线性和不确定性等特征，本研究首先利用CNN［18］二维遍历自变量特征，再结合BHO 方法［19］，构建基于Bi-LSTM 模型的集成学习模型.

使用CNN 逐一提取化合物分子描述符显著变量的多类关键信息特征，构造二维特征矩阵，使模型更适用于高维多属性PIC50的定量预测.具体过程如下.

3.1.1 构建特征矩阵

考虑到ERα 生物活性定量预测模型的评价指标为生物活性，以及不同分子描述符信息的有效性，本研究构建包含“参数-权重”信息的二维特征矩阵K作为模型输入，即

其中，xi，j为第i个化合物的第j个特征信息所占权重，i=1，2，…，1 974，j=1，2，…，n；为第i个化合物的生物活性值PIC50.

3.1.2 遍历特征矩阵

CNN是一种基于多层监督学习的典型的深度学习网络模型，可用来处理类似上述网格结构的数据.常用的CNN模型有AlexNet、VGG16、GoogLeNet 和LeNet 等［18］，已被广泛用于图像分类、语音识别和自然语言处理等领域.在处理“参数-权重”信息时，CNN可以很好地识别出数据的简单空间模式，并据此在更高级的层中生成更复杂的模式.CNN 主要由卷积层、池化层和全连接层构成，3个级联层描述为

其中，cl，j为第l层卷积层第j个神经元的特征输出；xl-1，i为第l-1层卷积层（激活层）第i个神经元的特征输入；wl，ij为第l层第j个神经元和上一层第i个神经元之间的权重；b为偏置项；xl，j为激活层的输出；l为网络层序号，i和j为神经元序号，φ(cl，j)为非线性的激活函数；pool(xl，j)为池化函数.

生物活性值PIC50具有随机性，因此在构建特征矩阵之后，将特征矩阵K输入到CNN中，获取分子描述符信息数据的权重属性特征，即横向遍历二维特征矩阵K.采用CNN提取数据属性特征后，还需提取状态特征，本研究采用LSTM 网络采集状态特征，即纵向遍历二维特征矩阵.

3.1.3 BHO

Bi-LSTM 需调整每层的神经元数和学习率等参数以保证算法获得最优性能.考虑到Mul-Bi-LSTM模型的训练时间较长，人工设置训练参数对预测性能的影响较大，引入BHO 方法［19］，利用其迭代次数少，运行速度快的优势，实现更有效地搜索可能的超参数空间，优化预测模型性能.BHO算法主要由概率代理模型和采集函数组成，算法目标x*可表示为产生最小均方根误差（root mean square error，RMSE）值和最小损失值的超参数组合，即

其中，f(x)为判断输入x表现优劣的一种测度，即在验证集上评估的RMSE 值和损失值；x为n维超参数决策向量；Ω为超参数决策域空间.

3.2 GBDT集成学习预测模型

生物活性值PIC50 预测属于非线性复杂问题，但线性模型或单一模型可能忽略生物活性的性质影响因素.集成学习通过把多种学习器组合在一起，协作形成一个具有更强学习性能的集成学习器，再将训练好的模型以不同的方式进行融合，从而实现更优的基于分子描述信息预测生物活性能力.本研究拟采用GBDT方法深度挖掘分子描述信息与生物活性值PIC50变量之间的隐性关系.

GBDT 模型［20］是一种分类回归树集成算法，其原理为通过寻找最佳划分特征，进而学习样本路径实现分类，是一种模型复杂度较高、参数随机性较强的学习器.GBDT 模型的基础是对决策树中回归树的迭代优化，基于梯度提升学习策略，在每次迭代时通过最小化损失函数L(xi，f(xi))，在减少残差的梯度方向新建立1棵弱决策树，最后将所有树的结论累加起来得到最终预测结果.记D={(x1，y1)，(x2，y2)，…，(xN，yN)}为包含N个训练样本的数据集，每个样本由d个特征属性描述，即xi=[xi1，xi2，…，xid]，建模过程如下.

1）初始化学习器.估计一个使损失函数极小化的常数值，构建只有一个根节点的树为

其中，c为估计使损失函数极小化的常数值；L(yi，c)为损失函数，i为样本量，i=1，2，…，N.本研究建立的GBDT 模型采用均方误差损失函数作为损失函数，即

其中，f(xi)为机器学习模型预测值.

2）迭代.构建M棵树，设迭代次数m=1，2，…，M，对样本i=1，2，…，N，计算损失函数的负梯度为

根据负梯度方向(xi，γmi)拟合第m棵回归树，得到m棵由J个叶子节点组成的决策树，其对应的叶子节点区域为Rmj，j=1，2，…，J，计算每个叶子节点的最佳残差拟合值，使得损失函数极小化为

其中，cm j为第j个叶子节点的yi与fm-1(xi)之间的最小误差；c为迭代轮数为m时的节点残差拟合值；yi为第j个叶子节点的样本xi观测值；fm-1(xi)为第j个叶子节点的样本xi在上一棵树上的预测值.

更新本轮机器学习模型的预测值为

3）经过M轮迭代，直到达到所预期的基学习器个数，得到最终的强化学习器为

4 模型验证与测试

4.1 评价指标

为评估模型的预测性能，以6种误差评价指标为切入点，建立融合关联指标与误差指标的评价方案.

1）关联指标.采用拟合优度R2与秩相关程度r来判断训练集预测ERα 生物活性PIC50与原始ERα生物活性PIC50之间的拟合程度，R2，r∈［0，1］，且二者越接近1，表明预测ERα生物活性PIC50与原始ERα生物活性PIC50之间拟合程度越好.

2）误差指标.采用误差平均值（error mean）和误差标准差（error standard deviation，记为error std）来评价误差的整体水平和离散程度，误差标准差越大，表明误差数据越离散；以均方误差（meansquare error，MSE）和归一化均方根误差（normalized root mean square error，NRMSE）作为度量指标判断预测结果的精度，两种误差的值越小，说明生物活性PIC50的预测精度越高.

4.2 模型预测效果分析

本研究将构建的Mul-BHO-Bi-LSTM集成模型与经典集成模型GBDT 相对比，统计2 个模型的预测生物活性PIC50数据，以R2、r作为判断Mul-BHOBi-LSTM模型预测值与实际值交汇点拟合程度的指标；以MSE、NRMSE、误差平均值和误差标准差作为误差度量指标判断模型的匹配程度，结果如表3.由表3可见：

表3 Mul-BHO-Bi-LSTM模型与GBDT模型ERα生物活性定量预测结果1）Table 3 Quantitative prediction results of ERα biological activity between Mul-BHO-Bi-LSTM model and GBDT model

1）Mul-BHO-Bi-LSTM 模型与GBDT 模型的综合性整体来看差异较大，说明采用贝叶斯超参数优化的深度学习模型结构在生物活性PIC50预测方面比GBDT集成学习预测模型更具优势，表现出更好的预测精度.

2）拟合程度.Mul-BHO-Bi-LSTM模型的R2和r较GBDT 模型整体提升了31.87%和19.95%，表明前者的预测值与原始值的相关性更高，即PIC50预测值更加接近实际值.

3）匹配程度.Mul-BHO-Bi-LSTM模型较GBDT模型的error mean、error std、MSE 和NRMSE 分别提升98.56%、30.38%、98.82%和88.95%，说明Mul-BHO-Bi-LSTM模型性能更稳定，鲁棒性更佳.

为进一步直观讨论Mul-BHO-Bi-LSTM模型的性能，绘制R2、整体误差与逐样本误差分布和r对比图，结果如图3、图4和图5.

1）拟合优度.Mul-BHO-Bi-LSTM 模型R2值为99.25%，结合图3 可见，PIC50预测值与实际值的交汇点（蓝色圆圈）紧密分布于Y=T（Y为实际值，T为预测值）上下，说明Mul-BHO-Bi-LSTM 模型的预测值十分接近实际值，预测效果理想.

图3 Mul-BHO-Bi-LSTM模型拟合优度Fig.3 Goodness of fit of Mul-BHO-Bi-LSTM model.Circles are data points，solid line is fitting curve，dash-dotted line means predicted value equals to actual value.

2）秩相关程度.Mul-BHO-Bi-LSTM模型r值达99.55%，体现为图4中实际值与预测值两条线基本重合，说明从实际数据与预测数据分布层面来看，模型表现出良好的预测性能，适用于多类分子描述变量预测生物活性PIC50问题.

图4 Mul-BHO-Bi-LSTM模型秩相关性Fig.4 Mul-BHO-Bi-LSTM model rank correlation.Blue line is actual value and the red line is predicted value.

3）预测误差.如图5，模型误差分布整体相对集中，但也不乏个别异常值.但总体来看，在对PIC50进行预测时，Mul-BHO-Bi-LSTM 模型的error mean 为0.002 5，error std 为0.122 6，MSE 值 为0.014 9，NRMSE值为0.018 7，误差相关指标均小于0.15.

图5 Mul-BHO-Bi-LSTM模型PIC50预测误差分布Fig.5 PIC50 prediction error distribution of Mul-BHO-Bi-LSTM model.

由图3 至图5 可见，数据驱动的方法虽可有效解决生物活性预测问题，但在面对实际中小样本、多特征、随机复杂的条件时，集成模型通常需要进行大量的参数调节实验来提升模型的结构性能，因此，当面对已知样本和计算单元受限制的情况时，此类模型对复杂问题的泛化能力及计算能力同样难以满足要求，难以学习更有用的特征.综合上述指标，说明本研究构建参数自寻优Mul-BHO-Bi-LSTM模型的总体效果良好，可用作化合物的定量结构-活性关系模型，具有较好鲁棒性和泛化性.

4.3 优化效果分析

为对Mul-BHO-Bi-LSTM模型的贝叶斯超参数自动寻优效果进行评价，令模型迭代45 次，得到贝叶斯调参可视化结果如图6.由图6 可见，贝叶斯超参数优化器迭代至第6、12和19次时，预测误差骤减约5%、3%和13%；迭代至第19次时，得到最小误差及最佳超参数点；迭代至39 次时，模型收敛，说明贝叶斯优化在目标函数评估成本高的任务时，超参数可在迭代中快速收敛，同时也说明Mul-BHO-Bi-LSTM模型可有效应用于ERα生物活性PIC50的预测.

图6 贝叶斯超参数优化器超参数优化过程Fig.6 Hyperparameter optimization process of BHO.The black circle is minimum classification error of observation.The red squre is estimated minimum classification error.The red cirde is minimum error superparameter.The green squre is optimal hyperparamder point.

结语

针对乳腺癌治疗靶标生物活性预测实际面临的多特征条件筛选问题，考虑到化合物的物理化学性质、拓扑结构特征等参数异质性，为实现生物活性PIC50指标特征的有效提取，保证能够及时且高效地识别分子描述符参数，采用RF 回归方法，测度分子描述符变量的重要度，有效筛选出生物活性预测所需的关键特征参数.

为实现ERα 生物活性的精准预测，基于CNN模型的二维特征矩阵输入结果，采用贝叶斯算法和Adam 优化算法对Bi-LSTM 模型进行超参数寻优.采用贝叶斯算法优化超参数后，集成机器学习预测模型Mul-BHO-Bi-LSTM可较快收敛，有效解决浅层机器学习存在的局部最小化和过拟合等问题.与经典GBDT 模型相比，Mul-BHO-Bi-LSTM 预测模型表现出较高的预测精度与性能，能改善混合模型调参导致计算时间较久、精度不高的缺陷，能够适用于小样本、多特征条件下分子描述变量预测生物活性PIC50问题.本研究方法可进一步拓展，应用于图像识别、自然语言处理和语音识别等领域.