基于機器學習方法的ERα 抑制劑活性預測

杜雪平*

（湖北工業大學理學院，湖北 武漢 430068）

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

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

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

1 數據集劃分與特征篩選

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

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

2 基于機器學習方法對ERα 抑制劑活性的預測

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

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

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

2.1.1 隨機森林算法

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

其中，cleft為數據集Dleft的樣本輸出均值，cright為數據集的樣Dright本輸出均值。

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

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

（2）使用訓練樣本集構建回歸樹，在節點的所有分子描述符中隨機選取部分分子描述符，依據最小均方差選擇最優分子描述符和劃分點，將當前節點劃分為兩個子節點，遞歸劃分直至滿足終止條件；

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

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

2.1.2 隨機森林調參與結果分析

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

圖1（a）訓練集預測效果

圖1（b）測試集預測效果

由圖1 可以看出，隨機森林回歸模型的預測效果較好，且在測試集上的預測效果與訓練集上的預測效果相似，說明調參后的隨機森林回歸模型具有良好的穩健性。利用網格搜索得到的最優參數組合和隨機森林默認參數分別構建隨機森林回歸預測模型得到的均方誤差結果如表1 所示。

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

2.2 基于其他機器學習方法對ERα 抑制劑活性的預測

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

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

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

2.3 各機器學習方法預測能力的比較

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

表2 三種模型預測效果比較

3 結論

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