基于機器學(xué)習(xí)方法的丙型肝炎病毒聚合酶NS5B 非核苷抑制劑的定量構(gòu)效關(guān)系研究

叢 湧 薛 英,2,*

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

1 引言

丙型肝炎病毒(HCV)是通過血液傳染的非甲、非乙型肝炎(non-A,non-B viral hepatitis)的主要致病因子.全世界約有1.7億慢性丙肝病毒感染者,每年有35萬余人死于與丙肝相關(guān)的肝臟疾病.1,2由于缺少特異有效的藥物和疫苗,目前臨床上采用α-2-聚乙二醇化干擾素(pegylated interferon-α-2a)與利巴韋林(ribavirin)聯(lián)合用藥治療丙型肝炎,其治療有效率僅為50%,治療周期較長,不良反應(yīng)發(fā)生頻率高且比較嚴(yán)重.3HCV的高度變異性使治療和預(yù)防面臨巨大困難,急需研發(fā)有效的抗HCV藥物,以補充完善目前現(xiàn)有的治療方法.近年來,尋找HCV特定靶向抗病毒治療藥物(STAT-C)是抗HCV研究的重要方向,這些特定靶向的靶點包括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ù)制進程中必不可少.研究發(fā)現(xiàn)通過對NS5B的抑制可以有效阻斷丙肝病毒復(fù)制,以NS5B為靶點的抑制劑設(shè)計、合成及生物測試越來越成為該領(lǐng)域的研究熱點,一系列具有不同骨架結(jié)構(gòu)的核苷型(NIs)、非核苷型(NNIs)NS5B抑制劑被不斷的合成及報道,6其中包括苯并咪唑,7吲哚,8噻吩,9吖啶酮衍生物,10苯并噻二嗪11等類型抑制劑.然而,這些研究主要集中在合成新的抑制劑化合物及測試其生物抑制活性,實驗周期較長且資金消耗巨大.因此,在早期藥物研發(fā)過程中,為了縮短藥物開發(fā)周期和節(jié)約成本,我們急需大力發(fā)展計算機模擬技術(shù)對小分子抑制劑生物活性做出初步定性和定量的預(yù)測,并從建立的數(shù)據(jù)統(tǒng)計模型中挖掘大分子蛋白靶點與小分子配體的可能作用機理,進一步指導(dǎo)和輔助新型抗HCV藥物分子的設(shè)計和發(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)與其所表達生物活性之間的相關(guān)性,在廣泛實踐中已被證明是一種非常有效的計算機輔助藥物設(shè)計工具.與其他方法相比,QSAR方法的優(yōu)點在于其建立的數(shù)據(jù)統(tǒng)計模型可以直觀指示在生物抑制活性確定中發(fā)揮重要作用的小分子抑制劑結(jié)構(gòu)性因素;構(gòu)建模型所需的分子描述符獨立于任何實驗條件,可以通過小分子抑制劑三維結(jié)構(gòu)直接計算得出;通過QSAR方法得到的構(gòu)效關(guān)系能夠提供非常有用的信息,這些信息可以進一步指導(dǎo)及輔助基于配體(ligand-based)和基于受體(receptor-based)的藥物設(shè)計.然而,遺憾的是,基于機器學(xué)習(xí)的二維定量構(gòu)效關(guān)系研究并不能明確揭示生物大分子蛋白受體與小分子抑制劑配體結(jié)合的三維構(gòu)象;在這樣的二維定量構(gòu)效關(guān)系研究中,研究人員經(jīng)常不會考慮樣本分子集三維構(gòu)象的疊合和取向.基于大分子靶蛋白和小分子抑制劑三維構(gòu)象的3D-QSAR仍然是定量構(gòu)效關(guān)系研究的發(fā)展趨勢,也是我們工作組未來的研究重點.傳統(tǒng)的化學(xué)信息學(xué)方法多采用多元線性回歸(MLR)、啟發(fā)式方法(HM)、主成分回歸(PCR)和偏最小二乘方法(PLS)等線性定量構(gòu)效關(guān)系方法建模.近年來,支持向量機(SVM)和徑向基函數(shù)網(wǎng)絡(luò)(RBF)等非線性回歸方法在多樣性分子結(jié)構(gòu)樣本集中對藥效學(xué)、藥代動力學(xué)和毒理學(xué)性質(zhì)不斷顯示其優(yōu)異的預(yù)測性能,12,13正越來越多的受到關(guān)注.Melagraki等14對98個苯并噻二嗪衍生物抗HCV抑制劑進行了定量構(gòu)效研究,通過采用線性逐步回歸特征消除選擇方法(ES-SWR)從ChemSar和Topix軟件計算的分子描述符集中篩選出5個重要的分子描述符建立QSAR模型.他們的QSAR模型對訓(xùn)練集的平方相關(guān)系數(shù)(R2)僅為0.74,對驗證集的R2為0.81.Li研究小組15運用最佳多元線性回歸方法(BMLR)選擇建模最優(yōu)描述符子集,并采用多元線性回歸、徑向基函數(shù)神經(jīng)網(wǎng)絡(luò)(RBFNN)和支持向量機方法建立線性和非線性QSAR模型對118個苯并噻二嗪衍生物抗HCV抑制劑進行了定量構(gòu)效關(guān)系研究,15其中非線性RBFNN和SVM模型給出更為準(zhǔn)確和理想的預(yù)測結(jié)果,兩種非線性模型對于訓(xùn)練集的R2分別為0.850和0.875,對于測試集的R2分別為0.893和0.854.雖然苯并噻二嗪衍生物抗HCV抑制劑的定量構(gòu)效關(guān)系研究已經(jīng)取得了一些有成效的結(jié)果,然而如何選擇合適的特征選擇和建模機器學(xué)習(xí)方法仍然是影響QSAR模型預(yù)測能力的兩個關(guān)鍵因素.本工作研究了89個最新報道的苯并異噻唑(benzoisothiazole)和苯并噻嗪(benzothiazine)類抗HCV抑制劑的定量構(gòu)效關(guān)系.我們嘗試使用遺傳算法組合偏最小二乘(GA-PLS)和線性逐步回歸分析(LSRA)特征選擇方法選擇最優(yōu)描述符子集.對每種特征選擇方法選擇的描述符子集分別采用多元線性回歸、偏最小二乘、遺傳算法組合支持向量機(GA-SVM)三種方法用訓(xùn)練集建立QSAR模型,并用這些模型預(yù)測了測試集中化合物的NS5B抑制活性.

2 材料與方法

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

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

2.2 分子描述符的計算

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

2.3 特征選擇方法

2.3.1 逐步回歸分析法

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

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

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

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

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

2.4 遺傳算法組合支持向量機方法

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

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

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

其中b為方程的偏置項.

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

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

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

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

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

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

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

3 結(jié)果和討論

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

逐步回歸分析法共選出6個分子描述符,其中包括2個簡單分子性質(zhì)描述符、1個分子連接性和形狀描述符、2個量子化學(xué)性質(zhì)描述符和1個分子幾何特性描述符;遺傳算法組合偏最小二乘方法共選出7個分子描述符,其中包括3個簡單分子性質(zhì)描述符、1個分子連接性和形狀描述符、1個電拓?fù)鋺B(tài)分子描述符、2個量子化學(xué)性質(zhì)描述符和1個分子幾何特性描述符.如表2所示,兩種特征選擇方法同時選出4個相同的分子描述符,分別為分子中雜原子數(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位引入極性基團如羥基或甲基磺酰胺基會提高抑制劑對酶NS5B的抑制效力.RNA聚合酶NS5B和其抑制劑的復(fù)合物晶體結(jié)構(gòu)進一步證實在1,3-二羰基功能基團上的負(fù)電荷易在環(huán)上發(fā)生離域與聚合酶NS5B骨架上Tyr448殘基的N-H基團及兩個保守水分子反生靜電相互作用;與此同時,苯并噻嗪環(huán)中的苯環(huán)正好與Phe193殘基面對面接觸,苯環(huán)之間發(fā)生π-π堆積作用.在C-7位引入甲基磺酰胺基的苯并噻嗪類抑制劑中,砜氧基團與橋接Ser556殘基和甲基磺酰胺基團的結(jié)構(gòu)水分子形成氫鍵靜電相互作用;磺酰胺基團也會與NS5B的Asn291殘基形成靜電相互作用;氨磺?；鶊F中酸性的N-H與NS5B Asp318殘基產(chǎn)生更強有力的氫鍵相互作用力,磺酰胺基團的這種獨特的三點式相互作用有助于解釋抑制劑中該極性基團的引入對于抑制活性的改善.構(gòu)性關(guān)系還發(fā)現(xiàn)在苯并噻嗪環(huán)的C-2和C-5位摻入吸電子基團降低環(huán)電子密度,可以有效改善膽汁的轉(zhuǎn)運識別,并同時減少代謝負(fù)產(chǎn)物的生成.20

在所選描述符中,nhet(分子中所含N,O和S等雜原子數(shù)目)、noxy(分子中O原子數(shù)目)和nsulph(分子中S原子數(shù)目)這三個簡單性質(zhì)描述符間接表征了抑制劑分子中極性基團(如羥基或甲基磺酰胺基等)與NS5B發(fā)生靜電相互作用及形成氫鍵的能力;親水性指數(shù)(Hiwpl)描述抑制劑分子與NS5B酶的親水疏水相互作用;兩個量子化學(xué)描述符,分子最大負(fù)電荷(AQ,min)和分子平均負(fù)電荷(Mnc)可能與苯并噻嗪環(huán)的電子離域化及C-2,C-5位上吸電基團的引入有關(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 逐步回歸分析法所選描述符的建模回歸結(jié)果

我們采用逐步回歸分析法所選的6個最優(yōu)描述符分別建立多元線性、偏最小二乘和支持向量機回歸模型,訓(xùn)練集(含45個化合物)用來訓(xùn)練和優(yōu)化模型,測試集(含44個小分子抑制劑)評價回歸模型的預(yù)測能力.訓(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為回歸分析的顯著性檢驗.圖S1(Supporting Information)為多元線性回歸模型對于訓(xùn)練集和測試集的實驗與預(yù)測pIC50數(shù)值對比圖,模型對于訓(xùn)練集、測試集及整個數(shù)據(jù)集的均方誤差(MSE)分別為0.121、0.122和0.122,相應(yīng)的相關(guān)系數(shù)分別為0.970、0.958和0.965.偏最小二乘方法抽選出3個主成分,訓(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為偏最小二乘回歸模型對于訓(xùn)練集和測試集的實驗與預(yù)測pIC50數(shù)值對比圖,模型對于訓(xùn)練集、測試集及整個數(shù)據(jù)集的MSE分別為0.140、0.122和0.131,相應(yīng)的相關(guān)系數(shù)分別為0.965、0.958和0.961.由于影響NS5B酶抑制劑活性的分子特征極其復(fù)雜,并不是所有的分子描述符都與活性存在絕對的線性關(guān)系,為了與上述線性回歸建模方法做比較,基于相同的最優(yōu)描述符子集,我們使用遺傳算法組合支持向量機方法建立非線性QSAR模型.首先,我們在訓(xùn)練集上使用GASVM方法同時優(yōu)化SVM模型的三個參數(shù),當(dāng)SVM模型參數(shù)取值為:C=7.7387,σ=0.3546,ε=0.04664時,回歸模型給出最好的留一法交叉驗證誤差(MSECV=0.135304);然后,我們采用這組優(yōu)化參數(shù)建立SVM模型.圖S3為SVM回歸模型對于訓(xùn)練集和測試集的實驗與預(yù)測pIC50數(shù)值對比圖,模型對于訓(xùn)練集、測試集及整個數(shù)據(jù)集的MSE分別為0.113、0.108和0.111,相應(yīng)的相關(guān)系數(shù)分別為0.972、0.962和0.968.在上述三種方法建立的QSAR模型中,支持向量機給出最好的回歸模型,而多元線性回歸給出預(yù)測效果最好的線性回歸模型,三種機器學(xué)習(xí)模型對抑制劑活性的預(yù)測結(jié)果詳見表S1.

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

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

4 結(jié)論

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

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.