徐靜安 賀少鵬 商照聰

多年來參與一些項目的評審，研究生論文的答辯，相關(guān)專業(yè)文獻的查閱等，筆者感受到隨著計算機軟硬件技術(shù)的發(fā)展，20世紀80年代以來，化合物的定量結(jié)構(gòu)-活性/性質(zhì)相關(guān)性（簡稱構(gòu)效關(guān)系，英文縮寫QSAR/QSPR）逐漸成為研究的熱點。2018年2月7日上午，筆者應(yīng)邀參加在華誼集團大廈舉辦的煤基多聯(lián)產(chǎn)工程中心和計算化學化工工程中心技術(shù)委員會的年度會議，會上，計算化學也受到了工程界的重視。

上?；ぱ芯吭嚎蒲泄ぷ饕策M行了相應(yīng)的安排、探索。2006級碩士研究生賀少鵬，其學位論文“有機污染物的正辛醇/水分配系數(shù)預(yù)測及QSPR研究”，導師是徐大剛、劉剛二位教授級高工。筆者作為研究生辦公室顧問一直跟蹤項目和參與討論。2009年3月筆者還從賀少鵬處借閱《化學計量學方法》，閱讀后2011年網(wǎng)購了幾本贈與有關(guān)專業(yè)人員學習、應(yīng)用。2012年又閱讀了商照聰、賀少鵬的論文“有機污染物分配系數(shù)（正辛醇/水）預(yù)測軟件比較研究”。

2013年院部采購了IBM高性能小型機；2014年，購置了DPS軟件，還和上海應(yīng)用技術(shù)大學共建共享配置了VASP軟件；2015年院部購置了Gaussian軟件。此外，研究院還在2013年和2014年針對性地招聘了量子化學軟件應(yīng)用的專業(yè)人員。

在材料、生物、環(huán)境等科學領(lǐng)域，構(gòu)效關(guān)系研究在宏觀、介觀及微觀層面展開，本文討論的是分子尺度的化合物構(gòu)效關(guān)系。

一、建立構(gòu)效關(guān)系模型的主要步驟

根據(jù)《有機污染化學》（王連生編著），建立QSPR/QSAR模型的主要步驟見圖1。具體如下：

1 訓練集化合物的選擇

圖1 建立QSAR/QSPR模型的主要步驟

根據(jù)一定的統(tǒng)計標準和結(jié)構(gòu)標準選擇類系化合物，作為建模的訓練集?；衔镞x擇的條件為：統(tǒng)計上的隨機性、結(jié)構(gòu)上的代表性和全面性，以及性質(zhì)/活性數(shù)據(jù)的可獲得性。

2 性質(zhì)/活性數(shù)據(jù)的獲得

數(shù)據(jù)收集的方法主要有3種：科學文獻的查閱、收集；性質(zhì)/活性數(shù)據(jù)庫中獲取；實驗室中測試。一般來說，性質(zhì)/活性數(shù)據(jù)有兩種：連續(xù)響應(yīng)數(shù)據(jù)（例如水溶解度等）和非連續(xù)分類、分級響應(yīng)數(shù)據(jù)（例如致癌性的陽性/陰性、毒性等級等）。

3 分子結(jié)構(gòu)描述符

首先應(yīng)用分子模擬方法，構(gòu)建正確的二維或三維分子結(jié)構(gòu)，采用構(gòu)象分析、分子力學等方法獲得最優(yōu)化的構(gòu)象；采用拓撲學方法、量子力學方法等計算化合物的分子結(jié)構(gòu)特征參數(shù)，獲取分子的結(jié)構(gòu)信息，進行分子結(jié)構(gòu)描述。

4 QSAR模型構(gòu)建

應(yīng)用特征變量篩選方法篩選描述符，應(yīng)用統(tǒng)計方法或其他建模方法將訓練集有機化合物的性質(zhì)/活性數(shù)據(jù)與分子結(jié)構(gòu)參數(shù)數(shù)據(jù)建立QSAR模型。常用的統(tǒng)計分析方法有回歸分析、偏最小二乘分析、因子分析、主成分分析、聚類分析等，其中多元逐步回歸分析是應(yīng)用最多的方法。近年來，人工神經(jīng)網(wǎng)絡(luò)和遺傳算法等高級建模方法也得到了關(guān)注和應(yīng)用。

5 QSAR模型的評價和驗證

模型的評價主要針對：模型的擬合優(yōu)度、穩(wěn)健性和預(yù)測能力。在回歸分析中，模型的擬合優(yōu)度采用回歸系數(shù)的平方（R2）或自由度校正的R2（R2edj）、顯著性水平、檢驗值F、標準偏差s等參數(shù)來評判。模型的穩(wěn)健性一般采用交叉驗證方法來進行檢驗，通常有兩種方式：逐一剔除法（即留一法）和分組剔除法（即留多法）。得到的交叉驗證的R2（q2）和標準預(yù)測誤差（SEP）用來評價模型的穩(wěn)健性。模型預(yù)測的驗證是構(gòu)建一個測試集，用訓練集建立的擬合模型來預(yù)測測試集化合物的性質(zhì)。只有具有統(tǒng)計上的顯著性、穩(wěn)健以及具有高度預(yù)測能力的模型才能夠進行應(yīng)用。

6 QSAR模型的應(yīng)用

應(yīng)用有三個方面：一利用模型對其他未知化合物的相關(guān)性質(zhì)/活性進行預(yù)測，在效應(yīng)評價和暴露評價等方面可彌補缺失的數(shù)據(jù)，對有機化合物進行篩選和評價；二根據(jù)模型的組成與形式，結(jié)合已有的化學、生物學知識，探求有機化合物的毒理性質(zhì)、環(huán)境過程和生態(tài)效應(yīng)等機理分析；三根據(jù)所闡明的結(jié)構(gòu)-性質(zhì)關(guān)系結(jié)果，為設(shè)計目標化合物指明方向。

二、案例一

選取自《化學計量學方法》?；衔锏慕Y(jié)構(gòu)是個三維圖像，對于化合物分子尺度構(gòu)效關(guān)系的研究，就是將結(jié)構(gòu)圖特征參數(shù)與其性質(zhì)/活性去構(gòu)造相關(guān)的數(shù)學模型。主要參數(shù)類型有（1）拓撲類參數(shù)；（2）幾何類參數(shù)；（3）電子類參數(shù)；（4）理化性質(zhì)類參數(shù)；（5）綜合類參數(shù)；等。構(gòu)效關(guān)系模型建立的目的是對新的未知樣本進行預(yù)測。

案例為50個酚類化合物麻醉毒性的QSAR研究，樣本量N=50，提取可能影響麻醉毒性的結(jié)構(gòu)特征參數(shù)有M=135個變量，經(jīng)過變量篩選獲得構(gòu)效關(guān)系的統(tǒng)計模型、多元線性回歸方程：

式中：SHDWi——分子投影面積；

按全回歸分析方法，自變量M=135，樣本量N=50，此回歸模型無法求解。由于有的自變量——結(jié)構(gòu)特征參數(shù)對其響應(yīng)參數(shù)不顯著，特別是分子描述符參數(shù)之間普遍存在共線性關(guān)系，所以結(jié)構(gòu)特征參數(shù)的提取、篩選是構(gòu)效關(guān)系研究的關(guān)鍵。應(yīng)用數(shù)據(jù)處理方法，本例經(jīng)篩選進入模型的自變量m=9，大大簡化了模型。本例M=135，選用最簡單的多元線性回歸，可能構(gòu)成的模型有2M-1=2135-1；如果選用二次多項式回歸，僅考慮一次項和二次項，可能構(gòu)成模型有（22M-1）個。采用窮舉法時計算工作難以操作，由此研究產(chǎn)生若干變量篩選算法。

傳統(tǒng)的變量篩選方法有前進法、后退法、逐步回歸法、最佳子集法；常用的逐步回歸法已可有效篩選變量。在多元回歸分析中，亦有使用主成分分析法、正交變換法篩選變量。如果多元線性回歸建模效果不好，研究對象較為復雜，基于MIV的人工神經(jīng)網(wǎng)絡(luò)法、自變量降維的遺傳算法，以及針對小樣本的支持向量回歸法（SVR）等值得關(guān)注。

作為應(yīng)用案例，本案例模型尚應(yīng)進一步進行預(yù)報測試的評估與驗證檢驗。

三、案例二

選取自“有機污染物的正辛醇/水分配系數(shù)預(yù)測及QSPR研究”。將美國國家環(huán)境保護局推薦的105種優(yōu)先毒性污染物作為考察對象，即樣本量N=105。按化學結(jié)構(gòu)可分為鹵代（烷、烯）烴類，苯系物，酚類，多環(huán)芳烴類，亞硝胺類等10個類系化合物。響應(yīng)分配系數(shù)實驗值選自SRC公司的PHYSPROP數(shù)據(jù)庫。

資料表明，對于分子描述符現(xiàn)已開發(fā)出上百種拓撲指數(shù)，原文列舉常用的拓撲指數(shù)38種，幾何結(jié)構(gòu)性質(zhì)描述符24種，量子化學描述符35種，溶劑化描述符13種及理化性質(zhì)10余種等等。原文用現(xiàn)有的采用不同特征參數(shù)的10種估算分配系數(shù)的軟件（ALOGPs，AC logP，AB/LogP，COSMOFrag，miLogP，ALOGP，MLOGP，KOWWIN，XLOGP 2，XLOGP 3） 對105 種污染物進行了計算分析。其中ALOGPs采用拓撲指數(shù)等描述符、KOWWIN采用碎片常數(shù)法進行預(yù)測估算，預(yù)測效果相對較好，其他軟件對此分配系數(shù)的性能預(yù)測偏差較大。為進一步改進提高有機污染物正辛醇/水分配系數(shù)logP的預(yù)測性能，原文探索篩選新的特征參數(shù)構(gòu)建QSPR新的模型。

分子描述符的計算提取。先使用ChemBio3D Ultra11.0軟件將105種有機物的二維分子結(jié)構(gòu)轉(zhuǎn)化為三維結(jié)構(gòu)式，并使用內(nèi)置的MM2分子優(yōu)化軟件計算出能量最低的狀態(tài)，進行能量優(yōu)化，使用MOPAC軟件進行幾何構(gòu)型優(yōu)化。再使用Chem3D軟件計算29種分子結(jié)構(gòu)描述符，其中有：分子面積、橢圓度等7種結(jié)構(gòu)性質(zhì)描述符；總能量、生成熱等8種量子化學描述符；分子拓撲指數(shù)、分子半徑等14種拓撲描述符。此外，根據(jù)研究文獻表明有機物的分配系數(shù)與其水溶解度logSw和相對分子質(zhì)量Mw有很大關(guān)系，相對分子質(zhì)量根據(jù)結(jié)構(gòu)式計算可得，水溶解度數(shù)據(jù)取自SRC公司的PHYSPROP數(shù)據(jù)庫。共提取特征參數(shù)M=31個。

對特征參數(shù)進行初篩選一般采用變量零值檢測、偏差測試、相關(guān)性測試、共線性測試。由于多元線性逐步回歸同時完成變量剔除和模型建立，本例采用此主體技術(shù)。

log Pow的QSPR模型：

（1）對7種結(jié)構(gòu)性質(zhì)描述符建模；

（2）對8種量子化學描述符建模；

（3）對14種拓撲描述符建模；

（4）對水溶解度、相對分子質(zhì)量建模；

（5）對31個描述符綜合建模。

經(jīng)多元線性逐步回歸，模型統(tǒng)計量匯總見表1。

表1 模型統(tǒng)計量匯總

由表1可見，綜合參數(shù)模型擬合性能最優(yōu)：

模型中：Ovality——分子橢圓度，

PM-Y——慣性主矩的y分量，

Bindx——Balaban指數(shù)，

PSA——分子極化表面面積，

logSw——水溶解度對數(shù)值，

Mw——相對分子質(zhì)量。

從數(shù)理統(tǒng)計角度，此構(gòu)效關(guān)系模型擬合性能具有顯著的統(tǒng)計意義，尚應(yīng)對模型預(yù)報性能進行評估驗證。

內(nèi)部驗證。采用留一法預(yù)測標準偏差。

預(yù)測相關(guān)系數(shù)平方：

預(yù)測和擬合相關(guān)系數(shù)及標準偏差變化不顯著，模型具有穩(wěn)定性。

外部驗證。對新考察的①萘、②苯甲酸甲酯、③苯甲酸乙酯、④苯乙酮、⑤二苯醚、⑥肉桂醇、⑦溴苯、⑧苯甲酸芐酯8種有機物，使用高效液相色譜實驗測定正辛醇/水分配系數(shù)，使用上述軟件計算6項參數(shù)。參數(shù)計算值、響應(yīng)預(yù)測值、實驗值見表2。

外部驗證的標準偏差

驗證標準偏差相對穩(wěn)定，具有統(tǒng)計意義。2004年QSAR國際會議正式形成經(jīng)濟合作與發(fā)展組織（英文簡稱OECD）規(guī)則，明確必須使用外部驗證集（即測試集）來評價模型的預(yù)測能力。如果樣本量足夠大，也可以從105個樣本中隨機取8個樣本作為測試集，97個樣本作為訓練集。本案例執(zhí)行該規(guī)范。

四、案例三

選自上海交通大學環(huán)境科學與工程學院的“Fenton法降解有機物的熱力學及構(gòu)效關(guān)系研究”一文。2017年11月29日筆者在院內(nèi)彭東輝教授級高工實驗室討論高濃度有機廢水處理項目時閱讀到此文，值得推薦。

表2 外部驗證數(shù)據(jù)匯總

案例針對有機廢水中的①3,4-二氯苯胺、②對氨基苯磺酸、③ 2,4-二硝基苯肼、④ 雙酚A、⑤ 酸性橙、⑥間甲酚紫6種有機污染物，實驗探究了不同溫度下Fenton氧化降解有機物的去除率及反應(yīng)動力學；在Fenton試劑過量、假定反應(yīng)初期為一級反應(yīng)動力學速率常數(shù)的基礎(chǔ)上，通過Arrhenius方程計算獲得6種有機物的Fenton反應(yīng)活化能Ea。

本文側(cè)重討論了結(jié)構(gòu)參數(shù)與活化能之間的構(gòu)效關(guān)系。根據(jù)現(xiàn)有資料，案例利用Hyperchem8.0軟件初步優(yōu)化6個有機物的結(jié)構(gòu)，再用Gaussian09軟件進行深度優(yōu)化，利用Materials Studio軟件進行參數(shù)計算。選取了19個量子化學結(jié)構(gòu)參數(shù)，其中16個由軟件計算獲得，3 個為組合參數(shù)（EGAP，E2GAP，EHOMO+ELOMO），各參數(shù)的物理意義見表3，各參數(shù)值及與Ea的相關(guān)性見表4。

由于論文報道的僅僅是研究工作的階段內(nèi)容，樣本量為N=6，結(jié)構(gòu)特征參數(shù)M=19，還無法構(gòu)建構(gòu)效關(guān)系模型，但可以進行初步分析：

（1）先計算單一特征參數(shù)和活化能的相關(guān)關(guān)系，見表4。

表3 量子化學結(jié)構(gòu)參數(shù)

表4參數(shù)與Ea間相關(guān)系數(shù)大于0.5的項數(shù)11/19=58%，而且大于0.7的項數(shù)有2項，粗略判斷初選的結(jié)構(gòu)參數(shù)是有效的。如果相關(guān)性普遍較低，就應(yīng)考慮調(diào)整、擴充特征參數(shù)的選擇范圍。

（2）從案例一、二及相關(guān)資料分析可知，在構(gòu)效關(guān)系研究中，由于化合物分子特征參數(shù)對化合物性能響應(yīng)不同，很多參數(shù)對某一響應(yīng)是不顯著的；更為普遍的現(xiàn)象是顯著項特征參數(shù)間還存在共線性現(xiàn)象，所以M項特征參數(shù)經(jīng)篩選僅有限的m項進入構(gòu)效關(guān)系模型。

（3）所謂共線性是指成對自變量（特征參數(shù)）之間的相關(guān)性。當相關(guān)性較高時，表示一個變量的信息含有對應(yīng)變量的信息，可以剔除一個變量。如同時引入統(tǒng)計模型，除了增加計算工作量外，還會使模型計算性能變差。

表4 各參數(shù)值及與Ea的相關(guān)性

本講義第四講回歸分析中的變量篩選技術(shù)及統(tǒng)計檢驗（《上?；ぁ?、2016年第8期）已做詳細討論?，F(xiàn)再簡單引用，原案例是4個變量的多元線性回歸，計算機在逐步回歸建模時，自動計算輸出變量間的相關(guān)性，如表5所示。

x1與 x3相關(guān)系數(shù) R=-0.824，P=0.001＜0.05；x2與x4的 R=-0.973，P=0.000。可見：x1與 x3，x2與 x4均為顯著相關(guān)，建模時將作剔除處理。

（4）由于構(gòu)效關(guān)系研究中樣本選擇是隨機的，多維空間結(jié)構(gòu)參數(shù)的數(shù)值分布也是隨機的，構(gòu)建構(gòu)效關(guān)系模型一般要求N/m≥5，一定的樣本量保證模型的穩(wěn)定性，所以完成本案例構(gòu)效關(guān)系研究尚需一定的實驗工作量。

表5 四個變量間的相關(guān)性

后記：在本講義編寫中，賀少鵬提交了第一節(jié)的文稿并提供同濟大學的“QSAR模型內(nèi)部和外部驗證方法綜述”等3篇資料。筆者2017年已欣然為賀少鵬寫推薦信，現(xiàn)其正在此領(lǐng)域選題攻讀在職博士。