抗癌候選藥物ERα抑制劑活性預(yù)測

夏玉蘭，謝濟(jì)銘，王雅婧，盧夢媛，王錦銳，秦雅琴

1）昆明理工大學(xué)交通工程學(xué)院，云南昆明 650504；2）溫州醫(yī)科大學(xué)第一臨床醫(yī)學(xué)院，浙江溫州325006

乳腺癌是目前世界上最常見且致死率較高的癌癥之一.根據(jù)世界衛(wèi)生組織（Wold Health Organization，WHO）國際癌癥研究機構(gòu)（International Agency for Research on Cancer，IARC）2020年發(fā)布的全球癌癥負(fù)擔(dān)數(shù)據(jù)，乳腺癌新發(fā)病例及死亡率均列首位，成為全球“第一癌”［1］.然而，現(xiàn)有的乳腺癌臨床治療方法不能完全阻止乳腺癌的復(fù)發(fā)及轉(zhuǎn)移擴(kuò)散.因此，研發(fā)新型抗乳腺癌藥物是治療乳腺癌和改善乳腺癌患者預(yù)后的有效策略［2］.

在抗乳腺癌藥物研發(fā)過程中，可通過對藥物化合物分子結(jié)構(gòu)特征、性質(zhì)進(jìn)行研究，給合數(shù)理統(tǒng)計的方法建立定量構(gòu)效關(guān)系模型來實現(xiàn)藥物活性預(yù)測.KIDERA 等［3］運用主成分分析（principal component analysis，PCA）從氨基酸的188 個特征參數(shù)中篩選出10 個最佳主成分變量，并采用多元線性回歸對多肽藥物進(jìn)行定量表征與活性預(yù)測，取得了較好的結(jié)果.王青艷等［4］提出一種雙層定量結(jié)構(gòu)-活性關(guān)系（quantitative structure-activity relationship，QSAR）的數(shù)學(xué)模型，發(fā)現(xiàn)二層次的預(yù)測網(wǎng)絡(luò)比單層次的預(yù)測精度更高.

隨著計算機數(shù)據(jù)挖掘技術(shù)的不斷發(fā)展，機器學(xué)習(xí)因其對高維數(shù)據(jù)深層次有較好的表達(dá)能力，已成為藥物化學(xué)領(lǐng)域用來解決復(fù)雜化合物分類問題和預(yù)測新的活性分子的一個重要方法［5-6］.黃斌［7］建立了基于支持向量機（support vector machine，SVM）的藥物透血腦屏障活性分類模型.SARDARI等［8］為了確定對分枝桿菌有效的化合物，首先采用k近鄰算法（k-nearest neighbor，kNN）和k-means 聚類方法對化合物的分子描述符進(jìn)行聚類，再利用人工神經(jīng)網(wǎng)絡(luò)（artificial neural network，ANN）建模方法建立了最小抑菌濃度預(yù)測模型.DUTT等［9］通過決策樹（decision tree，DT）算法和移動平均分析法預(yù)測G蛋白偶聯(lián)受體-40 激動活性.陸家興等［10］采用極端梯度提升（extreme gradient boosting，XGBoost）算法開發(fā)了一種基于LINCS-L1000擾動信號的藥物引導(dǎo)下細(xì)胞活性預(yù)測算法，并與PCA、局部線性嵌入（locally linear embedding，LLE）、核主成分分析（kernel principal component analysis，KPCA）和獨立成分分析（independent component analysis，ICA）進(jìn)行了比較，證明XGBoost 模型預(yù)測效果較優(yōu).然而，這些方法使用的都是淺層機器學(xué)習(xí)算法，且需要海量數(shù)據(jù)進(jìn)行訓(xùn)練測試以提高模型的預(yù)測精度，面對已知樣本和計算單元受到限制的情況，其對復(fù)雜問題的泛化能力難以滿足預(yù)測精度的要求，無法學(xué)習(xí)更有用的特征.為解決機器學(xué)習(xí)算法中因存在超參數(shù)導(dǎo)致的預(yù)測時間久，準(zhǔn)確率較低的問題，BERGSTRA等［11］指出當(dāng)超參數(shù)規(guī)模較大時，隨機搜索（random search，RS）在超參數(shù)尋優(yōu)方面更高效.李玉娟［12］提出一種改進(jìn)粒子群算法的深度學(xué)習(xí)超參數(shù)優(yōu)化方法，明顯提升了原算法的收斂速度和搜索能力.WU等［13］提出一種基于高斯過程的貝葉斯超參數(shù)優(yōu)化算法，只要少量樣本就可找到最優(yōu)值.但是，上述文獻(xiàn)并未對小樣本量且數(shù)據(jù)特征不明顯的數(shù)據(jù)集展開研究.

本研究通過對抗癌候選藥物化合物的分子描述符進(jìn)行統(tǒng)計分析，消除變量間的多重共線性，采用隨機森林（random forest，RF）方法［15］對所有潛在影響因素進(jìn)行重要度篩選，確定作用于乳腺癌治療靶標(biāo)雌激素受體α 亞型（estrogen receptor alpha subtype，ERα）生物活性的關(guān)鍵影響因素，利用現(xiàn)有研究中機器學(xué)習(xí)方法具有較高預(yù)測精度的特點，同時考慮機器學(xué)習(xí)算法中超參數(shù)的設(shè)置對預(yù)測精度的影響，采用貝葉斯超參數(shù)優(yōu)化（Bayesian hyperparametric optimization，BHO）方法，使用雙向長短期記憶（bi-directional long short-term memory network，Bi-LSTM）神經(jīng)網(wǎng)絡(luò)［16-17］，改進(jìn)傳統(tǒng)長短期記憶（long short-term memory，LSTM）網(wǎng)絡(luò)信息單向傳遞的不足，借助卷積神經(jīng)網(wǎng)絡(luò)（convolutional neural networks，CNN）［18］逐一提取分子描述符顯著變量的關(guān)鍵特征，獲取豐富的狀態(tài)信息，建立ERα 生物活性集成機器學(xué)習(xí)定量預(yù)測模型（multivariate-Bayesian hyperparametric optimized bi-directional long short-term memory，Mul-BHO-Bi-LSTM），并與基于集成學(xué)習(xí)的梯度提升決策樹（gradient boosting decision tree，GBDT）［20］預(yù)測模型結(jié)果進(jìn)行對比驗證.

1 實驗方法

針對乳腺癌治療靶標(biāo)ERα，首先收集一系列作用于該靶標(biāo)的化合物生物活性數(shù)據(jù)，然后以諸多分子結(jié)構(gòu)描述符作為自變量，選取生物活性值作為因變量，構(gòu)建ERα生物活性定量預(yù)測模型.預(yù)測模型的構(gòu)建如圖1，具體步驟如下：

圖1 ERα生物活性定量預(yù)測模型的構(gòu)建過程Fig.1 Construction of an quantitative prediction model for ERα biological activity.

1）首先對實驗數(shù)據(jù)進(jìn)行描述性統(tǒng)計分析，剔除實驗數(shù)據(jù)全部為0的分子描述符（自變量），然后診斷自變量間的多重共線性，逐個剔除不顯著的自變量，獲得模型的輸入變量.

2）構(gòu)建RF 回歸模型，測度自變量的重要性，依據(jù)重要性評分（variable importance measure，VIM，記為MVIM）對變量重要性進(jìn)行篩選，獲得MVIM＞0.01的變量共計20個.

3）充分考慮多類別分子描述符的理化性質(zhì)及拓?fù)浣Y(jié)構(gòu)特征，采用BHO方法調(diào)整Bi-LSTM模型的結(jié)構(gòu)參量，建立ERα 生物活性集成機器學(xué)習(xí)定量預(yù)測模型Mul-BHO-Bi-LSTM.同時，構(gòu)建集成預(yù)測模型GBDT作為對比模型.

4）將MVIM＞0.01的20個變量輸入到所構(gòu)建的Mul-BHO-Bi-LSTM 與GBDT 模型中，經(jīng)過訓(xùn)練和測試，對兩模型的預(yù)測結(jié)果進(jìn)行分析評價.

2 數(shù)據(jù)分析及處理

2.1 數(shù)據(jù)介紹

實驗所用數(shù)據(jù)源自公開數(shù)據(jù)集“2021年中國研究生數(shù)學(xué)建模競賽”.該數(shù)據(jù)集包含：①生物活性值：1 974 個化合物作用于ERα 的生物活性值PIC50，該值通常用于描述化合物生物活性，能夠拮抗ERα 活性的化合物可能是治療乳腺癌的候選藥物.②分子描述符，包括1 974 個化合物的729 個分子描述符信息.分子描述符是一系列用于描述化合物結(jié)構(gòu)和性質(zhì)特征的參數(shù)，包括物理化學(xué)性質(zhì)（如分子量和LogP等）與拓?fù)浣Y(jié)構(gòu)特征（如氫鍵供體數(shù)量和氫鍵受體數(shù)量等）等.

本研究通過統(tǒng)計發(fā)現(xiàn)，該數(shù)據(jù)集中化合物樣本量較小，PIC50分布無明顯規(guī)律，分子描述符數(shù)據(jù)多且雜，同時具有異常值較多、有效特征不明顯、難以直接準(zhǔn)確預(yù)測等特點.因此，需先對數(shù)據(jù)做深層次的處理分析，挖掘最影響乳腺癌治療靶標(biāo)ERα的PIC50變量，為構(gòu)建化合物的定量結(jié)構(gòu)-活性關(guān)系模型提供數(shù)據(jù)基礎(chǔ).

2.2 描述性統(tǒng)計分析

對1 974 個化合物的ERα 的生物活性數(shù)據(jù)及729 個分子描述符數(shù)據(jù)進(jìn)行描述性統(tǒng)計分析，結(jié)果如表1.由表1 可見，分子描述符的統(tǒng)計值整體差異較大，數(shù)據(jù)結(jié)構(gòu)規(guī)律不明顯，難以直接進(jìn)行分析預(yù)測，且有225個自變量所對應(yīng)的樣本數(shù)據(jù)全部為0，為減少數(shù)據(jù)冗余，刪除這225個自變量數(shù)據(jù).

表1 描述性統(tǒng)計分析結(jié)果（N=1 974）Table 1 Results of descriptive statistical analysis（N=1 974）

2.3 多重共線性診斷

為提升模型預(yù)測準(zhǔn)確率，需對1 974 個化合物的504 個分子描述符（自變量）進(jìn)行多重共線性診斷［14］，并逐個檢驗選入的變量.如此，得到線性回歸模型輸出結(jié)果自變量和因變量之間的相關(guān)系數(shù)R=0.89；擬合優(yōu)度指標(biāo)R2=0.79；調(diào)整后R2=0.75；R2的變化量ΔR2=0.7；標(biāo)準(zhǔn)誤差為0.72；回歸方程顯著性的統(tǒng)計量F的變化量ΔF=19.25，與ΔF關(guān)聯(lián)的自由度（degree of freedom，DF）的分子Df1=319；與ΔF關(guān)聯(lián)的自由度的分母Df2=1 652.表2 為方差分析結(jié)果.由表2 可見，在顯著性水平為0.05 的情況下，回歸模型的整體線性關(guān)系顯著成立.因此，本研究最終剔除AMR、apol、naAromAtom 和nAromBond 等185 個變量，余下的nAcid、AlogP、ALogp2 和ATSc2 等319 個變量為模型的影響自變量.

表2 方差分析Table 2 Analysis of variance

2.4 變量篩選

由于共線性診斷后的變量仍具有維度高、信息冗余復(fù)雜等特征，考慮到物理化學(xué)性質(zhì)和拓?fù)浣Y(jié)構(gòu)特征等參數(shù)的異質(zhì)性，為實現(xiàn)生物活性PIC50指標(biāo)特征的有效提取，保證能夠及時高效識別分子描述符參數(shù)，避免因為豐富的數(shù)據(jù)維度與屬性造成的“維數(shù)災(zāi)難”，因此在建立模型過程中，仍需對變量做進(jìn)一步篩選.

考慮到RF 模型在特征變量選擇方面具有模型精度高、效率快的優(yōu)點［15］，且能反向評估各輸入?yún)?shù)對目標(biāo)值的相對重要性，對于變量間的多元共線性具有一定的包容性，同時對噪聲數(shù)據(jù)和非平衡樣本數(shù)據(jù)的魯棒性較好，符合本研究的數(shù)據(jù)類型及模型構(gòu)建的思想，因此，本研究采用RF 回歸方法進(jìn)行變量重要度篩選.

經(jīng)描述性統(tǒng)計分析及多重共線性診斷后，將余下的變量數(shù)據(jù)輸入至RF模型中并進(jìn)行重要度分析，最終評價結(jié)果如圖2.其中，變量VC-5、maxssO、MDEC-23、minsOH 和BCUTc-1l 的重要度指數(shù)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 隨機森林模型特征變量篩選結(jié)果Fig.2 The screening results of the characteristic variables by the random forest model.

3 ERα生物活性定量預(yù)測模型

在深度學(xué)習(xí)模型中，循環(huán)神經(jīng)網(wǎng)絡(luò)（recurrent neural network，RNN）及其變體具有良好的非線性擬合特性與深層特征挖掘能力，在定量數(shù)據(jù)建模方面表現(xiàn)優(yōu)異，可應(yīng)用于類別不一、拓?fù)涮卣鲝?fù)雜的分子描述符特征建模.但隨著序列長度的增加，RNN梯度消失問題愈發(fā)嚴(yán)重，特征提取效率也會大幅降低.LSTM也是提取序列特征的神經(jīng)網(wǎng)絡(luò)［16-17］，能有效解決梯度消失問題，提高序列特征提取效率，但單向的LSTM 網(wǎng)絡(luò)難以全面捕捉變量信息等特征，因此，本研究采用性能更優(yōu)的Bi-LSTM網(wǎng)絡(luò)模型作為基本預(yù)測模型，將篩選出的20 個變量輸入至預(yù)測模型中，借助貝葉斯網(wǎng)絡(luò)進(jìn)行超參數(shù)調(diào)優(yōu)，建立ERα生物活性集成機器學(xué)習(xí)定量預(yù)測模型.

3.1 基于Mul-BHO-Bi-LSTM的定量預(yù)測模型

考慮到分子描述信息具有非線性和不確定性等特征，本研究首先利用CNN［18］二維遍歷自變量特征，再結(jié)合BHO 方法［19］，構(gòu)建基于Bi-LSTM 模型的集成學(xué)習(xí)模型.

使用CNN 逐一提取化合物分子描述符顯著變量的多類關(guān)鍵信息特征，構(gòu)造二維特征矩陣，使模型更適用于高維多屬性PIC50的定量預(yù)測.具體過程如下.

3.1.1 構(gòu)建特征矩陣

考慮到ERα 生物活性定量預(yù)測模型的評價指標(biāo)為生物活性，以及不同分子描述符信息的有效性，本研究構(gòu)建包含“參數(shù)-權(quán)重”信息的二維特征矩陣K作為模型輸入，即

其中，xi，j為第i個化合物的第j個特征信息所占權(quán)重，i=1，2，…，1 974，j=1，2，…，n；為第i個化合物的生物活性值PIC50.

3.1.2 遍歷特征矩陣

CNN是一種基于多層監(jiān)督學(xué)習(xí)的典型的深度學(xué)習(xí)網(wǎng)絡(luò)模型，可用來處理類似上述網(wǎng)格結(jié)構(gòu)的數(shù)據(jù).常用的CNN模型有AlexNet、VGG16、GoogLeNet 和LeNet 等［18］，已被廣泛用于圖像分類、語音識別和自然語言處理等領(lǐng)域.在處理“參數(shù)-權(quán)重”信息時，CNN可以很好地識別出數(shù)據(jù)的簡單空間模式，并據(jù)此在更高級的層中生成更復(fù)雜的模式.CNN 主要由卷積層、池化層和全連接層構(gòu)成，3個級聯(lián)層描述為

其中，cl，j為第l層卷積層第j個神經(jīng)元的特征輸出；xl-1，i為第l-1層卷積層（激活層）第i個神經(jīng)元的特征輸入；wl，ij為第l層第j個神經(jīng)元和上一層第i個神經(jīng)元之間的權(quán)重；b為偏置項；xl，j為激活層的輸出；l為網(wǎng)絡(luò)層序號，i和j為神經(jīng)元序號，φ(cl，j)為非線性的激活函數(shù)；pool(xl，j)為池化函數(shù).

生物活性值PIC50具有隨機性，因此在構(gòu)建特征矩陣之后，將特征矩陣K輸入到CNN中，獲取分子描述符信息數(shù)據(jù)的權(quán)重屬性特征，即橫向遍歷二維特征矩陣K.采用CNN提取數(shù)據(jù)屬性特征后，還需提取狀態(tài)特征，本研究采用LSTM 網(wǎng)絡(luò)采集狀態(tài)特征，即縱向遍歷二維特征矩陣.

3.1.3 BHO

Bi-LSTM 需調(diào)整每層的神經(jīng)元數(shù)和學(xué)習(xí)率等參數(shù)以保證算法獲得最優(yōu)性能.考慮到Mul-Bi-LSTM模型的訓(xùn)練時間較長，人工設(shè)置訓(xùn)練參數(shù)對預(yù)測性能的影響較大，引入BHO 方法［19］，利用其迭代次數(shù)少，運行速度快的優(yōu)勢，實現(xiàn)更有效地搜索可能的超參數(shù)空間，優(yōu)化預(yù)測模型性能.BHO算法主要由概率代理模型和采集函數(shù)組成，算法目標(biāo)x*可表示為產(chǎn)生最小均方根誤差（root mean square error，RMSE）值和最小損失值的超參數(shù)組合，即

其中，f(x)為判斷輸入x表現(xiàn)優(yōu)劣的一種測度，即在驗證集上評估的RMSE 值和損失值；x為n維超參數(shù)決策向量；Ω為超參數(shù)決策域空間.

3.2 GBDT集成學(xué)習(xí)預(yù)測模型

生物活性值PIC50 預(yù)測屬于非線性復(fù)雜問題，但線性模型或單一模型可能忽略生物活性的性質(zhì)影響因素.集成學(xué)習(xí)通過把多種學(xué)習(xí)器組合在一起，協(xié)作形成一個具有更強學(xué)習(xí)性能的集成學(xué)習(xí)器，再將訓(xùn)練好的模型以不同的方式進(jìn)行融合，從而實現(xiàn)更優(yōu)的基于分子描述信息預(yù)測生物活性能力.本研究擬采用GBDT方法深度挖掘分子描述信息與生物活性值PIC50變量之間的隱性關(guān)系.

GBDT 模型［20］是一種分類回歸樹集成算法，其原理為通過尋找最佳劃分特征，進(jìn)而學(xué)習(xí)樣本路徑實現(xiàn)分類，是一種模型復(fù)雜度較高、參數(shù)隨機性較強的學(xué)習(xí)器.GBDT 模型的基礎(chǔ)是對決策樹中回歸樹的迭代優(yōu)化，基于梯度提升學(xué)習(xí)策略，在每次迭代時通過最小化損失函數(shù)L(xi，f(xi))，在減少殘差的梯度方向新建立1棵弱決策樹，最后將所有樹的結(jié)論累加起來得到最終預(yù)測結(jié)果.記D={(x1，y1)，(x2，y2)，…，(xN，yN)}為包含N個訓(xùn)練樣本的數(shù)據(jù)集，每個樣本由d個特征屬性描述，即xi=[xi1，xi2，…，xid]，建模過程如下.

1）初始化學(xué)習(xí)器.估計一個使損失函數(shù)極小化的常數(shù)值，構(gòu)建只有一個根節(jié)點的樹為

其中，c為估計使損失函數(shù)極小化的常數(shù)值；L(yi，c)為損失函數(shù)，i為樣本量，i=1，2，…，N.本研究建立的GBDT 模型采用均方誤差損失函數(shù)作為損失函數(shù)，即

其中，f(xi)為機器學(xué)習(xí)模型預(yù)測值.

2）迭代.構(gòu)建M棵樹，設(shè)迭代次數(shù)m=1，2，…，M，對樣本i=1，2，…，N，計算損失函數(shù)的負(fù)梯度為

根據(jù)負(fù)梯度方向(xi，γmi)擬合第m棵回歸樹，得到m棵由J個葉子節(jié)點組成的決策樹，其對應(yīng)的葉子節(jié)點區(qū)域為Rmj，j=1，2，…，J，計算每個葉子節(jié)點的最佳殘差擬合值，使得損失函數(shù)極小化為

其中，cm j為第j個葉子節(jié)點的yi與fm-1(xi)之間的最小誤差；c為迭代輪數(shù)為m時的節(jié)點殘差擬合值；yi為第j個葉子節(jié)點的樣本xi觀測值；fm-1(xi)為第j個葉子節(jié)點的樣本xi在上一棵樹上的預(yù)測值.

更新本輪機器學(xué)習(xí)模型的預(yù)測值為

3）經(jīng)過M輪迭代，直到達(dá)到所預(yù)期的基學(xué)習(xí)器個數(shù)，得到最終的強化學(xué)習(xí)器為

4 模型驗證與測試

4.1 評價指標(biāo)

為評估模型的預(yù)測性能，以6種誤差評價指標(biāo)為切入點，建立融合關(guān)聯(lián)指標(biāo)與誤差指標(biāo)的評價方案.

1）關(guān)聯(lián)指標(biāo).采用擬合優(yōu)度R2與秩相關(guān)程度r來判斷訓(xùn)練集預(yù)測ERα 生物活性PIC50與原始ERα生物活性PIC50之間的擬合程度，R2，r∈［0，1］，且二者越接近1，表明預(yù)測ERα生物活性PIC50與原始ERα生物活性PIC50之間擬合程度越好.

2）誤差指標(biāo).采用誤差平均值（error mean）和誤差標(biāo)準(zhǔn)差（error standard deviation，記為error std）來評價誤差的整體水平和離散程度，誤差標(biāo)準(zhǔn)差越大，表明誤差數(shù)據(jù)越離散；以均方誤差（meansquare error，MSE）和歸一化均方根誤差（normalized root mean square error，NRMSE）作為度量指標(biāo)判斷預(yù)測結(jié)果的精度，兩種誤差的值越小，說明生物活性PIC50的預(yù)測精度越高.

4.2 模型預(yù)測效果分析

本研究將構(gòu)建的Mul-BHO-Bi-LSTM集成模型與經(jīng)典集成模型GBDT 相對比，統(tǒng)計2 個模型的預(yù)測生物活性PIC50數(shù)據(jù)，以R2、r作為判斷Mul-BHOBi-LSTM模型預(yù)測值與實際值交匯點擬合程度的指標(biāo)；以MSE、NRMSE、誤差平均值和誤差標(biāo)準(zhǔn)差作為誤差度量指標(biāo)判斷模型的匹配程度，結(jié)果如表3.由表3可見：

表3 Mul-BHO-Bi-LSTM模型與GBDT模型ERα生物活性定量預(yù)測結(jié)果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 模型的綜合性整體來看差異較大，說明采用貝葉斯超參數(shù)優(yōu)化的深度學(xué)習(xí)模型結(jié)構(gòu)在生物活性PIC50預(yù)測方面比GBDT集成學(xué)習(xí)預(yù)測模型更具優(yōu)勢，表現(xiàn)出更好的預(yù)測精度.

2）擬合程度.Mul-BHO-Bi-LSTM模型的R2和r較GBDT 模型整體提升了31.87%和19.95%，表明前者的預(yù)測值與原始值的相關(guān)性更高，即PIC50預(yù)測值更加接近實際值.

3）匹配程度.Mul-BHO-Bi-LSTM模型較GBDT模型的error mean、error std、MSE 和NRMSE 分別提升98.56%、30.38%、98.82%和88.95%，說明Mul-BHO-Bi-LSTM模型性能更穩(wěn)定，魯棒性更佳.

為進(jìn)一步直觀討論Mul-BHO-Bi-LSTM模型的性能，繪制R2、整體誤差與逐樣本誤差分布和r對比圖，結(jié)果如圖3、圖4和圖5.

1）擬合優(yōu)度.Mul-BHO-Bi-LSTM 模型R2值為99.25%，結(jié)合圖3 可見，PIC50預(yù)測值與實際值的交匯點（藍(lán)色圓圈）緊密分布于Y=T（Y為實際值，T為預(yù)測值）上下，說明Mul-BHO-Bi-LSTM 模型的預(yù)測值十分接近實際值，預(yù)測效果理想.

圖3 Mul-BHO-Bi-LSTM模型擬合優(yōu)度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）秩相關(guān)程度.Mul-BHO-Bi-LSTM模型r值達(dá)99.55%，體現(xiàn)為圖4中實際值與預(yù)測值兩條線基本重合，說明從實際數(shù)據(jù)與預(yù)測數(shù)據(jù)分布層面來看，模型表現(xiàn)出良好的預(yù)測性能，適用于多類分子描述變量預(yù)測生物活性PIC50問題.

圖4 Mul-BHO-Bi-LSTM模型秩相關(guān)性Fig.4 Mul-BHO-Bi-LSTM model rank correlation.Blue line is actual value and the red line is predicted value.

3）預(yù)測誤差.如圖5，模型誤差分布整體相對集中，但也不乏個別異常值.但總體來看，在對PIC50進(jìn)行預(yù)測時，Mul-BHO-Bi-LSTM 模型的error mean 為0.002 5，error std 為0.122 6，MSE 值 為0.014 9，NRMSE值為0.018 7，誤差相關(guān)指標(biāo)均小于0.15.

圖5 Mul-BHO-Bi-LSTM模型PIC50預(yù)測誤差分布Fig.5 PIC50 prediction error distribution of Mul-BHO-Bi-LSTM model.

由圖3 至圖5 可見，數(shù)據(jù)驅(qū)動的方法雖可有效解決生物活性預(yù)測問題，但在面對實際中小樣本、多特征、隨機復(fù)雜的條件時，集成模型通常需要進(jìn)行大量的參數(shù)調(diào)節(jié)實驗來提升模型的結(jié)構(gòu)性能，因此，當(dāng)面對已知樣本和計算單元受限制的情況時，此類模型對復(fù)雜問題的泛化能力及計算能力同樣難以滿足要求，難以學(xué)習(xí)更有用的特征.綜合上述指標(biāo)，說明本研究構(gòu)建參數(shù)自尋優(yōu)Mul-BHO-Bi-LSTM模型的總體效果良好，可用作化合物的定量結(jié)構(gòu)-活性關(guān)系模型，具有較好魯棒性和泛化性.

4.3 優(yōu)化效果分析

為對Mul-BHO-Bi-LSTM模型的貝葉斯超參數(shù)自動尋優(yōu)效果進(jìn)行評價，令模型迭代45 次，得到貝葉斯調(diào)參可視化結(jié)果如圖6.由圖6 可見，貝葉斯超參數(shù)優(yōu)化器迭代至第6、12和19次時，預(yù)測誤差驟減約5%、3%和13%；迭代至第19次時，得到最小誤差及最佳超參數(shù)點；迭代至39 次時，模型收斂，說明貝葉斯優(yōu)化在目標(biāo)函數(shù)評估成本高的任務(wù)時，超參數(shù)可在迭代中快速收斂，同時也說明Mul-BHO-Bi-LSTM模型可有效應(yīng)用于ERα生物活性PIC50的預(yù)測.

圖6 貝葉斯超參數(shù)優(yōu)化器超參數(shù)優(yōu)化過程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.

結(jié)語

針對乳腺癌治療靶標(biāo)生物活性預(yù)測實際面臨的多特征條件篩選問題，考慮到化合物的物理化學(xué)性質(zhì)、拓?fù)浣Y(jié)構(gòu)特征等參數(shù)異質(zhì)性，為實現(xiàn)生物活性PIC50指標(biāo)特征的有效提取，保證能夠及時且高效地識別分子描述符參數(shù)，采用RF 回歸方法，測度分子描述符變量的重要度，有效篩選出生物活性預(yù)測所需的關(guān)鍵特征參數(shù).

為實現(xiàn)ERα 生物活性的精準(zhǔn)預(yù)測，基于CNN模型的二維特征矩陣輸入結(jié)果，采用貝葉斯算法和Adam 優(yōu)化算法對Bi-LSTM 模型進(jìn)行超參數(shù)尋優(yōu).采用貝葉斯算法優(yōu)化超參數(shù)后，集成機器學(xué)習(xí)預(yù)測模型Mul-BHO-Bi-LSTM可較快收斂，有效解決淺層機器學(xué)習(xí)存在的局部最小化和過擬合等問題.與經(jīng)典GBDT 模型相比，Mul-BHO-Bi-LSTM 預(yù)測模型表現(xiàn)出較高的預(yù)測精度與性能，能改善混合模型調(diào)參導(dǎo)致計算時間較久、精度不高的缺陷，能夠適用于小樣本、多特征條件下分子描述變量預(yù)測生物活性PIC50問題.本研究方法可進(jìn)一步拓展，應(yīng)用于圖像識別、自然語言處理和語音識別等領(lǐng)域.