叢 湧 薛 英,2,*

(1四川大學(xué)化學(xué)學(xué)院,教育部綠色化學(xué)與技術(shù)重點(diǎn)實(shí)驗(yàn)室,成都610064;2西華大學(xué)四川省先進(jìn)科學(xué)計(jì)算重點(diǎn)實(shí)驗(yàn)室,成都610039)

1 引言

丙型肝炎病毒(HCV)是通過血液傳染的非甲、非乙型肝炎(non-A,non-B viral hepatitis)的主要致病因子.全世界約有1.7億慢性丙肝病毒感染者,每年有35萬余人死于與丙肝相關(guān)的肝臟疾病.1,2由于缺少特異有效的藥物和疫苗,目前臨床上采用α-2-聚乙二醇化干擾素(pegylated interferon-α-2a)與利巴韋林(ribavirin)聯(lián)合用藥治療丙型肝炎,其治療有效率僅為50%,治療周期較長(zhǎng),不良反應(yīng)發(fā)生頻率高且比較嚴(yán)重.3HCV的高度變異性使治療和預(yù)防面臨巨大困難,急需研發(fā)有效的抗HCV藥物,以補(bǔ)充完善目前現(xiàn)有的治療方法.近年來,尋找HCV特定靶向抗病毒治療藥物(STAT-C)是抗HCV研究的重要方向,這些特定靶向的靶點(diǎn)包括NS3/NsS4A蛋白酶/解旋酶和非結(jié)構(gòu)蛋白5B(NS5B)RNA依賴的聚合酶(NS5B RNA-dependent RNA polymerase)等,其中以丙型肝炎病毒NS5B RNA聚合酶為靶標(biāo)的抗HCV藥物研究頗受關(guān)注.4

NS5B是一種病毒編碼的RNA依賴性RNA聚合酶,可以調(diào)控丙型肝炎病毒RNA模板(-)鏈的合成及(+)鏈基因組RNA的再生,5在丙型肝炎病毒復(fù)制進(jìn)程中必不可少.研究發(fā)現(xiàn)通過對(duì)NS5B的抑制可以有效阻斷丙肝病毒復(fù)制,以NS5B為靶點(diǎn)的抑制劑設(shè)計(jì)、合成及生物測(cè)試越來越成為該領(lǐng)域的研究熱點(diǎn),一系列具有不同骨架結(jié)構(gòu)的核苷型(NIs)、非核苷型(NNIs)NS5B抑制劑被不斷的合成及報(bào)道,6其中包括苯并咪唑,7吲哚,8噻吩,9吖啶酮衍生物,10苯并噻二嗪11等類型抑制劑.然而,這些研究主要集中在合成新的抑制劑化合物及測(cè)試其生物抑制活性,實(shí)驗(yàn)周期較長(zhǎng)且資金消耗巨大.因此,在早期藥物研發(fā)過程中,為了縮短藥物開發(fā)周期和節(jié)約成本,我們急需大力發(fā)展計(jì)算機(jī)模擬技術(shù)對(duì)小分子抑制劑生物活性做出初步定性和定量的預(yù)測(cè),并從建立的數(shù)據(jù)統(tǒng)計(jì)模型中挖掘大分子蛋白靶點(diǎn)與小分子配體的可能作用機(jī)理,進(jìn)一步指導(dǎo)和輔助新型抗HCV藥物分子的設(shè)計(jì)和發(fā)現(xiàn).

結(jié)構(gòu)-活性關(guān)系研究(SAR)已經(jīng)成功應(yīng)用于抗HCV抑制劑的發(fā)現(xiàn).在這些方法當(dāng)中,定量結(jié)構(gòu)-活性關(guān)系(QSAR)研究分子結(jié)構(gòu)與其所表達(dá)生物活性之間的相關(guān)性,在廣泛實(shí)踐中已被證明是一種非常有效的計(jì)算機(jī)輔助藥物設(shè)計(jì)工具.與其他方法相比,QSAR方法的優(yōu)點(diǎn)在于其建立的數(shù)據(jù)統(tǒng)計(jì)模型可以直觀指示在生物抑制活性確定中發(fā)揮重要作用的小分子抑制劑結(jié)構(gòu)性因素;構(gòu)建模型所需的分子描述符獨(dú)立于任何實(shí)驗(yàn)條件,可以通過小分子抑制劑三維結(jié)構(gòu)直接計(jì)算得出;通過QSAR方法得到的構(gòu)效關(guān)系能夠提供非常有用的信息,這些信息可以進(jìn)一步指導(dǎo)及輔助基于配體(ligand-based)和基于受體(receptor-based)的藥物設(shè)計(jì).然而,遺憾的是,基于機(jī)器學(xué)習(xí)的二維定量構(gòu)效關(guān)系研究并不能明確揭示生物大分子蛋白受體與小分子抑制劑配體結(jié)合的三維構(gòu)象;在這樣的二維定量構(gòu)效關(guān)系研究中,研究人員經(jīng)常不會(huì)考慮樣本分子集三維構(gòu)象的疊合和取向.基于大分子靶蛋白和小分子抑制劑三維構(gòu)象的3D-QSAR仍然是定量構(gòu)效關(guān)系研究的發(fā)展趨勢(shì),也是我們工作組未來的研究重點(diǎn).傳統(tǒng)的化學(xué)信息學(xué)方法多采用多元線性回歸(MLR)、啟發(fā)式方法(HM)、主成分回歸(PCR)和偏最小二乘方法(PLS)等線性定量構(gòu)效關(guān)系方法建模.近年來,支持向量機(jī)(SVM)和徑向基函數(shù)網(wǎng)絡(luò)(RBF)等非線性回歸方法在多樣性分子結(jié)構(gòu)樣本集中對(duì)藥效學(xué)、藥代動(dòng)力學(xué)和毒理學(xué)性質(zhì)不斷顯示其優(yōu)異的預(yù)測(cè)性能,12,13正越來越多的受到關(guān)注.Melagraki等14對(duì)98個(gè)苯并噻二嗪衍生物抗HCV抑制劑進(jìn)行了定量構(gòu)效研究,通過采用線性逐步回歸特征消除選擇方法(ES-SWR)從ChemSar和Topix軟件計(jì)算的分子描述符集中篩選出5個(gè)重要的分子描述符建立QSAR模型.他們的QSAR模型對(duì)訓(xùn)練集的平方相關(guān)系數(shù)(R2)僅為0.74,對(duì)驗(yàn)證集的R2為0.81.Li研究小組15運(yùn)用最佳多元線性回歸方法(BMLR)選擇建模最優(yōu)描述符子集,并采用多元線性回歸、徑向基函數(shù)神經(jīng)網(wǎng)絡(luò)(RBFNN)和支持向量機(jī)方法建立線性和非線性QSAR模型對(duì)118個(gè)苯并噻二嗪衍生物抗HCV抑制劑進(jìn)行了定量構(gòu)效關(guān)系研究,15其中非線性RBFNN和SVM模型給出更為準(zhǔn)確和理想的預(yù)測(cè)結(jié)果,兩種非線性模型對(duì)于訓(xùn)練集的R2分別為0.850和0.875,對(duì)于測(cè)試集的R2分別為0.893和0.854.雖然苯并噻二嗪衍生物抗HCV抑制劑的定量構(gòu)效關(guān)系研究已經(jīng)取得了一些有成效的結(jié)果,然而如何選擇合適的特征選擇和建模機(jī)器學(xué)習(xí)方法仍然是影響QSAR模型預(yù)測(cè)能力的兩個(gè)關(guān)鍵因素.本工作研究了89個(gè)最新報(bào)道的苯并異噻唑(benzoisothiazole)和苯并噻嗪(benzothiazine)類抗HCV抑制劑的定量構(gòu)效關(guān)系.我們嘗試使用遺傳算法組合偏最小二乘(GA-PLS)和線性逐步回歸分析(LSRA)特征選擇方法選擇最優(yōu)描述符子集.對(duì)每種特征選擇方法選擇的描述符子集分別采用多元線性回歸、偏最小二乘、遺傳算法組合支持向量機(jī)(GA-SVM)三種方法用訓(xùn)練集建立QSAR模型,并用這些模型預(yù)測(cè)了測(cè)試集中化合物的NS5B抑制活性.

2 材料與方法

2.1 丙型肝炎病毒聚合酶NS5B抑制劑數(shù)據(jù)集的選取

本工作從最近發(fā)表的文獻(xiàn)中總共收集了89個(gè)苯并異噻唑16,17和苯并噻嗪18-20類丙型肝炎病毒聚合酶NS5B抑制劑.這些化合物的IC50值(1-868000 nmol·L-1)均用相同的實(shí)驗(yàn)測(cè)量方法在相同的實(shí)驗(yàn)條件下測(cè)定,被測(cè)抑制劑中大部分是高效NS5B抑制劑.首先使用ChemDraw軟件繪制每個(gè)小分子抑制劑的二維結(jié)構(gòu),隨后通過CORINA軟件將小分子二維平面結(jié)構(gòu)轉(zhuǎn)換成三維結(jié)構(gòu),再采用量子化學(xué)AM1計(jì)算方法優(yōu)化每一個(gè)抑制劑分子的三維構(gòu)型,然后對(duì)計(jì)算結(jié)果進(jìn)行手動(dòng)檢查以確保每個(gè)優(yōu)化分子都生成正確的手性結(jié)構(gòu)且沒有重復(fù).為了便于研究,我們將抑制劑的IC50值轉(zhuǎn)換成pIC50值(lg(109/IC50))作為QSAR模型的因變量.根據(jù)化合物的結(jié)構(gòu)與化學(xué)性質(zhì)在化學(xué)空間中的相似性和分布,21將所有的化合物分為訓(xùn)練集和測(cè)試集兩大類.訓(xùn)練集(含45個(gè)小分子抑制劑)訓(xùn)練和優(yōu)化回歸模型,測(cè)試集(含44個(gè)小分子抑制劑)評(píng)價(jià)回歸模型的預(yù)測(cè)能力.

2.2 分子描述符的計(jì)算

分子描述符在定量構(gòu)效關(guān)系研究中經(jīng)常被用于定量描述分子的結(jié)構(gòu)和物理化學(xué)特性.在我們的研究中,通過手動(dòng)方法從相關(guān)文獻(xiàn)中找出一千多個(gè)分子描述符,并根據(jù)化合物拓?fù)浣Y(jié)構(gòu)、電子結(jié)構(gòu)和幾何結(jié)構(gòu)等性質(zhì),剔除明顯冗余和與預(yù)測(cè)藥物性質(zhì)不相關(guān)的描述符,22最終篩選出189個(gè)與化合物性質(zhì)密切相關(guān)的分子描述符(見表1).其中包括18個(gè)簡(jiǎn)單分子性質(zhì)描述符(如分子量、可旋轉(zhuǎn)的鍵數(shù)),27個(gè)分子連接性和形狀描述符(如分子連接性指數(shù)和分子卡伯形狀指數(shù)),97個(gè)電拓?fù)鋺B(tài)分子描述符(如電子拓?fù)鋺B(tài)指數(shù)),22個(gè)量子化學(xué)性質(zhì)分子描述符(如原子電荷和分子的偶極矩)和25個(gè)分子幾何特性描述符(如溶劑可及表面積和疏水區(qū)域).我們應(yīng)用實(shí)驗(yàn)室自編的分子描述符計(jì)算程序,根據(jù)AM1方法優(yōu)化過后的化合物三維結(jié)構(gòu)計(jì)算所有的分子描述符.為了減少建模中描述符之間高度自相關(guān)帶來的多重共線性干擾,在特征選擇之前對(duì)這189個(gè)分子描述符集做了預(yù)處理,步驟如下:(1)移除在90%的樣本分子中具有相同數(shù)值的描述符;(2)移除相對(duì)標(biāo)準(zhǔn)偏差小于0.05的描述符;(3)對(duì)于Pearson相關(guān)系數(shù)超過0.95的一對(duì)描述符,留下與生物活性相關(guān)性較高的描述符,剔除另外一個(gè).23經(jīng)過上述預(yù)處理,我們最終保留了85個(gè)分子描述符進(jìn)行下一步的特征選擇.

2.3 特征選擇方法

2.3.1 逐步回歸分析法

我們借助SPSS軟件自帶的逐步回歸分析程序選擇最優(yōu)描述符子集,逐步回歸分析法中每步有兩個(gè)過程即引進(jìn)變量和剔除變量,且引進(jìn)變量和剔除變量均需作F檢驗(yàn)后方可繼續(xù)進(jìn)行,故又稱為雙重檢驗(yàn)回歸分析法.其具體步驟如下:(1)引入變量,引入變量的原則是未引進(jìn)變量中偏回歸平方和最大者并經(jīng)過F顯著性檢驗(yàn),若顯著則引進(jìn),否則終止.(2)剔除變量,剔除原則是在引進(jìn)的自變量中偏回歸平方和最小者,并經(jīng)過F檢驗(yàn)不顯著,則剔除.(3)終止條件即最優(yōu)條件,再無顯著自變量引進(jìn),也沒有不顯著自變量可以剔除.

2.3.2 遺傳算法組合偏最小二乘方法

我們借助MATLAB偏最小二乘-遺傳算法工具箱24實(shí)現(xiàn)最優(yōu)描述符子集的選取.GA-PLS是一種基于遺傳算法的優(yōu)化工具,25,26其算法過程描述如下:(1)定義和編碼染色體;(2)種群的初始化;(3)評(píng)價(jià)每個(gè)染色體的適應(yīng)度;(4)保護(hù)染色體;(5)保留最好的染色體;(6)對(duì)種群進(jìn)行交叉和變異遺傳操作;(7)若滿足終止條件停止程序,否則轉(zhuǎn)入步驟3.本實(shí)驗(yàn)采用了GA-PLS工具箱里的三個(gè)函數(shù),它們分別是GAPLSOPT(dataset,1),GAPLSOPT(dataset,2)和GAPLS(dataset,the number of evaluation,precision).

GAPLSOPT(dataset,1)函數(shù)測(cè)試樣本數(shù)據(jù)集是否適用GA-PLS方法選取特征,根據(jù)GA-PLS軟件設(shè)計(jì)者介紹,如果GAPLSOPT函數(shù)對(duì)于樣本數(shù)據(jù)集的測(cè)試輸出結(jié)果在0到5之間,使用GA-PLS方法對(duì)數(shù)據(jù)集進(jìn)行特征選擇是安全穩(wěn)定的.圖1中顯示了本工作89個(gè)小分子抑制劑對(duì)應(yīng)的85個(gè)分子描述符樣本數(shù)據(jù)集的GAPLSOPT(dataset,1)輸出結(jié)果,樣本數(shù)據(jù)集的隨機(jī)測(cè)試結(jié)果在0到4.9468之間,這說明采用GA-PLS方法對(duì)該數(shù)據(jù)集進(jìn)行特征選擇是可靠的.為了避免GA-PLS方法在訓(xùn)練過程中產(chǎn)生過擬合,我們使用GAPLSOPT(dataset,2)函數(shù)估計(jì)GAPLS函數(shù)所需的最優(yōu)評(píng)價(jià)次數(shù)(the number of evaluation)參數(shù),如圖2所示,GAPLSOPT(dataset,2)的差異曲線在評(píng)價(jià)次數(shù)為115處有全局最大值,因此,評(píng)價(jià)次數(shù)被置為115作為GAPLS函數(shù)的最優(yōu)控制參數(shù).經(jīng)過上述準(zhǔn)備工作,我們運(yùn)行GAPLS函數(shù)對(duì)數(shù)據(jù)集進(jìn)行特征選擇.為了減少隨機(jī)誤差,我們重復(fù)了10次GAPLS實(shí)驗(yàn)得到平均結(jié)果.圖3顯示了交叉驗(yàn)證響應(yīng)和每個(gè)描述符的被選擇頻率,在85個(gè)分子描述符中,GAPLS函數(shù)最終選出7個(gè)分子描述符用于下一步的QSAR建模.

表1 所用的分子描述符Table 1 Molecular descriptors used in this work

2.4 遺傳算法組合支持向量機(jī)方法

圖2 GAPLSOPT(2)差異曲線Fig.2 GAPLSOPT(2)difference curve

支持向量機(jī)是基于結(jié)構(gòu)風(fēng)險(xiǎn)最小化(SRM)理論的機(jī)器學(xué)習(xí)方法,其兩分類理論經(jīng)過多年發(fā)展已日趨成熟.通過引入ε-不敏感損失函數(shù),支持向量機(jī)也可以擴(kuò)展用來解決小樣本數(shù)據(jù)集的回歸問題.支持向量機(jī)經(jīng)由核函數(shù)K(xi,x)將輸入空間的X向量映射到高維希爾伯特空間H,其常用的核函數(shù)包括線性、多項(xiàng)式、徑向基和S形等函數(shù)形式,其中徑向基函數(shù)由于良好的非線性映像能力已在諸多領(lǐng)域得到了廣泛應(yīng)用.27本研究使用高斯徑向基函數(shù)(RBF)構(gòu)建支持向量機(jī)回歸模型.高斯核函數(shù)具體表示如下:

我們?cè)贙(xi,x)特征空間構(gòu)造決策函數(shù)(方程2)作為ε-支持向量機(jī)的最優(yōu)解:

其中b為方程的偏置項(xiàng).

圖3 GAPLS函數(shù)描述符選擇頻率估計(jì)Fig.3 Selected frequency figure by GAPLS function

支持向量回歸機(jī)的泛化性能取決于最優(yōu)正則化參數(shù)C,不敏感參數(shù)ε和RBF核函數(shù)寬度σ的選取.28正則化參數(shù)C對(duì)回歸函數(shù)的復(fù)雜性和泛化能力進(jìn)行折衷.在確定的數(shù)據(jù)子空間中,參數(shù)C取得太小,則對(duì)樣本數(shù)據(jù)中超出ε不敏感帶的樣本懲罰就越小,使訓(xùn)練誤差變大,系統(tǒng)的泛化能力變差,會(huì)出現(xiàn)“欠學(xué)習(xí)”現(xiàn)象;C取得太大,相應(yīng)的權(quán)重就小,系統(tǒng)的泛化能力變差,會(huì)出現(xiàn)“過學(xué)習(xí)”現(xiàn)象.不敏感參數(shù)ε控制著ε不敏感帶的寬度,影響著支持向量的數(shù)目.ε值選得太小,回歸估計(jì)精度高,但支持向量數(shù)目增多,ε選的太大,回歸估計(jì)精度降低,支持向量數(shù)目減少,支持向量機(jī)的稀疏性大.RBF核函數(shù)寬度σ反映了訓(xùn)練樣本數(shù)據(jù)的分布或范圍特性,它確定了局部領(lǐng)域的寬度,較大的σ意味著較低的方差.

遺傳算法是模擬達(dá)爾文生物進(jìn)化論自然選擇和遺傳學(xué)機(jī)理的生物進(jìn)化過程的計(jì)算模型,是一種通過模擬自然進(jìn)化過程搜索最優(yōu)化解的方法.29本工作中我們采用遺傳算法同時(shí)對(duì)支持向量回歸機(jī)的三個(gè)參數(shù)(正則化參數(shù)C,不敏感參數(shù)ε和RBF核函數(shù)寬度參數(shù)σ)進(jìn)行優(yōu)化.新型進(jìn)化計(jì)算框體—Python環(huán)境下分布式進(jìn)化算法(DEAP)30被用來實(shí)現(xiàn)SVM參數(shù)優(yōu)化的遺傳算法框架,Libsvm程序31被用于構(gòu)建ε-SVM回歸模型,通過引入Python多路處理模塊(The multi-processing techniques in Python),我們實(shí)現(xiàn)了GA-SVM算法的并發(fā),大大提高了算法的計(jì)算速度.在遺傳算法中,染色體使用二進(jìn)制編碼,每個(gè)染色體由C:(Ci,i=1-nc),ε:(εj,j=1-nε)和σ:(σk,k=1-nσ)三部分組成,Ci表示染色體中表征正則化參數(shù)C二進(jìn)制位串中第i位的數(shù)值,σk表示染色體中表征RBF核函數(shù)寬度參數(shù)σ二進(jìn)制位串中第j位的數(shù)值,εk表示染色體中表征不敏感參數(shù)ε二進(jìn)制位串中第k位的數(shù)值,nc、nσ和nε分別表示染色體中表征C、σ和ε的二進(jìn)制位串的長(zhǎng)度(在本工作中nc=10,nσ=10,nε=10).根據(jù)解碼公式(5),這三個(gè)二進(jìn)制編碼位串被分別轉(zhuǎn)換成正則化參數(shù)C、不敏感參數(shù)ε和RBF核函數(shù)寬度參數(shù)σ的實(shí)數(shù)值.

其中R表示二進(jìn)制位串所表征相應(yīng)參數(shù)的真實(shí)值,maxR和minR為相應(yīng)參數(shù)最大和最小指數(shù)冪取值范圍(參數(shù)取值在2minR和2maxR

之間),d為二進(jìn)制位串轉(zhuǎn)十進(jìn)制數(shù)值,l為二進(jìn)制位串的長(zhǎng)度.我們?cè)O(shè)置C的取值范圍為2-10-215;ε的取值范圍為:2-10

-210;σ的取值范圍為2-10-28;n重交叉驗(yàn)證均方誤差(MSECV)作為遺傳算法的適應(yīng)度函數(shù),定義如下:

其中,yi是實(shí)驗(yàn)真實(shí)值,?i是模型的計(jì)算值,ntest是n重交叉驗(yàn)證集的樣本化合物數(shù)目.遺傳算法的選擇、交叉和變異操作,我們通過調(diào)用DEAP內(nèi)置函數(shù)——selTournament(individuals,k,tournsize)、cxTwo Points(ind1,ind2)和mutFlipBit(individual,indpb)來實(shí)現(xiàn),設(shè)置種群規(guī)模為500,總共迭代100代,交叉率和變異率分別設(shè)置為0.6和0.2.當(dāng)某代種群中95%以上的個(gè)體得到相同的均方誤差適應(yīng)度數(shù)值停止算法,作為GA-SVM的停機(jī)準(zhǔn)則.

3 結(jié)果和討論

3.1 兩種特征選擇方法選出的最優(yōu)描述符子集及討論

逐步回歸分析法共選出6個(gè)分子描述符,其中包括2個(gè)簡(jiǎn)單分子性質(zhì)描述符、1個(gè)分子連接性和形狀描述符、2個(gè)量子化學(xué)性質(zhì)描述符和1個(gè)分子幾何特性描述符;遺傳算法組合偏最小二乘方法共選出7個(gè)分子描述符,其中包括3個(gè)簡(jiǎn)單分子性質(zhì)描述符、1個(gè)分子連接性和形狀描述符、1個(gè)電拓?fù)鋺B(tài)分子描述符、2個(gè)量子化學(xué)性質(zhì)描述符和1個(gè)分子幾何特性描述符.如表2所示,兩種特征選擇方法同時(shí)選出4個(gè)相同的分子描述符,分別為分子中雜原子數(shù)目(nhet)、氧原子數(shù)目(noxy)、分子最大負(fù)電荷(AQ,min)和親水性指數(shù)(Hiwpl).

起初苯并噻嗪類抑制劑的構(gòu)性關(guān)系研究集中在苯并噻嗪環(huán)的C-7位,18通過在該位置引入取代基來改善抑制活性.研究發(fā)現(xiàn)在苯并噻嗪環(huán)的C-7位引入極性基團(tuán)如羥基或甲基磺酰胺基會(huì)提高抑制劑對(duì)酶NS5B的抑制效力.RNA聚合酶NS5B和其抑制劑的復(fù)合物晶體結(jié)構(gòu)進(jìn)一步證實(shí)在1,3-二羰基功能基團(tuán)上的負(fù)電荷易在環(huán)上發(fā)生離域與聚合酶NS5B骨架上Tyr448殘基的N-H基團(tuán)及兩個(gè)保守水分子反生靜電相互作用;與此同時(shí),苯并噻嗪環(huán)中的苯環(huán)正好與Phe193殘基面對(duì)面接觸,苯環(huán)之間發(fā)生π-π堆積作用.在C-7位引入甲基磺酰胺基的苯并噻嗪類抑制劑中,砜氧基團(tuán)與橋接Ser556殘基和甲基磺酰胺基團(tuán)的結(jié)構(gòu)水分子形成氫鍵靜電相互作用;磺酰胺基團(tuán)也會(huì)與NS5B的Asn291殘基形成靜電相互作用;氨磺?；鶊F(tuán)中酸性的N-H與NS5B Asp318殘基產(chǎn)生更強(qiáng)有力的氫鍵相互作用力,磺酰胺基團(tuán)的這種獨(dú)特的三點(diǎn)式相互作用有助于解釋抑制劑中該極性基團(tuán)的引入對(duì)于抑制活性的改善.構(gòu)性關(guān)系還發(fā)現(xiàn)在苯并噻嗪環(huán)的C-2和C-5位摻入吸電子基團(tuán)降低環(huán)電子密度,可以有效改善膽汁的轉(zhuǎn)運(yùn)識(shí)別,并同時(shí)減少代謝負(fù)產(chǎn)物的生成.20

在所選描述符中,nhet(分子中所含N,O和S等雜原子數(shù)目)、noxy(分子中O原子數(shù)目)和nsulph(分子中S原子數(shù)目)這三個(gè)簡(jiǎn)單性質(zhì)描述符間接表征了抑制劑分子中極性基團(tuán)(如羥基或甲基磺酰胺基等)與NS5B發(fā)生靜電相互作用及形成氫鍵的能力;親水性指數(shù)(Hiwpl)描述抑制劑分子與NS5B酶的親水疏水相互作用;兩個(gè)量子化學(xué)描述符,分子最大負(fù)電荷(AQ,min)和分子平均負(fù)電荷(Mnc)可能與苯并噻嗪環(huán)的電子離域化及C-2,C-5位上吸電基團(tuán)的引入有關(guān);S(1)(羥基H原子電拓?fù)鋺B(tài)指數(shù))和QH,Max(分子中氫原子上的最大正電荷)給出氫鍵給體的質(zhì)子提供能力.從上面的討論可以看出,兩種方法所選擇的描述符能夠反映抑制劑分子的電荷分布、氫鍵相互作用、分子連接性、親水疏水相互作用等性質(zhì).

表2 LSRA和GA-PLS特征選擇方法選取的分子描述符Table 2 Molecular descriptors selected from the LSRAand GA-PLS feature selection methods

3.2 逐步回歸分析法所選描述符的建?；貧w結(jié)果

我們采用逐步回歸分析法所選的6個(gè)最優(yōu)描述符分別建立多元線性、偏最小二乘和支持向量機(jī)回歸模型,訓(xùn)練集(含45個(gè)化合物)用來訓(xùn)練和優(yōu)化模型,測(cè)試集(含44個(gè)小分子抑制劑)評(píng)價(jià)回歸模型的預(yù)測(cè)能力.訓(xùn)練得到的多元線性回歸方程為:pIC50=0.329×nhet+0.383×noxy-0.16723.903×Mnc+0.051×Hiwpl-5.588,0.931,s2=0.144,F=99.473,為調(diào)整相關(guān)系數(shù)平方,s2為回歸分析的標(biāo)準(zhǔn)偏差平方,即回歸分析方差;F為回歸分析的顯著性檢驗(yàn).圖S1(Supporting Information)為多元線性回歸模型對(duì)于訓(xùn)練集和測(cè)試集的實(shí)驗(yàn)與預(yù)測(cè)pIC50數(shù)值對(duì)比圖,模型對(duì)于訓(xùn)練集、測(cè)試集及整個(gè)數(shù)據(jù)集的均方誤差(MSE)分別為0.121、0.122和0.122,相應(yīng)的相關(guān)系數(shù)分別為0.970、0.958和0.965.偏最小二乘方法抽選出3個(gè)主成分,訓(xùn)練得到的回歸方程為:pIC50=0.2153×nhet+0.3652×noxy-0.1463×Mnc+0.1519×Hiwpl-7.3150s2=0.166,F=85.454,N=45,圖S2為偏最小二乘回歸模型對(duì)于訓(xùn)練集和測(cè)試集的實(shí)驗(yàn)與預(yù)測(cè)pIC50數(shù)值對(duì)比圖,模型對(duì)于訓(xùn)練集、測(cè)試集及整個(gè)數(shù)據(jù)集的MSE分別為0.140、0.122和0.131,相應(yīng)的相關(guān)系數(shù)分別為0.965、0.958和0.961.由于影響NS5B酶抑制劑活性的分子特征極其復(fù)雜,并不是所有的分子描述符都與活性存在絕對(duì)的線性關(guān)系,為了與上述線性回歸建模方法做比較,基于相同的最優(yōu)描述符子集,我們使用遺傳算法組合支持向量機(jī)方法建立非線性QSAR模型.首先,我們?cè)谟?xùn)練集上使用GASVM方法同時(shí)優(yōu)化SVM模型的三個(gè)參數(shù),當(dāng)SVM模型參數(shù)取值為:C=7.7387,σ=0.3546,ε=0.04664時(shí),回歸模型給出最好的留一法交叉驗(yàn)證誤差(MSECV=0.135304);然后,我們采用這組優(yōu)化參數(shù)建立SVM模型.圖S3為SVM回歸模型對(duì)于訓(xùn)練集和測(cè)試集的實(shí)驗(yàn)與預(yù)測(cè)pIC50數(shù)值對(duì)比圖,模型對(duì)于訓(xùn)練集、測(cè)試集及整個(gè)數(shù)據(jù)集的MSE分別為0.113、0.108和0.111,相應(yīng)的相關(guān)系數(shù)分別為0.972、0.962和0.968.在上述三種方法建立的QSAR模型中,支持向量機(jī)給出最好的回歸模型,而多元線性回歸給出預(yù)測(cè)效果最好的線性回歸模型,三種機(jī)器學(xué)習(xí)模型對(duì)抑制劑活性的預(yù)測(cè)結(jié)果詳見表S1.

3.3 遺傳算法組合偏最小二乘法所選描述符的建?；貧w結(jié)果

我們采用GA-PLS所選的7個(gè)最優(yōu)描述符子集分別建立偏最小二乘、多元線性和支持向量機(jī)回歸模型.訓(xùn)練得到的多元線性回歸方程為:pIC50=0.207×nhet+0.295×nsulph+0.304×noxy-0.035×S(1)+F=76.566,N=45,圖S4為多元線性回歸模型對(duì)于訓(xùn)練集和測(cè)試集的實(shí)驗(yàn)與預(yù)測(cè)pIC50數(shù)值對(duì)比圖,模型對(duì)于訓(xùn)練集、測(cè)試集及整個(gè)數(shù)據(jù)集的MSE分別為0.131、0.243和0.186,相應(yīng)的相關(guān)系數(shù)分別為0.967、0.918和0.946.偏最小二乘方法抽選出3個(gè)主成分,訓(xùn)練得到的回歸方程為:pIC50=0.1667×nhet+0.5267×nsulph+0.2915×noxy+0.0220×S(1)+8.1236×QH,Max-8.7093×AQ,min+0.0827×Hiwpl-3.7535,R2=0.933,R2adjust=0.920,s2=0.163,F=73.606,N=45,圖S5為偏最小二乘回歸模型對(duì)于訓(xùn)練集和測(cè)試集的實(shí)驗(yàn)與預(yù)測(cè)pIC50數(shù)值對(duì)比圖.模型對(duì)于訓(xùn)練集、測(cè)試集及整個(gè)數(shù)據(jù)集的MSE分別為0.134、0.114和0.124,相應(yīng)的相關(guān)系數(shù)分別為0.966、0.960和0.964.基于相同的最優(yōu)描述符子集,使用遺傳算法組合支持向量機(jī)方法建立非線性QSAR模型,當(dāng)SVM模型參數(shù)取值為:C=11.620579,σ=0.182558,ε=0.224924時(shí),回歸模型給出最好的留一法交叉驗(yàn)證誤差(MSECV=0.145978),圖S6為SVM回歸模型對(duì)于訓(xùn)練集和測(cè)試集的實(shí)驗(yàn)與預(yù)測(cè)pIC50數(shù)值對(duì)比圖,模型對(duì)于訓(xùn)練集、測(cè)試集及整個(gè)數(shù)據(jù)集的MSE分別為0.123、0.118和0.121,相應(yīng)的相關(guān)系數(shù)分別為0.970、0.958和0.965.在上述三種方法建立的QSAR模型中,支持向量機(jī)回歸模型對(duì)于訓(xùn)練集和整個(gè)數(shù)據(jù)集給出最好的預(yù)測(cè)結(jié)果;兩種線性回歸模型中,偏最小二乘模型對(duì)測(cè)試集和整個(gè)數(shù)據(jù)集給出最好的預(yù)測(cè)結(jié)果,而多元線性回歸模型對(duì)于訓(xùn)練集給出最好的預(yù)測(cè)效果,三種機(jī)器學(xué)習(xí)模型對(duì)抑制劑活性的預(yù)測(cè)結(jié)果詳見表S2.

4 結(jié)論

本工作研究了89個(gè)苯并異噻唑和苯并噻嗪類抗HCV抑制劑的定量構(gòu)效關(guān)系.線性逐步回歸分析和遺傳算法組合偏最小二乘方法被用來選取最優(yōu)描述符子集,對(duì)于如上兩種特征選擇方法所選描述符,我們分別采用多元線性回歸、偏最小二乘、遺傳算法組合支持向量機(jī)方法建模并得到了比較滿意的預(yù)測(cè)結(jié)果.在采用LSRA所選描述符建立的三個(gè)QSAR模型中,支持向量機(jī)給出最好的回歸模型,多元線性回歸給出預(yù)測(cè)效果最好的線性模型;在采用GA-PLS所選描述符建立的三個(gè)QSAR模型中,支持向量機(jī)同樣給出最好的回歸模型,而偏最小二乘回歸給出效果最好的線性模型.研究結(jié)果表明,非線性支持向量機(jī)方法的建模回歸效果好于傳統(tǒng)的多元線性回歸和偏最小二乘方法;采用LSRA和GA-PLS特征選擇方法所選描述符建立的模型都得到比較滿意的回歸效果,這表明兩種特征選擇方法都能篩選合適的描述符子集用于QSAR建模;多元線性回歸和偏最小二乘方法建立的線性QSAR方程可以直觀指示在生物抑制活性確定中發(fā)揮重要作用的小分子結(jié)構(gòu)和物理化學(xué)特征信息,這些信息可以進(jìn)一步指導(dǎo)及輔助基于配體和受體的抗HCV藥物設(shè)計(jì).

Supporting Information: The plots of experimental vs computational lg(109/IC50)values are given in Figs.S1-S6.The information of the investigated dataset is provided in Tables S1 and S2.This information is available free of charge via the internet at http://www.whxb.pku.edu.cn.

(1) Choo,Q.L.;Weiner,A.J.;Overby,L.R.;Bradley,D.W.;Houghton,M.Science 1989,244,359.doi:10.1126/science.2523562

(2)(a)Lauer,G.M.;Walker,B.D.N.Engl.J.Med.2001,345,41.doi:10.1056/NEJM200107053450107(b)Di Bisceglie,A.M.Lancet 1998,351,351.(c)Alter,M.J.;Kruszon-Moran,D.;Nainan,O.V.;McQuillan,G.M.;Gao,F.;Moyer,L.A.;Kaslow,R.A.;Margolis,H.S.N.Engl.J.Med.1999,341,556.

(3)Manns,M.P.;McHutchison,J.G.;Gordon,S.C.;Rustgi,V.K.;Shiffman,M.;Reindollar,R.;Goodman,Z.D.;Koury,K.;Ling,M.H.;Albrecht,J.K.Lancet 2002,347,975.

(4) (a)Koch,U.;Narjes,F.Curr.Top.Med.Chem.2007,7,1302.doi:10.2174/156802607781212211(b)R?nn,R.;Sandstr?m,A.Curr.Top.Med.Chem.2008,8,533.(c)Zapf,C.W.;Bloom,J.D.;Levin,J.I.Ann.Rep.Med.Chem.2007,42,281.

(5) Appel,N.;Schaller,T.;Penin,F.;Bartenschlager,R.J.Biol.Chem.2006,281,9833.doi:10.1074/jbc.R500026200

(6) Ni,Z.J.;Wagman,A.S.Curr.Opin.Drug Discov.Dev.2004,7,446.

(7) Beaulieu,P.L.;Bos,M.;Bousquet,Y.;Fazal,G.;Gauthier,J.;Gillard,J.;Goulet,S.;LaPlante,S.;Poupart,M.A.;Lefebvre,S.;McKercher,G.;Pellerin,C.;Austel,V.;Kukolj,G.Bioorg.Med.Chem.Lett.2004,14,119.doi:10.1016/j.bmcl.2003.10.023

(8) Stansfield,I.;Ercolani,C.;Mackay,A.;Conte,I.;Pompei,M.;Koch,U.;Gennari,N.;Giuliano,C.;Rowley,M.;Narjes,F.Bioorg.Med.Chem.Lett.2009,19,627.doi:10.1016/j.bmcl.2008.12.068

(9)Louise-May,S.;Yang,W.;Nie,X.;Liu,D.;Deshpande,M.S.;Phadke,A.S.;Huang,M.;Agarwal,A.Bioorg.Med.Chem.Lett.2007,17,3905.doi:10.1016/j.bmcl.2007.04.103

(10) Stankiewicz-Drogon,A.;Palchykovska,L.G.;Kostina,V.G.;Alexeeva,I.V.;Shved,A.D.;Boguszewska-Chachulska,A.M.Bioorg.Med.Chem.2008,16,8846.doi:10.1016/j.bmc.2008.08.074

(11) Bosse,T.D.;Larson,D.P.;Wagner,R.;Hutchinson,D.K.;Rockway,T.W.;Kati,W.M.;Liu,Y.;Masse,S.;Middleton,T.;Mo,H.;Montgomery,D.;Jiang,W.;Koev,G.;Kempf,D.J.;Molla,A.Bioorg.Med.Chem.Lett.2008,18,568.doi:10.1016/j.bmcl.2007.11.088

(12) Lü,W.J.;Chen,Y.L.;Ma,W.P.;Zhang,X.Y.;Luan,F.;Liu,M.C.;Chen,X.G.;Hu,Z.D.Euro.J.Med.Chem.2008,43,569.doi:10.1016/j.ejmech.2007.04.011

(13)Luan,F.;Liu,H.T.;Ma,W.P.;Fan,B.T.Euro.J.Med.Chem.2008,43,43.doi:10.1016/j.ejmech.2007.03.002

(14) Melagraki,G.;Afantitis,A.;Sarimveis,H.;Koutentis,P.A.;Markopoulos,J.;Igglessi-Markopoulou,O.Bioorg.Med.Chem.2007,15,7237.doi:10.1016/j.bmc.2007.08.036

(15) Su,L.;Li,L.;Li,Y.;Zhang,X.;Huang,X.;Zhai,H.Med.Chem.Res.2012,21,2079.doi:10.1007/s00044-011-9734-x

(16) deVicente,J.;Hendricks,R.T.;Smith,D.B.;Fell,J.B.;Fischer,J.;Spencer,S.R.;Stengel,P.J.;Mohr,P.;Robinson,J.E.;Blake,J.F.;Hilgenkamp,R.K.;Yee,C.;Adjabeng,G.;Elworthy,T.R.;Li,J.;Wang,B.;Bamberg,J.T.;Harris,S.F.;Wong,A.;Leveque,V.J.P.;Najera,I.;Pogam,S.L.;Rajyaguru,S.;Ao-Ieong,G.;Alexandrova,L.;Larrabee,S.;Brandl,M.;Briggs,A.;Sukhtankar,S.;Farrell,R.Bioorg.Med.Chem.Lett.2009,19,5652.doi:10.1016/j.bmcl.2009.08.022

(17) Hendricks,R.T.;Spencer,S.R.;Blake,J.F.;Fell,J.B.;Fischer,J.;Stengel,P.J.;Leveque,V.J.P.;Pogam,S.L.;Rajyaguru,S.;Najera,I.;Swallow,S.Bioorg.Med.Chem.Lett.2009,19,410.doi:10.1016/j.bmcl.2008.11.060

(18) deVicente,J.;Hendricks,R.T.;Smith,D.B.;Fell,J.B.;Fischer,J.;Spencer,S.R.;Stengel,P.J.;Mohr,P.;Robinson,J.E.;Blake,J.F.;Hilgenkamp,R.K.;Yee,C.;Adjabeng,G.;Elworthy,T.R.;Tracy,J.;Chin,E.;Li,J.;Wang,B.;Bamberg,J.T.;Stephenson,R.;Oshiro,C.;Harris,S.F.;Ghate,M.;Leveque,V.;Najera,I.;Pogam,S.L.;Rajyaguru,S.;Ao-Ieong,G.;Alexandrova,L.;Larrabee,S.;Brandl,M.;Briggs,A.;Sukhtankar,S.;Farrell,R.;Xu,B.Bioorg.Med.Chem.Lett.2009,19,3642.doi:10.1016/j.bmcl.2009.05.004

(19) Hendricks,R.T.;Fell,J.B.;Blake,J.F.;Fischer,J.P.;Robinson,J.E.;Spencer,S.R.;Stengel,P.J.;Bernacki,A.L.;Leveque,V.J.P.;Pogam,S.L.;Rajyaguru,S.;Najera,I.;Josey,J.A.;Harris,J.R.;Swallow,S.Bioorg.Med.Chem.Lett.2009,19,3637.doi:10.1016/j.bmcl.2009.04.119

(20) deVicente,J.;Hendricks,R.T.;Smith,D.B.;Fell,J.B.;Fischer,J.;Spencer,S.R.;Stengel,P.J.;Mohr,P.;Robinson,J.E.;Blake,J.F.;Hilgenkamp,R.K.;Yee,C.;Zhao,J.;Elworthy,T.R.;Tracy,J.;Chin,E.;Li,J.;Lui,A.;Wang,B.;Oshiro,C.;Harris,S.F.;Ghate,M.;Leveque,V.J.P.;Najera,I.;Pogam,S.L.;Rajyaguru,S.;Ao-Ieong,G.;Alexandrova,L.;Fitch,B.;Brandl,M.;Masjedizadeh,M.;Wua,S.Y.;de Keczer,S.;Voronin,T.Bioorg.Med.Chem.Lett.2009,19,5648.doi:10.1016/j.bmcl.2009.08.023

(21) Todeschini,R.;Consonni,V.Handbook of Molecular Descriptors;Wiley-VCH:New York,2000.

(22)Xue,Y.;Li,Z.R.;Yap,C.W.;Sun,L.Z.;Chen,X.;Chen,Y.Z.J.Chem.Inform.Comp.Sci.2004,44,1630.doi:10.1021/ci049869h

(23)Tan,N.X.;Rao,H.B.;Li,Z.R.;Li,X.Y.SAR QSAR Environ.Res.2009,20,27.doi:10.1080/10629360902724085

(24) http://www.models.kvl.dk/source/GAPLS/index.asp,accessed June 2008.

(25) Leardi,R.;Boggia,R.;Terrile,M.J.Chemom.1992,6,267.

(26) Leardi,R.J.Chemom.1994,8,65.

(27) Burbidge,R.;Trotter,M.;Buxton,B.;Holden,S.Comput.Chem.2001,26,5.doi:10.1016/S0097-8485(01)00094-8

(28) Cherkassky,V.;Ma,Y.Selection of Meta-parameters for Support Vector Regression.Proceedings of the International Conference onArtificial Neural Networks,Madrid,Spain,Aug 28-30,2002.

(29)Hao,M.;Li,Y.;Wang,Y.;Zhang,S.Anal.Chim.Acta 2011,690,53.doi:10.1016/j.aca.2011.02.004

(30) Rainville,F.M.D.;Fortin,F.A.;Gardner,M.A.;Parizeau,M.;Gagné,C.DEAP:APython Framework for Evolutionary Algorithms.In EvoSoft Workshop,Companion Proc.of the GeneticandEvolutionaryComputationConference,July 07-11,2012.

(31)Chang,C.C.;Lin,C.J.LIBSVM:ALibrary for Support Vector Machines,2001.Software available at http://www.csie.ntu.edu.tw/-cjlin/libsvm,accessed Jun 2008.