基于機器學習的香豆素衍生物的建模研究

張培儉, 夏潤澤, 高 湛, 劉 壯

(青島大學計算機科學技術學院, 山東 青島 266071)

乳腺癌是女性最常見的惡性腫瘤之一,也是女性癌癥死亡的主要原因,其發病率隨著年齡的增長急劇上升[1-2],全球癌癥統計發現,2018年乳腺癌新發病例約208.9萬例[3],其中包含死亡病例627 000例,約占全球女性癌癥死亡總數的15%[4]。自上世紀60年代,科學家研制出選擇性雌激素調制器(SERM),首次被批準用于輔助治療乳腺癌的SERM是他莫昔芬,該藥物已為許多患者提供治療。雖然SERM改善了ERα陽性乳腺癌患者的癥狀,但其價格昂貴,且伴有諸多副作用,嚴重限制了其治療效果[5],目前活性高但副作用低的SERM尚未發現。近幾年,香豆素發展成一種多功能的分子支架,廣泛分布在自然界,具有多種藥理及治療作用,如抗細菌,抗真菌,抗瘧疾和抗癌等[6]。研究表明,某些香豆素系的雜交體會抑制MCF-7乳腺癌細胞(ER陽性(ER +)人乳腺癌細胞系)的增殖[7],其衍生物對治療乳腺癌具有重大意義。香豆素衍生物通過對乳腺癌細胞系(MCF-7)半抑制濃度(IC50)值來判斷對癌細胞的抑制效果,傳統的化學測定方法耗費大量的時間和資源,利用定量構效關系(quantitative structure-activity relationship,QSAR)可以快速準確的預測化合物的IC50。QSAR模型是基于分子結構,通過數學方法預測未經測試的化合物的活性[8-9],GAO Z等[10]通過使用隨機森林驗證描述符可靠性,并利用混合核函數的支持向量機回歸對[1,2,3]三唑[4,5-d]嘧啶衍生物進行抗增殖預測,取得了較好效果;張克俊等人[11]利用基因表達是編程方法建立模型,有效地預測醛類化合物毒性;WANG Y等人[12]使用機器學習方法研究他克林衍生物對乙酰膽堿酯酶抑制劑的活性,為抗阿爾茲海默癥藥物研發提供了幫助;司宏宗等人[13]利用支持向量機建立預測模型,預測分子結構對藥物與血漿蛋白的結合率的影響,對藥物篩選提供參考;宋富成等人[14]通過基因表達式編程方法研究離子液體理化結構與青?；【鶴67毒性關系,取得良好的可靠性?；诖?本文基于啟發式算法篩選出的描述符,利用支持向量回歸(support vector regression,SVR),廣義回歸神經網絡(general regression neural network,GRNN)和K近鄰(K-nearest neighbor,KNN)建立了3種QSAR模型,擬合57種香豆素衍生物對MCF-7乳腺癌細胞的抑制作用并對其作出預測。實驗結果表明,利用SVR建立的模型具有較好的預測能力,且模型穩定性最強。該研究為研究治療乳腺癌的潛在藥物分子提供指導。

1 數據

文獻[4]、[13]及[15-18]均使用相同的化合物活性IC50測量方法, 即(3-(4,5-Dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide(簡稱MTT)方法測量化合物的活性值,獲取了57個香豆素衍生物的IC50值,其結構和對應的lg(IC50) 測量值和預測值如表1所示。表中化合物按照約3∶1的比例隨機劃分為訓練集和測試集,lg(IC50)值作為預測標簽,標注*符號的為測試集。

表1 香豆素衍生物及其lg(IC50)測量值和預測值

2 研究方法

2.1 分子描述符

分子描述符的計算步驟如下:

1) 利用ChemDraw軟件繪制57個香豆素衍生物的分子結構。

2) 在HyperChem中將分子結構利用MM+分子力學力場進行初步優化,通過半經驗AM1和PM3方法進一步優化,獲得穩定的分子結構[19]。

3) 將HyperChem中得到的zmt文件導入MOPAC,計算出mno文件[20]。

4) 通過CODESSA程序計算出5類描述符,即構成型、拓撲型、幾何型、靜電型和量子化學型[21-22]。

利用啟發式算法計算描述符,效率高且不受數據集大小的限制,并可自動從描述符空間中選擇特征最顯著的描述符組,其遵循的規則為:

1) 選擇化合物共有的描述符;

2) 排除所有化合物中數值變化很小的描述符;

3) 排除共線描述符,即相關系數大于0.8的描述符。

經計算,本研究最終確定使用3個描述符。

2.2 廣義回歸神經網絡

GRNN是F.S.Donaid在1991年提出的徑向基人工神經網絡的一種變體,主要用于預測和控制工廠過程建?；蛞话阌成鋯栴}[23-24]。GRNN是一種具有高度并行結構的單向學習算法,使多維測量空間中的數據稀疏,也能提供從一個觀測到另一個觀測的平穩過渡。GRNN結構如圖1所示。

圖1 GRNN結構

GRNN由輸入層、模式層、求和層和輸出層組成[25]。輸入層是接收輸入信號,神經元數目等于模型的獨立特征數量;模式層對輸入數據和訓練數據集進行必要的映射,模式層中神經元的數量通常等于訓練集中的樣本數量。模式層節點的輸出乘以適當的權重,然后在求和層將所有權重加在一起;輸出層節點負責在輸入數據集上提供所需的結果。由于預測值(香豆素衍生物IC50的lg值)的維度是1,所以求和層僅包含2個單元,且每個單元連接模式層的所有節點。第1個求和節點將模式層的所有輸出求和,并計算式(1)的分子。模式層中的第i個節點與第1個求和節點之間的連接權重等于yi,與第2個求和節點的權重等于1,輸出單元計算求和層的2個輸出的商,得出系統模型的期望輸出值。如果作為輸入向量,GRNN的輸出值為

(1)

(2)

式中,n是訓練數據的數量;x為輸入數據,y是期望輸出;σ是平滑因子,σ越大,函數近似越平滑,它的取值由具體情況確定[26]。

2.3 K近鄰

1991年,Dasarathy[27]提出KNN算法,該算法用于處理多維數據集分類或回歸監督學習,根據其最鄰近點估算未知點或缺失點的值。最近鄰通常被確定為與相鄰的未知點的距離最短的點,采用歐式距離度量2個樣本之間的距離。在高維數據中,為鄰近的變量分配不同的權重是合適的,所以本研究采用高斯函數加權估計的方法計算預測值。歐式距離函數為

(3)

式中,n為數據的特征維度。

2.4 支持向量回歸

SVR是支持向量機(support vector machine,SVM)的擴展,SVM成功應用于許多領域的回歸分析[34],它是將輸入向量映射到高維特征空間,然后在高維特征空間中進行線性回歸,利用核函數實現內積運算,降低高維空間中運算的復雜性。

在傳統的回歸問題中,給定訓練集D={(x1,y1),(x2,y2),…,(xm,ym)},其中xi∈R是輸入向量,yi∈R是對應的輸出值,預測值f(x)與真實值y的差值完全包含在預測結果的損失中。與此不同的是,SVR引入了ε不敏感損失函數,允許在f(x)和y之間最大化ε,提高模型的泛化能力,其中ε表示模型允許的預測誤差。因此,SVR的目標函數可表示為

(4)

上式中,ω為超平面系數,C為正則常數,lε為ε不敏感損失函數,其式為

(5)

(6)

在此基礎上引入拉格朗日乘子和核函數后,最終目標表達式為

(7)

2.5 評價指標

為了驗證GRNN、KNN和SVR構建回歸模型穩定性的定量指標,R2和均方根誤差ERMS的表達式為

(8)

(9)

3 結果

根據啟發式算法得到3個描述符特征MTICN(min total interaction for a C-N bond)、MBOC(max bond order of a C atom)及NOCl(number of Cl atoms),并由其建立模型,得到MTICN與MBOC對IC50產生正向影響,NOCl則產生負向影響,MBOC系數絕對值為7.852,較另外兩者大,故其對IC50具有更大的影響。其公式為

lg(IC50)=0.195d1+7.852d2-0.073d3-16.530

(10)

式中,d1、d2和d3表示3種描述符對應的數據。

3個描述符的物理-化學含義如表2所示,3個描述符的皮爾遜相關系數矩陣如表3所示。

表2 3個描述符的物理-化學含義

表3 描述符相關系數矩陣

本文隨機劃分數據并建立3個QSAR模型,采用四折交叉驗證模型的魯棒性[39]。

3.1 廣義回歸神經網絡模型

利用GRNN建立QSAR模型,采用高斯函數作為模式層的傳遞函數,由于高斯函數參數σ的取值影響GRNN的精度,所以必須選取合適的值。為了優化參數,嘗試從0.1到1均勻遞增變化,最終確定σ為0.1。GRNN訓練集和測試集的R2分別為0.949和0.911,RMSE分別為0.092和0.121,GRNN預測結果如圖2所示。

圖2 GRNN預測結果

由圖2可以看出,利用GRNN構建的模型具有良好的預測能力,其四折交叉驗證的訓練集和測試集的R2分別為0.957和0.766,RMSE分別為0.083和0.187,表現出過擬合,模型魯棒性較差。

3.2 K近鄰模型

本研究中K取值為3,距離函數為高斯函數,KNN模型預測結果如圖3所示。

由圖3可以看出,訓練集和測試集的R2分別為0.963和0.950,RMSE分別為0.077和0.088,KNN構建的模型也具有良好的預測能力,其四折交叉驗證的訓練集和測試集的R2分別為0.969和0.794,RMSE分別為0.069和0.185,在交叉驗證時,模型表現出過擬合現象,魯棒性較差。

3.3 支持向量回歸模型

SVR建立的模型包括懲罰系數C、? 和 ?-不敏感函數、核函數 κ 以及 κ 的相應參數。C為誤差的容忍度,大小取決于數據的噪聲,而噪聲通常是未知的?-不敏感函數允許存在稀疏解,訓練后?的最優值為0.1;核函數選用高斯核,形式為F(u,v)=exp(-γ|u-v|2),其中γ為常數,決定了映射到高維空間的特征向量的數量,對訓練模型的速度有顯著影響。γ越小,支持向量越多,反之亦然,最后采用網格搜索法確定最佳的C和γ:C=1.43,γ=499.768。SVR預測結果如圖4所示。

圖4 SVR預測結果

由圖4可以看出,訓練集和測試集的R2分別為0.861和0.829,RMSE分別為0.023和0.068。四折交叉驗證訓練集和測試集的R2分別為0.886和0.801,RMSE分別為0.142和0.184。因此,SVR構建的模型具有較好的預測能力和最強的魯棒性,與真實數據較吻合。

3.4 結果分析對比

通過交叉驗證對比,SVR模型預測效果良好,具有最好的魯棒性。KNN模型和GRNN模型雖然表現出很好的預測能力,但是在交叉驗證時也表現出了過擬合的趨勢,這是由于本研究數據量較小,且特征維度較低導致。

4 結束語

本文通過收集已知化合物結構性質信息建立QSAR模型,更準確地預測未知化合物結構的活性和毒性。在啟發式方法下,通過支持向量機、廣義回歸神經網絡和k近鄰方法建立3種QSAR模型,使用57種香豆素衍生物對MCF-7細胞的抑制作用進行預測,3種模型預測結果均與實際值吻合較好,且SVR模型最具魯棒性,表明本文構建的模型對未知乳腺癌藥物研發能夠提供可靠的活性預測支持。此外,本文采用的MTICN、MBOC和NOCl 3個關鍵分子描述符對香豆素衍生物活性具有重要影響,為該類衍生物的研究提供方向指引,降低藥物研發成本。然而,由于數據集數量限制及機器學習方法缺陷,無法將結論泛化到更多化合物上使用,并難以對香豆素衍生物活性變化進行合理解釋,后續將圍繞此問題進行進一步研究,為藥物篩選提供更多理論指導。