彭奕萍

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

引言

在藥物研發(fā)中，一個(gè)化合物要想成為候選藥物，必須要在人體內(nèi)具備良好的藥代動(dòng)力學(xué)性質(zhì)和安全性，合稱為ADMET（Absorption 吸收、Distribution 分布、Metabolism 代謝、Excretion 排泄、Toxicity 毒性）性質(zhì)[1]。ADME 主要代表了化合物的藥代動(dòng)力學(xué)性質(zhì)，描述了化合物在生物體內(nèi)的濃度隨時(shí)間變化的規(guī)律。一個(gè)化合物在擁有良好生物活性的同時(shí)，還需要有良好的ADMET 性質(zhì)。如果化合物很難被人體吸收、在人體內(nèi)代謝太快、或者具有某種毒性，即便擁有良好的活性也無法使其成為藥物。因此在藥物研發(fā)早期，對化合物成藥性和安全性進(jìn)行一定的評估，可以幫助提高藥物研發(fā)成功率、降低研發(fā)成本。

本文先使用低方差過濾和帶L1 正則的邏輯回歸進(jìn)行特征選擇，接著使用邏輯斯蒂、線性支持向量機(jī)、非線性支持向量機(jī)、隨機(jī)森林和梯度提升樹模型五類常用機(jī)器學(xué)習(xí)算法建模，對每個(gè)化合物的ADMET 性質(zhì)進(jìn)行分類預(yù)測。最后選擇AUC 作為評估指標(biāo)，對各個(gè)模型的分類預(yù)測結(jié)果進(jìn)行對比評估。

1 數(shù)據(jù)

本文在藥物ADMET 的重要性質(zhì)中，選擇小腸上皮細(xì)胞滲透性（Caco-2）、代謝穩(wěn)定性（CYP3A4）、心臟毒性（hERG）、生物利用度（HOB）和遺傳毒性（MN），這5 個(gè)藥物研發(fā)中關(guān)鍵的性質(zhì)進(jìn)行建模[2]。

數(shù)據(jù)來源于“華為杯”第十八屆中國研究生數(shù)學(xué)建模競賽D 題。數(shù)據(jù)集一給出了1 974 個(gè)化合物的729 個(gè)分子描述符信息，數(shù)據(jù)集的大小為1 974×729?；衔锏姆肿用枋龇且幌盗杏糜诿枋龌衔锝Y(jié)構(gòu)和性質(zhì)特征的參數(shù)，包括物理化學(xué)性質(zhì)和拓?fù)浣Y(jié)構(gòu)特征等。數(shù)據(jù)集二記錄了化合物結(jié)構(gòu)的SMILES式，和每個(gè)化合物的5 種ADMET 性質(zhì)，其中，ADMET性質(zhì)采用二分類記錄相應(yīng)取值。Caco-2 取值為“1”代表該化合物的小腸上皮細(xì)胞滲透性較好，“0”代表較差；CYP3A4 取值“1”代表該化合物能夠被CYP3A4代謝，“0”代表不能被代謝；hERG 取值為“1”代表該化合物具有心臟毒性，“0”代表不具有心臟毒性；HOB值為“1”代表該化合物的口服生物利用度較好，“0”代表該化合物口服的生物利用度差；MN 值為“1”代表該化合物具有遺傳毒性，為“0”代表該化合物不具有遺傳毒性?；衔顰DMET 數(shù)據(jù)集描述如表1 所示。

表1 ADMET 數(shù)據(jù)集描述

2 模型的建立

本文使用了Logistic 回歸[3]、線性支持向量機(jī)、非線性支持向量機(jī)[4]、隨機(jī)森林[5]、梯度提升樹[6]分別構(gòu)建了化合物ADMET 性質(zhì)的分類預(yù)測模型。實(shí)驗(yàn)建模流程如下頁圖1 所示。

圖1 建模流程圖

2.1 特征選擇

特征選擇是從大量特征中選出一組特征，刪除冗余或者不相關(guān)的信息，提高分類器的性能。本問題原始數(shù)據(jù)集中共有729 個(gè)分子描述符，數(shù)據(jù)維度過大會(huì)在一定程度上影響模型的收斂。為了方便后續(xù)的建模，需要對特征進(jìn)行選擇。

特征選擇算法，按照集合評價(jià)策略可分為，基于濾波（Filter）的特征選擇算法和基于嵌入式（Wrapper）的特征選擇算法[7]。濾波特征選擇的優(yōu)勢在于可以很快排除很大數(shù)量的非關(guān)鍵建模特征，縮小特征搜索范圍，但它并不能保證選擇出一個(gè)規(guī)模較小的優(yōu)特征子集，同時(shí)由于它與后續(xù)的機(jī)器學(xué)習(xí)模型無關(guān)，可能導(dǎo)致較差的模型性能，因此常用來作為特征選擇的預(yù)篩選器。嵌入式方法將特征選擇過程嵌入到機(jī)器學(xué)習(xí)模型中，即利用機(jī)器學(xué)習(xí)來為每一個(gè)特征打分，根據(jù)特征得分來篩選特征。與過濾法相比，速度上雖然要比濾波法慢很多，但其所選的特征子集規(guī)模往往較小，能更好地篩選出建模關(guān)鍵特征。

基于以上分析，本文在特征選擇上先使用了低方差過濾法作一輪特征預(yù)篩選，篩選出方差為零的特征。特征的方差為零，說明該特征在所有觀測點(diǎn)的值都相同，意味著該特征無法解釋目標(biāo)變量（ADMET 性質(zhì)）的任何變化，對建模沒有作用，故將其刪除。接著使用嵌入法再進(jìn)行一輪特征選擇（本文中使用的是帶L1 正則的邏輯回歸），進(jìn)一步剔除對化合物ADMET性質(zhì)影響較小的噪聲特征。在特征選擇環(huán)節(jié)，通過結(jié)合濾波法和嵌入法，在保證特征篩選的效率的同時(shí)，也兼顧了關(guān)鍵特征篩選的質(zhì)量。本文特征選擇過程如圖2 所示。

圖2 特征選擇流程圖

本文利用Python 語言編程，讀入數(shù)據(jù)并計(jì)算729個(gè)特征的方差，通過將方差的閾值設(shè)為零，剔除掉方差為零的特征（無法解釋化合物ADMET 性質(zhì)的特征）。經(jīng)過低方差過濾后，快速排除了225 個(gè)非關(guān)鍵建模特征，得到的容量為504 的特征集合，但是這個(gè)特征維度仍然較大，不利于模型的收斂，進(jìn)而影響模型的性能。為了進(jìn)一步降低數(shù)據(jù)維度尋找較優(yōu)的建模特征集合，在經(jīng)過第一輪的方差過濾法過濾掉一部分特征后，繼續(xù)使用帶L1 正則的邏輯回歸進(jìn)行特征選擇。通過在邏輯回歸的損失函數(shù)中引入L1 正則項(xiàng)防止過擬合，同時(shí)線性模型中的L1 正則項(xiàng)能夠有效地將某些特征的特征系數(shù)縮小為零，從而實(shí)現(xiàn)解的稀疏。而基于帶正則項(xiàng)線性模型的特征系數(shù)，可以實(shí)現(xiàn)為特征打分，系數(shù)越高往往也就意味著該特征在模型中越重要。利用Python 語言編程實(shí)現(xiàn)上述特征選擇算法，為了有利于邏輯回歸模型的收斂，先將數(shù)據(jù)作歸一化處理，使用scikit-learn 建立帶L1 正則的邏輯回歸模型，帶入歸一化后的數(shù)據(jù)訓(xùn)練模型，并在驗(yàn)證集上計(jì)算模型表現(xiàn)，根據(jù)模型的特征系數(shù)，篩選出對ADMET性質(zhì)影響最大的50 個(gè)分子描述符。

2.2 模型訓(xùn)練

本文使用特征選擇得到的化合物的50 個(gè)分子描述符作為建模特征，分別以化合物的5 個(gè)ADMET 性質(zhì)作為目標(biāo)變量，利用邏輯斯蒂回歸、線性支持向量機(jī)、非線性支持向量機(jī)、隨機(jī)森林、梯度提升樹算法構(gòu)建化合物ADMET 性質(zhì)的分類預(yù)測模型。

利用Python 語言，將數(shù)據(jù)集中的1 974 條數(shù)據(jù)，按照3∶1∶1 的比例劃分為訓(xùn)練集、測試集和驗(yàn)證集。編程實(shí)現(xiàn)模型的構(gòu)建和訓(xùn)練，并且在模型訓(xùn)練中使用網(wǎng)格搜索技術(shù)，通過交叉驗(yàn)證來評估超參數(shù)值的所有可能組合，尋找最優(yōu)的超參數(shù)組合。

以構(gòu)建Caco-2（小腸上皮細(xì)胞滲透性）這一化合物ADEMT 性質(zhì)的邏輯回歸分類預(yù)測模型為例，使用Python 的scikit-learn 包的GridSearchCV 函數(shù)實(shí)現(xiàn)網(wǎng)格搜索，尋找最優(yōu)超參數(shù)。對于本問題，邏輯回歸模型使用到的幾個(gè)超參數(shù)為：C（正則化系數(shù)λ 的倒數(shù)）、tol（停止迭代求解的閾值）、solver（優(yōu)化算法選擇參數(shù)）、max_iter（算法收斂最大迭代次數(shù)）。實(shí)驗(yàn)中設(shè)置C可選取值為 [0.4，0.6，0.8，1]；tol的取值為[1e-4，1e-3，1e-2]；solver 的取值為['newton-cg'，'sag'，'lbfgs']；max_iter 的取值為[5，10，15，20]。利用網(wǎng)格搜索，對所有的超參數(shù)組合（共114 種）分別訓(xùn)練邏輯回歸模型，同時(shí)通過3 折交叉驗(yàn)證進(jìn)行評估，最終得到最優(yōu)的超參數(shù)組合為C=0.4、tol=0.0001、solver='newton-cg'、max_iter=10。

本文針對化合物的5 個(gè)ADMET 性質(zhì)，使用5 種機(jī)器學(xué)習(xí)方法，共構(gòu)建并訓(xùn)練了25 個(gè)分類預(yù)測模型，將訓(xùn)練好的模型保存為.pkl 文件。

3 分類預(yù)測與結(jié)果評估

3.1 評估指標(biāo)的選擇

在模型評估過程中，二值分類問題常采用的評估指標(biāo)有準(zhǔn)確率和AUC。準(zhǔn)確率是指分類正確的樣本占總樣本的比例，是最簡單也是最直觀的評價(jià)指標(biāo)。AUC 指的是ROC 曲線下的面積大小，能夠量化反映出基于ROC 曲線衡量出的模型性能。

當(dāng)正負(fù)樣本不均衡時(shí)，占比大的樣本往往會(huì)成為影響準(zhǔn)確率的主要因素，因此不適合選擇準(zhǔn)確率作為評價(jià)指標(biāo)。而由于ROC 曲線的橫坐標(biāo)為假陽性率，縱坐標(biāo)為真陽性率。當(dāng)正負(fù)樣本分布發(fā)生變化時(shí)，ROC曲線的形狀能夠基本保持不變，這個(gè)特點(diǎn)讓ROC 曲線能夠盡量降低正負(fù)樣本不均衡帶來的影響，從而更加穩(wěn)定地反映出模型的好壞。

本問題中，通過對數(shù)據(jù)集的分析統(tǒng)計(jì)可以發(fā)現(xiàn)，數(shù)據(jù)存在一定程度上的正負(fù)樣本失衡的現(xiàn)象，所以本文最終選擇AUC 來作為分類預(yù)測結(jié)果的評估指標(biāo)。

3.2 預(yù)測結(jié)果評估

加載pkl 模型文件，分別計(jì)算25 個(gè)分類預(yù)測模型在驗(yàn)證集上的預(yù)測結(jié)果得到預(yù)測值，并編程計(jì)算分類預(yù)測模型在驗(yàn)證集上的AUC 值，如表2 所示。

表2 分類預(yù)測結(jié)果對比

以AUC 值作為評價(jià)指標(biāo)，在該數(shù)據(jù)集上隨機(jī)森林算法相較于其他算法來說具有較高的適用性，在5 個(gè)ADMET 性質(zhì)的預(yù)測中都取得了很好的預(yù)測效果。針對Caco-2 的分類預(yù)測中，邏輯回歸模型表現(xiàn)最好，取得了0.970 的曲線下面積?；诜蔷€性支持向量機(jī)的hERG 分類預(yù)測模型，相較于其他模型，預(yù)測效果更好，達(dá)到了0.976 的AUC 值。梯度提升樹模型在MN 性質(zhì)的預(yù)測中表現(xiàn)最優(yōu)，AUC 值達(dá)到了0.976。隨機(jī)森林算法在CPY3A4 和HOB 的分類預(yù)測中都取得了最優(yōu)的模型效果，AUC 值分別為0.973 和0.975。模型復(fù)雜度最小的線性支持向量機(jī)在所有模型中表現(xiàn)最差，但是在針對hERG 和HOB 性質(zhì)的分類預(yù)測中能夠逼近最優(yōu)效果，結(jié)合訓(xùn)練時(shí)間成本和模型泛化能力考慮，仍然有它的價(jià)值。