基于機(jī)器學(xué)習(xí)方法的ERα 抑制劑活性預(yù)測(cè)

杜雪平*

（湖北工業(yè)大學(xué)理學(xué)院，湖北 武漢 430068）

乳腺癌是世界上最常見且致死率較高的癌癥類型[1]。近十年來我國(guó)的乳腺癌發(fā)病率上升了47 %，發(fā)病率還在逐年增加，且乳腺癌發(fā)病逐漸呈年輕化。

雌激素受體（ER）在乳腺癌發(fā)展過程中起著非常重要的作用，是乳腺癌內(nèi)分泌療法最主要的靶點(diǎn)[2]。ER 分為ERα 和ERβ 兩種亞型[3]，ERα 主要分布在乳房和子宮內(nèi)膜中，ERβ 與神經(jīng)系統(tǒng)和免疫系統(tǒng)有關(guān)。ERα 在正常的乳腺組織中表現(xiàn)水平很低，但在乳腺癌患者的乳腺組織中表達(dá)水平異常增高，因此ERα 被認(rèn)為是治療乳腺癌的重要靶標(biāo)，抑制ERα 受體的活性是治療乳腺癌的重要手段，能夠抑制ERα 活性的化合物可能是治療乳腺癌的候選藥物。通過實(shí)驗(yàn)的方法來高通量篩選化合物費(fèi)時(shí)費(fèi)力，因此可以采用基于計(jì)算的虛擬篩選方法，其中基于機(jī)器學(xué)習(xí)來構(gòu)建化合物的定量結(jié)構(gòu)-活性關(guān)系（Quantitative Structure-Activity Relationship, QSAR）模型是最主流的方法。目前構(gòu)建化合物的QSAR 模型有如下要求[4]：（1）確定的目標(biāo)（化合物生物活性）；（2）明確的算法；（3）確定的應(yīng)用領(lǐng)域；（4）顯著的相關(guān)性、良好的穩(wěn)健性和預(yù)測(cè)能力；（5）模型易于解釋。Dadfar E 等人[5]利用人工神經(jīng)網(wǎng)絡(luò)（ANN）方法建立磺胺類藥物化合物的活性預(yù)測(cè)模型，雖有較好的預(yù)測(cè)能力，但ANN 方法存在黑箱。Kurunczi L 等人[6]在構(gòu)建QSAR 模型時(shí)利用偏最小二乘法（PLS）進(jìn)行變量選擇，Asikainen A H 等人[7]利用k-近鄰（KNN）方法進(jìn)行變量選擇，采用PLS 方法和KNN方法篩選出的變量較多，不易于對(duì)模型進(jìn)行有效解釋。

本文使用方差過濾法和Lasso 回歸對(duì)分子描述符進(jìn)行合理篩選，基于隨機(jī)森林、支持向量機(jī)和多元線性回歸三種機(jī)器學(xué)習(xí)方法構(gòu)建ERα 抑制劑的活性預(yù)測(cè)模型，其中使用隨機(jī)森林具有更好的預(yù)測(cè)能力和穩(wěn)健性。

1 數(shù)據(jù)集劃分與特征篩選

本文數(shù)據(jù)使用“華為杯”第十八屆中國(guó)研究生數(shù)學(xué)建模競(jìng)賽D 題中數(shù)據(jù)，包括1974 個(gè)化合物的729 個(gè)分子描述符和生物活性 pIC50 值。 使用 sklearn.model_selection 模塊中的train_test_split 函數(shù)來將1974個(gè)化合物以4:1 劃分為訓(xùn)練集和測(cè)試集，訓(xùn)練集樣本數(shù)為1579，測(cè)試集樣本數(shù)為395。在訓(xùn)練集上訓(xùn)練模型，再用測(cè)試集的數(shù)據(jù)來考察模型的預(yù)測(cè)效果。

本文數(shù)據(jù)集有729 個(gè)分子描述符，特征維度大，不利于模型的構(gòu)建，因此需要進(jìn)行變量篩選。結(jié)合數(shù)據(jù)集存在特征維度龐大的特點(diǎn)，本文將過濾法[8]與嵌入法[8]相結(jié)合，首先使用方差過濾法對(duì)分子描述符變量進(jìn)行初步篩選，方差過濾法簡(jiǎn)單，能夠快速剔除掉信息量很小的特征變量。再使用Lasso 回歸[9]消除噪聲特征（即對(duì)生物活性值影響很小的特征）和關(guān)聯(lián)特征（即特征之間相關(guān)性較強(qiáng)的特征），不僅能夠保證模型擁有良好的性能，還節(jié)省了大量的處理時(shí)間和計(jì)算能力。特征篩選具體步驟如下：（1）方差過濾法：本文首先基于方差過濾法利用Python 軟件對(duì)數(shù)據(jù)集中729 個(gè)分子描述符進(jìn)行初步篩選，將方差閾值設(shè)定為0.05。對(duì)任一分子描述符，遍歷所有樣本計(jì)算該分子描述符的方差，如果方差小于等于0.05 則將其剔除，即刪除取值變化不明顯的分子描述符，保留方差大于0.05 的分子描述符。經(jīng)過方差過濾法最終在729 個(gè)變量中剔除了369 個(gè)變量，保留了360 個(gè)變量。（2）Lasso 回歸算法：分子描述符經(jīng)過初步篩選之后，再使用Lasso 回歸進(jìn)一步篩選。以化合物活性pIC50 值作為目標(biāo)變量，360 個(gè)分子描述符作為自變量構(gòu)建Lasso 回歸模型，通過對(duì)損失函數(shù)加入懲罰項(xiàng)，使得訓(xùn)練求解參數(shù)過程中會(huì)考慮系數(shù)的大小，通過設(shè)置縮減系數(shù)(懲罰系數(shù)=0.001)，使得影響較小的特征的系數(shù)衰減到0。Lasso回歸系數(shù)代表了分子描述符變量對(duì)生物活性pIC50 值的重要性，Lasso 回歸系數(shù)絕對(duì)值越大，說明分子描述符對(duì)pIC50 值越重要，根據(jù)重要性排序，選擇對(duì)pIC50 值影響最大的50 個(gè)分子描述符。

2 基于機(jī)器學(xué)習(xí)方法對(duì)ERα 抑制劑活性的預(yù)測(cè)

本文分別用隨機(jī)森林、支持向量機(jī)和多元線性回歸等機(jī)器學(xué)習(xí)方法對(duì)ERα 抑制劑的活性進(jìn)行預(yù)測(cè)，并用均方誤差MSE 來評(píng)價(jià)模型預(yù)測(cè)效果。MSE 是預(yù)測(cè)值與真實(shí)值差的平方和的平均，即：

MSE 的范圍是[0,+∞)，當(dāng)預(yù)測(cè)值與真實(shí)值完全相同時(shí)，MSE 等于0，MSE 越大，代表預(yù)測(cè)誤差越大。

2.1 基于隨機(jī)森林對(duì)ERα 抑制劑活性的預(yù)測(cè)

2.1.1 隨機(jī)森林算法

隨機(jī)森林（Random Forest，簡(jiǎn)稱RF）是通過Bagging思想將多棵CART 回歸樹集成的一種有監(jiān)督學(xué)習(xí)算法[10]。Bagging 是根據(jù)Bootstrap 思想（有放回的隨機(jī)抽樣）構(gòu)建的一種集成學(xué)習(xí)算法[11]。CART 回歸樹最優(yōu)特征和劃分點(diǎn)的選擇依據(jù)是最小均方差，即對(duì)任意劃分特征A，其對(duì)應(yīng)的任意劃分點(diǎn)a 所劃分成的數(shù)據(jù)集和，找出使集合和的均方差最小，同時(shí)使和的均方差之和最小的劃分特征和劃分點(diǎn)，可以表達(dá)為：

其中，cleft為數(shù)據(jù)集Dleft的樣本輸出均值，cright為數(shù)據(jù)集的樣Dright本輸出均值。

本文利用隨機(jī)森林回歸模型進(jìn)行預(yù)測(cè)的步驟如下：

（1）從樣本量為N 的化合物訓(xùn)練集中有放回的隨機(jī)抽取n（n < N）個(gè)樣本，重復(fù)m 次，共生成m 個(gè)訓(xùn)練樣本集；

（2）使用訓(xùn)練樣本集構(gòu)建回歸樹，在節(jié)點(diǎn)的所有分子描述符中隨機(jī)選取部分分子描述符，依據(jù)最小均方差選擇最優(yōu)分子描述符和劃分點(diǎn)，將當(dāng)前節(jié)點(diǎn)劃分為兩個(gè)子節(jié)點(diǎn)，遞歸劃分直至滿足終止條件；

（3）重復(fù)步驟（2），構(gòu)建的m 棵回歸樹就組成了隨機(jī)森林回歸模型；

（4）輸入化合物測(cè)試樣本，m 棵樹預(yù)測(cè)值的平均值為最終預(yù)測(cè)結(jié)果，將其與真實(shí)值對(duì)比，來評(píng)價(jià)模型的預(yù)測(cè)效果。

2.1.2 隨機(jī)森林調(diào)參與結(jié)果分析

使用篩選出的50 個(gè)分子描述符作為自變量，以化合物的活性值作為因變量構(gòu)建隨機(jī)森林回歸預(yù)測(cè)模型。利用Python 的sklearn 包做隨機(jī)森林回歸預(yù)測(cè)時(shí)，主要涉及到三個(gè)重要超參數(shù)：n_estimators (回歸樹的個(gè)數(shù))、max_depth(回歸樹的最大深度)和min_samples_leaf(葉子結(jié)點(diǎn)最少樣本數(shù))。回歸樹的個(gè)數(shù)太小，模型容易欠擬合；回歸樹的個(gè)數(shù)太大會(huì)導(dǎo)致計(jì)算量過大，并且回歸樹個(gè)數(shù)增加到一定數(shù)量后，模型效果不再顯著提升?；貧w樹的最大深度過小容易導(dǎo)致模型欠擬合，過大容易導(dǎo)致模型過擬合。葉子結(jié)點(diǎn)最少樣本數(shù)涉及到回歸樹的剪枝，如果葉子結(jié)點(diǎn)數(shù)小于min_samples_leaf，則該葉子結(jié)點(diǎn)和兄弟節(jié)點(diǎn)都將被剪枝，剪枝過程可以提高隨機(jī)森林回歸模型的泛化能力。手工調(diào)制超參數(shù)需要耗費(fèi)大量時(shí)間來探索不同組合得到的效果，我們使用網(wǎng)格搜索來選擇最優(yōu)參數(shù)。分別設(shè)置n_estimators 的取值有50, 60, 70, 80,90,100，max_depth 的取值有8, 10, 12，min_samples_leaf的取值有20, 25, 30, 35, 40，同時(shí)使用5 折交叉驗(yàn)證，共有90 種n_estimators、max_depth 和min_samples_leaf的組合方式。而每一種組合方式要在訓(xùn)練集上訓(xùn)練5次，所以一共要訓(xùn)練450 次。利用網(wǎng)格搜索，進(jìn)行五折交叉驗(yàn)證訓(xùn)練隨機(jī)森林回歸模型，訓(xùn)練結(jié)束后得到的最優(yōu)超參數(shù)組合方式為n_estimators = 70、max_depth =12 和min_samples_leaf = 20。分別在訓(xùn)練集和測(cè)試集上截取40 個(gè)數(shù)據(jù)，預(yù)測(cè)效果如圖1 所示。

圖1（a）訓(xùn)練集預(yù)測(cè)效果

圖1（b）測(cè)試集預(yù)測(cè)效果

由圖1 可以看出，隨機(jī)森林回歸模型的預(yù)測(cè)效果較好，且在測(cè)試集上的預(yù)測(cè)效果與訓(xùn)練集上的預(yù)測(cè)效果相似，說明調(diào)參后的隨機(jī)森林回歸模型具有良好的穩(wěn)健性。利用網(wǎng)格搜索得到的最優(yōu)參數(shù)組合和隨機(jī)森林默認(rèn)參數(shù)分別構(gòu)建隨機(jī)森林回歸預(yù)測(cè)模型得到的均方誤差結(jié)果如表1 所示。

由結(jié)果可知，使用默認(rèn)參數(shù)構(gòu)建的隨機(jī)森林回歸預(yù)測(cè)模型，在訓(xùn)練集上的預(yù)測(cè)精度很高，但測(cè)試集均方誤差相對(duì)訓(xùn)練集均方誤差過大，產(chǎn)生了過擬合現(xiàn)象。通過網(wǎng)格搜索調(diào)整參數(shù)和使用交叉驗(yàn)證訓(xùn)練模型之后，訓(xùn)練集和測(cè)試集的預(yù)測(cè)效果都很好，均方誤差很接近，模型的泛化能力明顯提升，可以對(duì)ERα 抑制劑的活性進(jìn)行有效預(yù)測(cè)。

2.2 基于其他機(jī)器學(xué)習(xí)方法對(duì)ERα 抑制劑活性的預(yù)測(cè)

2.2.1 基于支持向量機(jī)對(duì)ERα 抑制劑活性的預(yù)測(cè)

支持向量機(jī)（Support Vector Machine,簡(jiǎn)稱SVM）由Corinna Cortes 等人于1995 年首次提出，屬于有監(jiān)督的機(jī)器學(xué)習(xí)方法，在解決非線性、小樣本和高維特征的分類和回歸問題時(shí)有很好的的效果[12]。支持向量機(jī)回歸（SVR）通過加入距離誤差epsilon 的損失函數(shù)來度量回歸精度。使用高斯函數(shù)作為支持向量機(jī)回歸模型的核函數(shù)，設(shè)置模型參數(shù)為：高斯核函數(shù)（懲罰系數(shù)C = 1.25,距離誤差epsilon = 0.1，核函數(shù)參數(shù)gamma = 0.1）,在訓(xùn)練集上和測(cè)試集上的均方誤差分別為0.653, 0.792，可知支持向量機(jī)回歸模型用于ERα 抑制劑的活性預(yù)測(cè)效果較好。

2.2.2 基于多元線性回歸對(duì)ERα 抑制劑活性的預(yù)測(cè)多元線性回歸（multiple linear regression, 簡(jiǎn)稱MLR）是QSAR 中最早采用和最經(jīng)典的數(shù)學(xué)建模方法[13]。用復(fù)相關(guān)系數(shù)R2來對(duì)多元線性回歸模型的擬合程度進(jìn)行評(píng)價(jià)。

2.3 各機(jī)器學(xué)習(xí)方法預(yù)測(cè)能力的比較

本文使用三種機(jī)器學(xué)習(xí)方法對(duì)ERα 抑制劑的活性進(jìn)行預(yù)測(cè)，對(duì)于隨機(jī)森林和支持向量機(jī)模型的建立，需要調(diào)整參數(shù)以得到更好的預(yù)測(cè)效果，對(duì)于多元線性回歸模型，需要進(jìn)行擬合優(yōu)度檢驗(yàn)來判斷模型的可用性，具體預(yù)測(cè)效果如表2 所示。三個(gè)模型均有良好的預(yù)測(cè)能力，且隨機(jī)森林方法在訓(xùn)練集和測(cè)試集上的均方誤差都比其他兩種方法的要小，表現(xiàn)出了更好的預(yù)測(cè)能力和泛化能力。

表2 三種模型預(yù)測(cè)效果比較

3 結(jié)論

本文分別使用隨機(jī)森林、支持向量機(jī)和多元線性回歸構(gòu)建了ERα 抑制劑生物活性預(yù)測(cè)模型，使用方差過濾法和Lasoo 回歸篩選出與ERα 抑制劑活性最相關(guān)的分子描述符。通過對(duì)分子描述符的合理篩選和模型參數(shù)的優(yōu)化，本文建立的ERα 抑制劑活性活性預(yù)測(cè)模型具有良好的預(yù)測(cè)效果，且隨機(jī)森林表現(xiàn)出了更好的預(yù)測(cè)能力和穩(wěn)健性，認(rèn)為隨機(jī)森林模型更適用于ERα 抑制劑的活性預(yù)測(cè)。