張培儉, 夏潤(rùn)澤, 高 湛, 劉 壯

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

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

1 數(shù)據(jù)

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

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

2 研究方法

2.1 分子描述符

分子描述符的計(jì)算步驟如下:

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

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

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

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

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

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

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

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

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

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

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

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

GRNN由輸入層、模式層、求和層和輸出層組成[25]。輸入層是接收輸入信號(hào),神經(jīng)元數(shù)目等于模型的獨(dú)立特征數(shù)量;模式層對(duì)輸入數(shù)據(jù)和訓(xùn)練數(shù)據(jù)集進(jìn)行必要的映射,模式層中神經(jīng)元的數(shù)量通常等于訓(xùn)練集中的樣本數(shù)量。模式層節(jié)點(diǎn)的輸出乘以適當(dāng)?shù)臋?quán)重,然后在求和層將所有權(quán)重加在一起;輸出層節(jié)點(diǎn)負(fù)責(zé)在輸入數(shù)據(jù)集上提供所需的結(jié)果。由于預(yù)測(cè)值(香豆素衍生物IC50的lg值)的維度是1,所以求和層僅包含2個(gè)單元,且每個(gè)單元連接模式層的所有節(jié)點(diǎn)。第1個(gè)求和節(jié)點(diǎn)將模式層的所有輸出求和,并計(jì)算式(1)的分子。模式層中的第i個(gè)節(jié)點(diǎn)與第1個(gè)求和節(jié)點(diǎn)之間的連接權(quán)重等于yi,與第2個(gè)求和節(jié)點(diǎn)的權(quán)重等于1,輸出單元計(jì)算求和層的2個(gè)輸出的商,得出系統(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ù)其最鄰近點(diǎn)估算未知點(diǎn)或缺失點(diǎn)的值。最近鄰?fù)ǔ１淮_定為與相鄰的未知點(diǎn)的距離最短的點(diǎn),采用歐式距離度量2個(gè)樣本之間的距離。在高維數(shù)據(jù)中,為鄰近的變量分配不同的權(quán)重是合適的,所以本研究采用高斯函數(shù)加權(quán)估計(jì)的方法計(jì)算預(yù)測(cè)值。歐式距離函數(shù)為

(3)

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

2.4 支持向量回歸

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

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

(4)

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

(5)

(6)

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

(7)

2.5 評(píng)價(jià)指標(biāo)

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

(8)

(9)

3 結(jié)果

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

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

(10)

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

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

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

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

本文隨機(jī)劃分?jǐn)?shù)據(jù)并建立3個(gè)QSAR模型,采用四折交叉驗(yàn)證模型的魯棒性[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)練集和測(cè)試集的R2分別為0.949和0.911,RMSE分別為0.092和0.121,GRNN預(yù)測(cè)結(jié)果如圖2所示。

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

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

3.2 K近鄰模型

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

由圖3可以看出,訓(xùn)練集和測(cè)試集的R2分別為0.963和0.950,RMSE分別為0.077和0.088,KNN構(gòu)建的模型也具有良好的預(yù)測(cè)能力,其四折交叉驗(yàn)證的訓(xùn)練集和測(cè)試集的R2分別為0.969和0.794,RMSE分別為0.069和0.185,在交叉驗(yàn)證時(shí),模型表現(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ù)量,對(duì)訓(xùn)練模型的速度有顯著影響。γ越小,支持向量越多,反之亦然,最后采用網(wǎng)格搜索法確定最佳的C和γ:C=1.43,γ=499.768。SVR預(yù)測(cè)結(jié)果如圖4所示。

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

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

3.4 結(jié)果分析對(duì)比

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

4 結(jié)束語

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