杜雪平*

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

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

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

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

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

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

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

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

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

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

2.1 基于隨機森林對ERα 抑制劑活性的預(yù)測

2.1.1 隨機森林算法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.2.1 基于支持向量機對ERα 抑制劑活性的預(yù)測

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

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

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

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

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

3 結(jié)論

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