徐靜安 賀少鵬 商照聰
多年來參與一些項目的評審,研究生論文的答辯,相關專業文獻的查閱等,筆者感受到隨著計算機軟硬件技術的發展,20世紀80年代以來,化合物的定量結構-活性/性質相關性(簡稱構效關系,英文縮寫QSAR/QSPR)逐漸成為研究的熱點。2018年2月7日上午,筆者應邀參加在華誼集團大廈舉辦的煤基多聯產工程中心和計算化學化工工程中心技術委員會的年度會議,會上,計算化學也受到了工程界的重視。
上海化工研究院科研工作也進行了相應的安排、探索。2006級碩士研究生賀少鵬,其學位論文“有機污染物的正辛醇/水分配系數預測及QSPR研究”,導師是徐大剛、劉剛二位教授級高工。筆者作為研究生辦公室顧問一直跟蹤項目和參與討論。2009年3月筆者還從賀少鵬處借閱《化學計量學方法》,閱讀后2011年網購了幾本贈與有關專業人員學習、應用。2012年又閱讀了商照聰、賀少鵬的論文“有機污染物分配系數(正辛醇/水)預測軟件比較研究”。
2013年院部采購了IBM高性能小型機;2014年,購置了DPS軟件,還和上海應用技術大學共建共享配置了VASP軟件;2015年院部購置了Gaussian軟件。此外,研究院還在2013年和2014年針對性地招聘了量子化學軟件應用的專業人員。
在材料、生物、環境等科學領域,構效關系研究在宏觀、介觀及微觀層面展開,本文討論的是分子尺度的化合物構效關系。
根據《有機污染化學》(王連生編著),建立QSPR/QSAR模型的主要步驟見圖1。具體如下:

圖1 建立QSAR/QSPR模型的主要步驟
根據一定的統計標準和結構標準選擇類系化合物,作為建模的訓練集。化合物選擇的條件為:統計上的隨機性、結構上的代表性和全面性,以及性質/活性數據的可獲得性。
數據收集的方法主要有3種:科學文獻的查閱、收集;性質/活性數據庫中獲取;實驗室中測試。一般來說,性質/活性數據有兩種:連續響應數據(例如水溶解度等)和非連續分類、分級響應數據(例如致癌性的陽性/陰性、毒性等級等)。
首先應用分子模擬方法,構建正確的二維或三維分子結構,采用構象分析、分子力學等方法獲得最優化的構象;采用拓撲學方法、量子力學方法等計算化合物的分子結構特征參數,獲取分子的結構信息,進行分子結構描述。
應用特征變量篩選方法篩選描述符,應用統計方法或其他建模方法將訓練集有機化合物的性質/活性數據與分子結構參數數據建立QSAR模型。常用的統計分析方法有回歸分析、偏最小二乘分析、因子分析、主成分分析、聚類分析等,其中多元逐步回歸分析是應用最多的方法。近年來,人工神經網絡和遺傳算法等高級建模方法也得到了關注和應用。
模型的評價主要針對:模型的擬合優度、穩健性和預測能力。在回歸分析中,模型的擬合優度采用回歸系數的平方(R2)或自由度校正的R2(R2edj)、顯著性水平、檢驗值F、標準偏差s等參數來評判。模型的穩健性一般采用交叉驗證方法來進行檢驗,通常有兩種方式:逐一剔除法(即留一法)和分組剔除法(即留多法)。得到的交叉驗證的R2(q2)和標準預測誤差(SEP)用來評價模型的穩健性。模型預測的驗證是構建一個測試集,用訓練集建立的擬合模型來預測測試集化合物的性質。只有具有統計上的顯著性、穩健以及具有高度預測能力的模型才能夠進行應用。
應用有三個方面:一利用模型對其他未知化合物的相關性質/活性進行預測,在效應評價和暴露評價等方面可彌補缺失的數據,對有機化合物進行篩選和評價;二根據模型的組成與形式,結合已有的化學、生物學知識,探求有機化合物的毒理性質、環境過程和生態效應等機理分析;三根據所闡明的結構-性質關系結果,為設計目標化合物指明方向。
選取自《化學計量學方法》。化合物的結構是個三維圖像,對于化合物分子尺度構效關系的研究,就是將結構圖特征參數與其性質/活性去構造相關的數學模型。主要參數類型有(1)拓撲類參數;(2)幾何類參數;(3)電子類參數;(4)理化性質類參數;(5)綜合類參數;等。構效關系模型建立的目的是對新的未知樣本進行預測。
案例為50個酚類化合物麻醉毒性的QSAR研究,樣本量N=50,提取可能影響麻醉毒性的結構特征參數有M=135個變量,經過變量篩選獲得構效關系的統計模型、多元線性回歸方程:

式中:SHDWi——分子投影面積;

按全回歸分析方法,自變量M=135,樣本量N=50,此回歸模型無法求解。由于有的自變量——結構特征參數對其響應參數不顯著,特別是分子描述符參數之間普遍存在共線性關系,所以結構特征參數的提取、篩選是構效關系研究的關鍵。應用數據處理方法,本例經篩選進入模型的自變量m=9,大大簡化了模型。本例M=135,選用最簡單的多元線性回歸,可能構成的模型有2M-1=2135-1;如果選用二次多項式回歸,僅考慮一次項和二次項,可能構成模型有(22M-1)個。采用窮舉法時計算工作難以操作,由此研究產生若干變量篩選算法。
傳統的變量篩選方法有前進法、后退法、逐步回歸法、最佳子集法;常用的逐步回歸法已可有效篩選變量。在多元回歸分析中,亦有使用主成分分析法、正交變換法篩選變量。如果多元線性回歸建模效果不好,研究對象較為復雜,基于MIV的人工神經網絡法、自變量降維的遺傳算法,以及針對小樣本的支持向量回歸法(SVR)等值得關注。
作為應用案例,本案例模型尚應進一步進行預報測試的評估與驗證檢驗。
選取自“有機污染物的正辛醇/水分配系數預測及QSPR研究”。將美國國家環境保護局推薦的105種優先毒性污染物作為考察對象,即樣本量N=105。按化學結構可分為鹵代(烷、烯)烴類,苯系物,酚類,多環芳烴類,亞硝胺類等10個類系化合物。響應分配系數實驗值選自SRC公司的PHYSPROP數據庫。
資料表明,對于分子描述符現已開發出上百種拓撲指數,原文列舉常用的拓撲指數38種,幾何結構性質描述符24種,量子化學描述符35種,溶劑化描述符13種及理化性質10余種等等。原文用現有的采用不同特征參數的10種估算分配系數的軟件(ALOGPs,AC logP,AB/LogP,COSMOFrag,miLogP,ALOGP,MLOGP,KOWWIN,XLOGP 2,XLOGP 3) 對105 種污染物進行了計算分析。其中ALOGPs采用拓撲指數等描述符、KOWWIN采用碎片常數法進行預測估算,預測效果相對較好,其他軟件對此分配系數的性能預測偏差較大。為進一步改進提高有機污染物正辛醇/水分配系數logP的預測性能,原文探索篩選新的特征參數構建QSPR新的模型。
分子描述符的計算提取。先使用ChemBio3D Ultra11.0軟件將105種有機物的二維分子結構轉化為三維結構式,并使用內置的MM2分子優化軟件計算出能量最低的狀態,進行能量優化,使用MOPAC軟件進行幾何構型優化。再使用Chem3D軟件計算29種分子結構描述符,其中有:分子面積、橢圓度等7種結構性質描述符;總能量、生成熱等8種量子化學描述符;分子拓撲指數、分子半徑等14種拓撲描述符。此外,根據研究文獻表明有機物的分配系數與其水溶解度logSw和相對分子質量Mw有很大關系,相對分子質量根據結構式計算可得,水溶解度數據取自SRC公司的PHYSPROP數據庫。共提取特征參數M=31個。
對特征參數進行初篩選一般采用變量零值檢測、偏差測試、相關性測試、共線性測試。由于多元線性逐步回歸同時完成變量剔除和模型建立,本例采用此主體技術。
log Pow的QSPR模型:
(1)對7種結構性質描述符建模;
(2)對8種量子化學描述符建模;
(3)對14種拓撲描述符建模;
(4)對水溶解度、相對分子質量建模;
(5)對31個描述符綜合建模。
經多元線性逐步回歸,模型統計量匯總見表1。

表1 模型統計量匯總
由表1可見,綜合參數模型擬合性能最優:


模型中:Ovality——分子橢圓度,
PM-Y——慣性主矩的y分量,
Bindx——Balaban指數,
PSA——分子極化表面面積,
logSw——水溶解度對數值,
Mw——相對分子質量。
從數理統計角度,此構效關系模型擬合性能具有顯著的統計意義,尚應對模型預報性能進行評估驗證。
內部驗證。采用留一法預測標準偏差。

預測相關系數平方:

預測和擬合相關系數及標準偏差變化不顯著,模型具有穩定性。
外部驗證。對新考察的①萘、②苯甲酸甲酯、③苯甲酸乙酯、④苯乙酮、⑤二苯醚、⑥肉桂醇、⑦溴苯、⑧苯甲酸芐酯8種有機物,使用高效液相色譜實驗測定正辛醇/水分配系數,使用上述軟件計算6項參數。參數計算值、響應預測值、實驗值見表2。
外部驗證的標準偏差

驗證標準偏差相對穩定,具有統計意義。2004年QSAR國際會議正式形成經濟合作與發展組織(英文簡稱OECD)規則,明確必須使用外部驗證集(即測試集)來評價模型的預測能力。如果樣本量足夠大,也可以從105個樣本中隨機取8個樣本作為測試集,97個樣本作為訓練集。本案例執行該規范。
選自上海交通大學環境科學與工程學院的“Fenton法降解有機物的熱力學及構效關系研究”一文。2017年11月29日筆者在院內彭東輝教授級高工實驗室討論高濃度有機廢水處理項目時閱讀到此文,值得推薦。

表2 外部驗證數據匯總
案例針對有機廢水中的①3,4-二氯苯胺、②對氨基苯磺酸、③ 2,4-二硝基苯肼、④ 雙酚A、⑤ 酸性橙、⑥間甲酚紫6種有機污染物,實驗探究了不同溫度下Fenton氧化降解有機物的去除率及反應動力學;在Fenton試劑過量、假定反應初期為一級反應動力學速率常數的基礎上,通過Arrhenius方程計算獲得6種有機物的Fenton反應活化能Ea。
本文側重討論了結構參數與活化能之間的構效關系。根據現有資料,案例利用Hyperchem8.0軟件初步優化6個有機物的結構,再用Gaussian09軟件進行深度優化,利用Materials Studio軟件進行參數計算。選取了19個量子化學結構參數,其中16個由軟件計算獲得,3 個為組合參數(EGAP,E2GAP,EHOMO+ELOMO),各參數的物理意義見表3,各參數值及與Ea的相關性見表4。
由于論文報道的僅僅是研究工作的階段內容,樣本量為N=6,結構特征參數M=19,還無法構建構效關系模型,但可以進行初步分析:
(1)先計算單一特征參數和活化能的相關關系,見表4。

表3 量子化學結構參數
表4參數與Ea間相關系數大于0.5的項數11/19=58%,而且大于0.7的項數有2項,粗略判斷初選的結構參數是有效的。如果相關性普遍較低,就應考慮調整、擴充特征參數的選擇范圍。
(2)從案例一、二及相關資料分析可知,在構效關系研究中,由于化合物分子特征參數對化合物性能響應不同,很多參數對某一響應是不顯著的;更為普遍的現象是顯著項特征參數間還存在共線性現象,所以M項特征參數經篩選僅有限的m項進入構效關系模型。
(3)所謂共線性是指成對自變量(特征參數)之間的相關性。當相關性較高時,表示一個變量的信息含有對應變量的信息,可以剔除一個變量。如同時引入統計模型,除了增加計算工作量外,還會使模型計算性能變差。

表4 各參數值及與Ea的相關性
本講義第四講回歸分析中的變量篩選技術及統計檢驗(《上海化工》、2016年第8期)已做詳細討論。現再簡單引用,原案例是4個變量的多元線性回歸,計算機在逐步回歸建模時,自動計算輸出變量間的相關性,如表5所示。
x1與 x3相關系數 R=-0.824,P=0.001<0.05;x2與x4的 R=-0.973,P=0.000。可見:x1與 x3,x2與 x4均為顯著相關,建模時將作剔除處理。
(4)由于構效關系研究中樣本選擇是隨機的,多維空間結構參數的數值分布也是隨機的,構建構效關系模型一般要求N/m≥5,一定的樣本量保證模型的穩定性,所以完成本案例構效關系研究尚需一定的實驗工作量。

表5 四個變量間的相關性
后記:在本講義編寫中,賀少鵬提交了第一節的文稿并提供同濟大學的“QSAR模型內部和外部驗證方法綜述”等3篇資料。筆者2017年已欣然為賀少鵬寫推薦信,現其正在此領域選題攻讀在職博士。