劉 英，尹京晨，計(jì)惠芬，賈青竹

(1. 天津科技大學(xué)海洋與環(huán)境學(xué)院，天津 300457；2. 浙江東洋環(huán)境工程有限公司，湖州 313000)

單端孢霉烯族毒素(trichothecenes)是一類化學(xué)性質(zhì)相近的真菌毒素，由鐮刀菌屬、單端孢屬和頭孢霉屬等真菌產(chǎn)生[1-2]，其結(jié)構(gòu)屬于倍半枯環(huán)氧化物，具有1個(gè)三環(huán)核心，在C12和C13位含有1個(gè)環(huán)氧化物，這是一個(gè)毒性必需基團(tuán)(核心結(jié)構(gòu)見圖 1)．該類毒素不僅污染小麥和玉米等禾谷類作物，也會(huì)通過食物鏈危害肉、奶、蛋等畜產(chǎn)品，WTO已將其列為天然存在的最危險(xiǎn)的食品污染物之一[3]．而且，由于該類毒素對(duì)真核細(xì)胞具有多重抑制作用，通過對(duì)蛋白質(zhì)、DNA 和 RNA 合成的抑制，對(duì)線粒體功能、細(xì)胞分裂和膜功能的抑制以及促進(jìn)凋亡基因的表達(dá)，對(duì)人畜健康產(chǎn)生免疫抑制，導(dǎo)致其具有典型的“致癌、致畸、致突變”作用，從而嚴(yán)重威脅人畜健康[4]．

圖1 單端孢霉烯核結(jié)構(gòu)Fig. 1 Structure of the trichothecene nucleus

定量構(gòu)效關(guān)系(quantitative structure-activity rela-tionship，QSAR)是研究化合物結(jié)構(gòu)與生物活性相關(guān)性的典型方法，廣泛應(yīng)用于新藥物設(shè)計(jì)以及對(duì)于毒性機(jī)制的理解探究領(lǐng)域[5-8]．盡管已有學(xué)者對(duì)單端孢霉烯族毒素的結(jié)構(gòu)與活性關(guān)系進(jìn)行了考察，但是至今關(guān)于這一系列毒素的 QSAR 研究依然很少[5-6].Steinmetz等[6]基于分子比較力場(chǎng)分析(comparative molecular field analysis，CoMFA)方法和比較分子相似因子分析(comparative molecular similarity indices analysis，CoMSIA)方法對(duì)32個(gè)大環(huán)內(nèi)酯物單端孢霉烯族毒素和 15個(gè)單端孢霉烯族毒素的毒性進(jìn)行了3D(3 dimension)QSAR研究，結(jié)果表明單端孢霉烯族毒素類物質(zhì)的取代基以及分子構(gòu)象對(duì)其生物活性都有一定影響．

本課題組之前提出系列范數(shù)指數(shù)描述符已經(jīng)成功地預(yù)測(cè)了有機(jī)物多種活性，包括麻醉性污染物的水生毒性、雜環(huán)化合物的藥理學(xué)和毒理學(xué)活性等[7-10].先前的研究工作表明該系列范數(shù)描述符也許可以從本質(zhì)上體現(xiàn)分子結(jié)構(gòu)，從而反映化合物的多個(gè)物化性質(zhì).

在本課題組之前工作的基礎(chǔ)上，本工作又提出新的描述符，并據(jù)此建立單端孢霉烯真菌毒素的毒性預(yù)測(cè)模型，對(duì) 35個(gè)單端孢霉烯真菌毒素的毒性進(jìn)行了計(jì)算，采用留一法交叉驗(yàn)證(LOO)、外部驗(yàn)證和 Y隨機(jī)驗(yàn)證手段對(duì)模型進(jìn)行了驗(yàn)證，并對(duì)模型的應(yīng)用域(AD)進(jìn)行了評(píng)價(jià)．本工作深入研究此類化合物分子結(jié)構(gòu)對(duì)其毒性的影響，將有助于人們采取更科學(xué)的毒素防控及清除機(jī)制．

1 材料與方法

1.1 樣本集

Jarvis等[11-12]測(cè)定了大環(huán)內(nèi)酯物真菌毒素在延長(zhǎng)接種p-388白血病小鼠生命中的能力，這些數(shù)據(jù)可以從根本上反映大環(huán)內(nèi)酯物真菌毒素對(duì)于小鼠白血病細(xì)胞的影響，并將測(cè)定活性定量為R，即R＝100×測(cè)試集動(dòng)物存活天數(shù)/對(duì)照組動(dòng)物存活天數(shù)．在此工作中，包含大環(huán)內(nèi)酯物真菌毒素毒性數(shù)據(jù)的化合物分子(代表性分子結(jié)構(gòu)見圖2)從文獻(xiàn)中[6,13-14]獲得．

圖2 兩種大環(huán)內(nèi)酯物真菌毒素結(jié)構(gòu)Fig. 2 Structure of two macrolide trichothecene mycotoxins

1.2 分子結(jié)構(gòu)優(yōu)化

利用軟件HyperChem 7.0構(gòu)建分子的3D結(jié)構(gòu)，對(duì)有機(jī)物分子結(jié)構(gòu)進(jìn)行優(yōu)化，具體采用量子化學(xué)中的從頭算法ab initio在ST0-3G進(jìn)行能量最低優(yōu)化，分子結(jié)構(gòu)信息中包括分子中原子的種類、數(shù)量、分支度，原子之間的鍵連接形式以及穩(wěn)定結(jié)構(gòu)下的原子電荷．

1.3 范數(shù)指數(shù)描述符和模型

在分子結(jié)構(gòu)優(yōu)化基礎(chǔ)上，為了區(qū)分原子的特點(diǎn)和種類，本工作提出了性質(zhì)矩陣．同時(shí)，歐氏空間距離矩陣是一種能夠全面、有效的描述分子的方法，因此，本工作提出了不同的矩陣．

下面列出上述具體矩陣．

為了討論上述矩陣的模，本工作使用了兩類范數(shù)指數(shù) norm(TP，1)和 norm(TP，2)，其定義式如下：norm(TP，1)：Fro函數(shù)

在式(8)中，操作符(.×)的含義是標(biāo)量相乘．通過范數(shù)計(jì)算得到48個(gè)范數(shù)描述符，為避免描述符過多而造成模型過度擬合，對(duì)建模所用描述符進(jìn)行優(yōu)化篩選，最后優(yōu)選了6個(gè)范數(shù)描述符，建立預(yù)測(cè)模型如下：

模型中系數(shù)b見表1．

表1 模型系數(shù)Tab. 1 Parameters of this model

1.4 模型評(píng)價(jià)

通過下列參數(shù)評(píng)價(jià)本工作建立模型的穩(wěn)定性：相關(guān)系數(shù)的平方(R2)、平均相對(duì)誤差值．為了進(jìn)一步驗(yàn)證模型的準(zhǔn)確性，還使用了留一交叉驗(yàn)證法(Q2)、外部驗(yàn)證法、Y隨機(jī)驗(yàn)證法進(jìn)行驗(yàn)證．

同時(shí)，本工作選擇了應(yīng)用域驗(yàn)證本模型的應(yīng)用范圍[15]，其定 義為式(15)：

在上述等式中，h是杠桿值，x是變量值．

通過與閾值(h*)比較，如果h＜h*，表明在這個(gè)區(qū)域內(nèi)模型的預(yù)測(cè)性是穩(wěn)定的；反之，這個(gè)模型的預(yù)測(cè)能力是較差的[16-17]．閾值的定義為式(16)

式中：p′是自變量數(shù)量加1；n是訓(xùn)練集數(shù)量．

2 結(jié)果與討論

2.1 大環(huán)內(nèi)酯物真菌毒素毒性的QSAR模型

利用新建模型式(9)對(duì) 35個(gè)大環(huán)內(nèi)酯物真菌毒素的毒性進(jìn)行了預(yù)測(cè)，圖3是大環(huán)內(nèi)酯物真菌毒素毒性的實(shí)驗(yàn)值log10R與預(yù)測(cè)值對(duì)比散點(diǎn)圖．

圖3 log10R預(yù)測(cè)值和實(shí)驗(yàn)值相關(guān)性Fig. 3 Correlation between log10R predicted by this model and the experimental data

由圖3可知：所有毒素的毒性預(yù)測(cè)點(diǎn)和實(shí)驗(yàn)點(diǎn)均位于對(duì)角線上及附近，表明本模型計(jì)算結(jié)果與實(shí)驗(yàn)值有很好的一致性．本預(yù)測(cè)模型相關(guān)統(tǒng)計(jì)數(shù)據(jù)、和 F值分別為 0.925,4、0.996,6和 78.95，說明本模型計(jì)算結(jié)果的精確性較好．

圖 4是殘差與實(shí)驗(yàn)值的對(duì)比圖．從圖 4可以看出所有毒素的毒性預(yù)測(cè)殘差都分布在-0.3～0.3，這表明本工作的系統(tǒng)誤差小，模型預(yù)測(cè)能力穩(wěn)定；同時(shí)，可以說明本工作提出的變量能有效地反映分子結(jié)構(gòu)與毒性之間的關(guān)系．

圖4 殘差與實(shí)驗(yàn)值對(duì)比圖Fig. 4 Plot of the residual vs．experimental data from this model

本工作與Steinmetz等[6]對(duì)比結(jié)果見表2．由表 2可以看出文獻(xiàn)中 CoMSIA方法和本工作得到的相關(guān)系數(shù) R2都大于 0.900，R2均在同一數(shù)量級(jí)，但是本工作交叉驗(yàn)證的相關(guān)系數(shù)遠(yuǎn)大于 CoMSIA方法，說明本工作的模型更穩(wěn)定．對(duì)于F值而言，本工作的F值同樣遠(yuǎn)大于文獻(xiàn)中兩種方法得出的 F值，更進(jìn)一步證明本模型的準(zhǔn)確性高、穩(wěn)定性好．

表2 本方法和文獻(xiàn)方法預(yù)測(cè)結(jié)果對(duì)比Tab. 2 Statistics comparison of reference methods and this work

本工作采用多元線性回歸分析法，建立關(guān)于大環(huán)內(nèi)酯物真菌毒素的結(jié)構(gòu)與毒性關(guān)系模型．Steinmetz等[6]利用CoMFA方法研究結(jié)果表明，空間位阻和靜電因素對(duì)毒素的活性有重大影響，其中靜電因素的貢獻(xiàn)更大；同時(shí)，相比之下包含氫鍵供體/受體的CoMSIA方法計(jì)算結(jié)果更為精確．本工作在計(jì)算描述符的優(yōu)化確定中，首先利用歐氏距離以及相關(guān)距離矩陣考慮了大環(huán)內(nèi)酯物真菌毒素分子中各個(gè)組成原子的空間分布；同時(shí)，為了實(shí)現(xiàn)大環(huán)內(nèi)酯物真菌毒素分子中不同原子對(duì)其毒性的貢獻(xiàn)影響，提出了原子屬性矩陣即增廣矩陣，其中包括電負(fù)性、原子質(zhì)量、原子支化度和原子電荷．由此，根據(jù)以上相關(guān)矩陣建立的范數(shù)描述符能從分子水平上反映真菌毒素分子結(jié)構(gòu)特點(diǎn)以及真菌毒素與生物體之間的相互作用力．綜上，與前人研究相比，本工作優(yōu)點(diǎn)在于：能給出確定的模型形式，方便其他研究者驗(yàn)證或使用；模型預(yù)測(cè)準(zhǔn)確性和穩(wěn)定性強(qiáng)．

在此以第13個(gè)化合物為例，計(jì)算該化合物的毒性.化合物Baccharin B8的結(jié)構(gòu)如圖5所示．

首先，按照式(1—3)的定義式，根據(jù)化學(xué)圖計(jì)算出該化合物的矩陣 T1、T2和性質(zhì)矩陣 PE；然后依據(jù)式(6—8)，計(jì)算出增廣矩陣 TP1,a,b、TP2,a,b和 TP3,a,b,.表 3是上述所有矩陣的范數(shù)描述符的值．應(yīng)用式(9)以及表2和表3中的參數(shù)，第13個(gè)化合物的毒性預(yù)測(cè)值計(jì)算過程如下：

第 13個(gè)化合物的 log10R預(yù)測(cè)值是 2.363，實(shí)驗(yàn)值是2.367，相對(duì)偏差是0.17%,．

圖5 Baccharin B8的結(jié)構(gòu)Fig. 5 Structure of Baccharin B8

表3 第13個(gè)化合物的相關(guān)變量值Tab. 3 Variable values of the thirteenth chemical

2.2 留一交叉驗(yàn)證

留一交叉驗(yàn)證法就是從所有的數(shù)據(jù)中每次選擇一個(gè)分子作為一次測(cè)試集，用剩余的分子建立模型預(yù)測(cè)測(cè)試集分子的毒性．如果留一交叉驗(yàn)證的相關(guān)系數(shù)大于0.5，則認(rèn)為模型是合理的[18]．圖6表明：留一交叉驗(yàn)證法的毒性預(yù)測(cè)值與實(shí)驗(yàn)值有較好吻合度，表明本工作建立的模型的內(nèi)部預(yù)測(cè)能力較高．

圖6 留一交叉驗(yàn)證預(yù)測(cè)值和實(shí)驗(yàn)值相關(guān)性Fig. 6 Correlation between the data predicted by LOOCV and the experimental data

圖 7是兩種模型預(yù)測(cè)結(jié)果的標(biāo)準(zhǔn)誤差(SE)分布情況對(duì)比圖．由圖 7可知：交叉驗(yàn)證的誤差趨勢(shì)和本方法建立的模型的誤差趨勢(shì)一致，說明本工作建立的模型穩(wěn)定、預(yù)測(cè)能力好．同時(shí)，留一交叉驗(yàn)證法的結(jié)果具有較高 Q2值(0.878,9)，以上所有驗(yàn)證結(jié)果均表明本工作提出的方法可以穩(wěn)定地預(yù)測(cè)大環(huán)內(nèi)酯物真菌毒素的毒性．

圖7 本模型和LOO模型毒性預(yù)測(cè)標(biāo)準(zhǔn)誤差分布Fig. 7 Standard error distributions of the toxicity predicted by this model and LOO-CV

2.3 外部驗(yàn)證

模型的內(nèi)部驗(yàn)證只能說明該方法對(duì)于類似物的活性預(yù)測(cè)能力好，如果用該方法預(yù)測(cè)新的不同種類的化合物的活性，則結(jié)果有可能是不可靠的．因此，需要進(jìn)行外部驗(yàn)證．外部驗(yàn)證就是將全部的分子分為訓(xùn)練集和測(cè)試集，訓(xùn)練集用于建立模型，測(cè)試集用來驗(yàn)證模型．有效的外部驗(yàn)證能確保本工作提出的方法對(duì)測(cè)試集分子預(yù)測(cè)的可預(yù)見性和適用性．

在本工作中，訓(xùn)練集和測(cè)試集的分子與文獻(xiàn)[6]相同，訓(xùn)練集和測(cè)試集的分子毒性預(yù)測(cè)值見圖 8，其結(jié)果是0.925,4、是 0.996,6．這表明預(yù)測(cè)值和實(shí)驗(yàn)值有很好的一致性，進(jìn)而說明本工作建立的模型不僅可以很好地預(yù)測(cè)同類化合物毒性，而且對(duì)多種類化合物具有較高的預(yù)測(cè)性．

圖8 訓(xùn)練集和測(cè)試集的實(shí)驗(yàn)值和預(yù)測(cè)值的散點(diǎn)圖Fig. 8 Scatter plot predicted and experimental data of the training set and test set

2.4 Y隨機(jī)驗(yàn)證

Y隨機(jī)驗(yàn)證的原理是將實(shí)驗(yàn)值隨機(jī)擾亂、變量 不變，再建立模型，得到的結(jié)果與最初建立模型的結(jié)果進(jìn)行對(duì)比．

在本工作中，進(jìn)行了 5次隨機(jī)驗(yàn)證，新模型預(yù)測(cè)R2和 Q2的結(jié)果與原模型的對(duì)比結(jié)果見表 4．由表 4可以得出，5次Y隨機(jī)驗(yàn)證得出的新模型的R2和Q2都很低，甚至為0，說明新模型的預(yù)測(cè)能力低．由此證明，本工作建立的原始模型是穩(wěn)定的，并非偶然建立.

表4 Y隨機(jī)驗(yàn)證結(jié)果Tab. 4 Results of Y randomization test of this model

2.5 應(yīng)用域驗(yàn)證

對(duì)于構(gòu)效關(guān)系模型統(tǒng)計(jì)參數(shù)的驗(yàn)證來說，模型的應(yīng)用域評(píng)價(jià)依然是最困難的一個(gè)評(píng)價(jià)標(biāo)準(zhǔn)．應(yīng)用域評(píng)價(jià)是模型能否準(zhǔn)確預(yù)測(cè)新化合物性質(zhì)或活性的重要指標(biāo)．目前，已經(jīng)有多種關(guān)于模型應(yīng)用域評(píng)價(jià)的方法[19-21]．本工作利用杠桿方法檢測(cè)計(jì)算模型的應(yīng)用域，結(jié)果如圖9所示．從圖 9可以看出，對(duì)于35個(gè)有機(jī)物，只有 1個(gè)有機(jī)物的杠桿值小于-3，其余所有分子的杠桿值均在規(guī)定的閾值內(nèi)，說明本方法建立的模型應(yīng)用域較廣泛．

圖9 本模型的應(yīng)用域Fig. 9 Application domain of this model

3 結(jié) 語

本工作基于有機(jī)毒素化學(xué)結(jié)構(gòu)構(gòu)建的距離矩陣和性質(zhì)矩陣，提出了一個(gè)新范數(shù)描述符，并據(jù)此建立了預(yù)測(cè)單端孢霉烯真菌毒素毒性的 QSAR模型．研究結(jié)果表明：本工作提出的范數(shù)描述符能夠有效描述化合物毒性與其結(jié)構(gòu)之間的關(guān)系；本模型毒性預(yù)測(cè)值與實(shí)驗(yàn)值有很好的一致性，其中模型的相關(guān)性系數(shù)R2均大于 0.9，F(xiàn)值為 78.95，留一交叉驗(yàn)證測(cè)試(Q2為 0.878,9)、Y隨機(jī)驗(yàn)證以及與文獻(xiàn)的對(duì)比結(jié)果均證明本模型穩(wěn)定性更好、預(yù)測(cè)能力更高、準(zhǔn)確性更好.同時(shí)，應(yīng)用域驗(yàn)證表明本模型計(jì)算結(jié)果準(zhǔn)確、穩(wěn)定、可靠、適用范圍廣．因此，本工作提出的方法將可能推廣用于化學(xué)品的環(huán)境風(fēng)險(xiǎn)評(píng)價(jià)、生態(tài)健康評(píng)價(jià)等領(lǐng)域.

參考文獻(xiàn)：

［1］ Zhou T，He J W，Gong J. Microbial transformation of trichothecene mycotoxins[J]. World Mycotoxin Journal，2008，1(1)：23-30.

［2］ 鄒忠義，賀稚非，李洪軍，等. 單端孢霉烯族毒素轉(zhuǎn)化降解研究進(jìn)展[J]. 食品科學(xué)，2010，31(19)：443-448.

［3］ 周闖，何成華，司慧民，等. 2012年國內(nèi)飼料及原料霉菌毒素污染調(diào)查分析[J]. 畜牧與獸醫(yī)，2014，46(1)：81-84.

［4］ 許偉，耿芳芳，范夢(mèng)雪，等. 脫氧雪腐鐮刀菌烯醇毒性的研究進(jìn)展[J]. 生物學(xué)雜志，2016，33(1)：78-81.

［5］ Betina V. Structure-activity relationships among mycotoxins[J]. Chemico-Biological Interactions，1989，71(2/3)：105-146.

［6］ Steinmetz W E，Rodarte C B，Lin A. 3D QSAR study of the toxicity of trichothecene mycotoxins[J]. European Journal of Medicinal Chemistry，2009，44(11)：4485-4489.

［7］ Zhu Z，Wang Q，Jia Q，et al. Structure-property relationship for the pharmacological and toxicological activity of heterocyclic compounds[J]. Acta Physico-chimica Sinica，2014，30(6)：1086-1090.

［8］ Qian H，Kanwal S，Jia Q，et al. Norm index-based quantitative structure-activity relationship to predict βcyclodextrin complex binding constants[J]. Acta Physico-chimica Sinica，2015，31(5)：893-898.

［9］ Wang Q，Jia Q，Yan L，et al. Quantitative structuretoxicity relationship of the aquatic toxicity for various narcotic pollutants using the norm indexes[J]. Chemosphere，2014，108：383-387.

［10］ Jia Q，Cui X，Li L，et al. A Quantitative structure-activity relationship for high affinity 5-HT1A receptor ligands based on norm indexes[J]. Journal of Physical Chemistry B，2015，119(51)：15561-15567.

［11］ Jarvis B B，Stahly G P，Pavanasasivam G. Isolation and characterization of the trichoverroids and new roridins and verrucarins[J]. Journal of Medicinal Chemistry，1982，47(6)：1117-1124.

［12］ Jarvis B B，Vrudhula V M，Midiwo J O，et al. New trichoverroids from myrothecium verrucaria：Verrol and 12，13-deoxytrichodermadiene[J]. Journal of Medicinal Chemistry，1983，48(15)：2576-2578.

［13］ Cramer R D，Patterson D E，Bunce J D. Comparative molecular field analysis(CoMFA). 1. Effect of shape on binding of steroids to carrier proteins[J]. Journal of the American Chemical Society，1988，110(18)：5959-5967.

［14］ Klebe G，Abraham U，Mietzner T. Molecular similarity indices in a comparative analysis(CoMSIA)of drug molecules to correlate and predict their biological activity[J]. Journal of Medicinal Chemistry，1994，37(24)：4130-4146.

［15］ Netzeva T I，Worth A P，Aldenberg T，et al. Current status of methods for defining the applicability domain of(quantitative)structure-activity relationships[J]. American Theological Library Association，2005，33(2)：155-173.

［16］ Gramatica P. Principles of QSAR models validation：Internal and external[J]. QSAR & Combinatorial Science，2007，26(5)：694-701.

［17］ Gramatica P，Giani E，Papa E. Statistical external validation and consensus modeling：A QSPR case study for Kocprediction[J]. Journal of Molecular Graphics and Modelling，2007，25(6)：755-766.

［18］ Mihali? Z，Trinajsti? N. A graph-theoretical approach to structure-property relationships[J]. Journal of Chemical Education，1992，69(9)：701.

［19］ Tropsha A. Best practices for QSAR model development，validation，and exploitation[J]. Molecular Informatics，2010，29(6/7)：476-488.

［20］ Tropsha A，Golbraikh A. Predictive QSAR modeling workflow，model applicability domains，and virtual screening[J]. Current Pharmaceutical Design，2007，13(34)：3494-3504.

［21］ Jaworska J，Nikolova-Jeliazkova N，Aldenberg T. QSAR applicabilty domain estimation by projection of the training set descriptor space：A review[J]. American Theological Library Association，2005，33(5)：445-459.