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

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

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

1 數(shù)據(jù)

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

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

2 研究方法

2.1 分子描述符

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

1) 利用ChemDraw軟件繪制57個香豆素衍生物的分子結(jié)構(gòu)。

2) 在HyperChem中將分子結(jié)構(gòu)利用MM+分子力學(xué)力場進行初步優(yōu)化,通過半經(jīng)驗AM1和PM3方法進一步優(yōu)化,獲得穩(wěn)定的分子結(jié)構(gòu)[19]。

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

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

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

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

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

3) 排除共線描述符,即相關(guān)系數(shù)大于0.8的描述符。

經(jīng)計算,本研究最終確定使用3個描述符。

2.2 廣義回歸神經(jīng)網(wǎng)絡(luò)

GRNN是F.S.Donaid在1991年提出的徑向基人工神經(jīng)網(wǎng)絡(luò)的一種變體,主要用于預(yù)測和控制工廠過程建模或一般映射問題[23-24]。GRNN是一種具有高度并行結(jié)構(gòu)的單向?qū)W習(xí)算法,使多維測量空間中的數(shù)據(jù)稀疏,也能提供從一個觀測到另一個觀測的平穩(wěn)過渡。GRNN結(jié)構(gòu)如圖1所示。

圖1 GRNN結(jié)構(gòu)

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

(1)

(2)

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

2.3 K近鄰

1991年,Dasarathy[27]提出KNN算法,該算法用于處理多維數(shù)據(jù)集分類或回歸監(jiān)督學(xué)習(xí),根據(jù)其最鄰近點估算未知點或缺失點的值。最近鄰?fù)ǔ１淮_定為與相鄰的未知點的距離最短的點,采用歐式距離度量2個樣本之間的距離。在高維數(shù)據(jù)中,為鄰近的變量分配不同的權(quán)重是合適的,所以本研究采用高斯函數(shù)加權(quán)估計的方法計算預(yù)測值。歐式距離函數(shù)為

(3)

式中,n為數(shù)據(jù)的特征維度。

2.4 支持向量回歸

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

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

(4)

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

(5)

(6)

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

(7)

2.5 評價指標

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

(8)

(9)

3 結(jié)果

根據(jù)啟發(fā)式算法得到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產(chǎn)生正向影響,NOCl則產(chǎn)生負向影響,MBOC系數(shù)絕對值為7.852,較另外兩者大,故其對IC50具有更大的影響。其公式為

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

(10)

式中,d1、d2和d3表示3種描述符對應(yīng)的數(shù)據(jù)。

3個描述符的物理-化學(xué)含義如表2所示,3個描述符的皮爾遜相關(guān)系數(shù)矩陣如表3所示。

表2 3個描述符的物理-化學(xué)含義

表3 描述符相關(guān)系數(shù)矩陣

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

3.1 廣義回歸神經(jīng)網(wǎng)絡(luò)模型

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

圖2 GRNN預(yù)測結(jié)果

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

3.2 K近鄰模型

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

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

3.3 支持向量回歸模型

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

圖4 SVR預(yù)測結(jié)果

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

3.4 結(jié)果分析對比

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

4 結(jié)束語

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