2014, Vol. 16
驱蚊剂可避免蚊虫袭扰人们日常生活，有效控制多种传染性疾病，开发高效新驱蚊剂意义重大[1]。美国于20世纪50年代研制合成了避蚊胺(diethyltoluamide) (N，N-二乙基间甲基苯甲酰胺)，具有广谱、高效和低毒等特点，随后被广泛生产和使用。此后，以酰胺类化合物为驱蚊剂的合成成为研究焦点[2]。
定量构效关系(Quantitative Structure-Activity Relationship，QSAR)研究对指导高效新驱蚊剂的分子设计与合成、阐明驱蚊剂的驱避机理有重要意义[3, 4, 5]。分子描述符的获取[6]、描述符的选择[7, 8]、回归模型的选择[9]及模型解释[10, 11, 12]是对已知生物活性的具有同类结构的化合物进行QSAR研究的4个关键步骤。Katritzky等[13]对DEET及其类似物共40个化合物的QSAR研究表明，沸点(蒸气压或挥发度)、分子的形状与亲脂性(疏水性)等是影响驱避活性的主要因素。Oliferenko等[14]从43个酰胺类似物出发，以埃及伊蚊Aedes aegypti的气味结合蛋白AaegOPB1为靶标，通过分子场拓扑分析和分子对接等筛选到了多个具有开发潜力的高驱避活性化合物。
笔者以40种酰胺类化合物为研究对象[15]，借助PCLIENT量子化学计算软件获取每个化合物的初始分子描述符1 773个，基于支持向量机(Support Vector Machine，SVM)，经二元矩阵重排过滤器(Binary Matrix Shuffling Filter，BMSF)与多轮末尾淘汰(Worst Descriptor Elimination Multi-round，WDEM)非线性筛选，获得8个保留分子描述符，建立了高精度的非线性SVR-QSAR模型，并进一步以SVR非线性解释体系分析了各保留描述符对化合物驱避活性的影响。
1 数据和方法 1.1 数据来源
40个酰胺类化合物的结构及其驱避活性引自文献[15, 16](见表 1)。活性指标为每个化合物对埃及伊蚊Aedes aegypti的有效保护时间。
 | 表 1 40个酰胺类化合物的保留描述符及生物活性 Table 1 Retained descriptors and biological activity of 40 amide compounds |
首先以软件PCLIENT(http://www.vcclab.org/lab/pclient/start.html)的JME编辑器画出每个化合物的分子结构式，然后将其导入任务窗口中，根据分子结构信息即可算得1 773个分子描述符。 1.3 分子描述符的非线性筛选
在所获取的1 773个分子描述符中存在大量冗余及无关的描述符，不仅影响建模精度且不利于模型解释。由于分子描述符与生物活性间常呈现为复杂的非线性关系，因此传统的逐步线性回归等线性特征筛选方法应用受限。前期，本研究室基于SVM建立了高维特征非线性选择新方法BMSF与低维特征非线性选择新方法WDEM，在基于芯片数据的癌信息基因选择和多肽定量序效建模等研究领域获得成功应用[7, 8]。
BMSF高维特征初筛算法如下：对于有n个样本，m个分子描述符的数据集定义为矩阵(yi，xij)，i=1…n，j=1…m。每个分子描述符有1(选取)和0(不选取)两种情况。产生一个包含0和1两个元素的随机矩阵K×m，规定每列0与1的个数相等。从随机矩阵中找出每行为1的元素所对应的原始训练集中的分子描述符，以SVR经10折交叉测试获得K个均方误差(Mean Square Error，MSE)值。K×m随机矩阵(自变量)与K个MSE (因变量)组成新训练集并训练建模，将随机矩阵的某列元素0、1互换后(其他列不变)作为测试集，预测得K个MSE0与K个MSE1，若均值MSE0≤均值MSE1，则剔除相应的分子描述符；遍历m次，得第一轮保留分子描述符。重复以上过程，经多轮筛选至没有分子描述符可被剔除为止[8]。
WDEM多轮末尾淘汰精筛算法如下：假设经BMSF筛选后还剩余m1个分子描述符，n个样本组成矩阵(yi，xij)，i=1…n，j=1…m1,使用这m1个分子描述符构建SVR回归模型，并以留一法做10折交叉测试计算出均方误差MSE0，然后剔除第j个分子描述符，用剩下的m1 -1个分子描述符再次建立回归模型，计算均方误差MSEj，比较MSEj与MSE0，若MSEj≤MSE0，则剔除分子描述符xj并进入下一轮筛选，反之筛选结束[7]。假定经筛选后剩余m2个分子描述符。 1.4 模型评估与解释
分子描述符选择完成后，需要选择一个具有足够解析能力的非线性回归模型，而常用的人工神经网络模型不适用于小样本，且易产生过拟合。本研究选用基于结构风险最小、非线性、适于小样本并能有效避免过拟合的支持向量回归(Support Vector Regression，SVR)模型[9]。SVR缺乏一个显性的表达式，可解释性差。前期，本研究室以F测验为基础，基于SVR建立了一套较完整的非线性解释性体系，包括模型回归显著性测验、单因子重要性显著性测验和单因子效应分析等，并验证了其合理性与有效性[10, 11, 12]。
1)模型回归显著性测验：以决定系数R2、平均相对误差(Mean Relative Error，MRE)、均方根误差(Root Mean Square Error，RMSE)及F检验评估模型优劣。
其中，n为样本个数，yi和
i分别为样本生物活性的实验值和预测值，y为所有样本生物活性实验值的均值。R2值愈大，MRE和RMSE值愈小，则模型回归性能愈优。
其中，U为回归平方和：
，Q为剩余离差平方和：
，m2为保留分子描述符数，模型显著性测验的自由度为(m2，n-m2-1)。若F >Fα(m2，n-m2-1)，则在α水平上SVR模型非线性回归显著[17]。
2)单因子重要性显著性测验：如果分子描述符xj(自变量)对驱蚊剂驱避活性y(因变量)有重要影响，则xj的变化将会引起驱避活性预测值的明显变化。由此推知，若分子描述符固定为其零水平表达值xj，代入回归模型，得到相应预测活性值后，基于回归平方和Uj、剩余平方和Qj，可计算获得分子描述符xj对回归平方和增加所做的贡献值(U-Uj)。对多元线性回归模型，离差平方和 
；而对于SVR模型，离差平方和则为SSy≠Q+U。为比较各保留描述符之间的重要性，引入公式Qj′=Qj/(Qj+Uj)×SSy和Uj′=Uj/(Qj+Uj)×SSy，对Qj、Uj进行规格化处理后可得：SSy=Qj′+Uj′。规定Vj=U′-Uj′，得到统计量Fj。若Fj>Fα(1，n-m2-1)，则在α水平上单因子重要性显著。
3)单因子效应分析：在驱蚊剂定量构效关系研究中，若能明确各保留分子描述符变化对驱蚊剂驱避活性的影响趋势，将对设计与合成高效新驱蚊剂具有指导意义。对单个保留分子描述符xj进行效应分析，可将除分子描述符xj外的其他保留分子描述符都固定为其平均值，使指定分子描述符xj在已规定的取值范围内按照一定步长取值，回代入SVR模型后得到预测驱避活性值j，并绘制相应的xj-j趋势图[17]。由此可知当其他保留分子描述符取其平均值时，驱避活性随特定描述符xj的变化规律。
本研究中BMSF高维特征初筛、WDEM多轮末尾淘汰精筛、SVR建模和非线性解释体系等均采用自编MATLAB程序通过调用LIBSVM3.1软件包实现[18]。核函数选用径向基核，核函数参数采用该软件包中的gridregression.py程序默认范围、步长经格点搜索自动获取。
2 结果与分析 2.1 供试化合物的SVR-QSAR模型
基于40个样本，1 773个初始分子描述符的SVR模型F=2.73；经BMSF非线性高维特征初筛后，16个描述符的SVR模型F=210.12；再经WDEM精筛后8个保留描述符的SVR模型F=8 465>F0.01(8,31)，R2=0.999 6，非线性回归达极显著。可见特征筛选效果明显。
Bhonsle等[15]基于这40个样本，30个分子描述符构建的多元线性回归方程R2=0.989 4。显然，本研究以较少的分子描述符构建的SVR模型精度更高。
基于8个保留分子描述符的SVR模型，进一步以留一法(Leave One Out，LOO)交叉验证给出各样本的预测值(表 1，图 1)。由图 1可见，预测值与实验值分布于坐标轴对角线附近，R2LOO=0.935 1，MRE=37.87%，RMSE=0.382，进一步表明本研究所建SVR模型稳定性好，可信度高。
 | 图 1 留一法检验的驱避剂生物活性实验值与预测值 Fig. 1 Experimental values and predicted values of repellent bio-activities with Leave-One-Out test |
基于SVR可解释性体系，对8个保留分子描述
符的单因子重要性显著性进行分析。结果表明，其F值均大于临界值F0.01/8(1,31)=7.53，达极显著(表 2)。8个保留分子描述符对化合物驱避活性影响的重要性依次为TPSA(Tot)> RDF035e>BLTF96>Eig1e>RDF055u>IC4>RDF075e>G(N..O)。
| 表 2 特征筛选后的8个保留分子描述符 Table 2 Eight retained descriptors after feature screening |
其单因子效应分析结果见图 2。在40个样本保留分子描述符取值范围内，当TPSA(Tot)＝20.31、RDF035e＝18.61、BLTF96＝-3.21、Eig1e＝25.01、RDF055u＝21.31、IC4＝4.68、RDF075e＝15.07和G(N..O)＝0时，化合物具有最高的预测驱避活性；保留分子描述符与驱避活性呈非线性关系，其中RDF035e和BLTF96与驱避活性呈开口向下抛物线变化，Eig1e与驱避活性呈开口向上抛物线变化，IC4、RDF055u和RDF075e与驱避活性呈正相关，G(N..O)和TPSA(Tot)与驱避活性呈负相关。
 | 图 2 8个保留分子描述符的单因子效应 Fig. 2 Single-factor effects of 8 retained descriptors |
对于一个好的QSAR模型，首先要确保所获得的分子描述符能全面表征化合物结构与活性间的复杂关系。通过量子化学软件，对每个化合物简捷地获取尽可能全面的、数以千计的初始分子描述符是较为理想的选择。然而，高维特征、小样本不仅导致“维数灾难”，且无关与冗余描述符还会影响建模精度，增加模型复杂性并使得模型解释困难，(非线性)特征选择此时变得尤为关键。本研究中，1 773个初始分子描述符的SVR模型F=2.73，8个保留分子描述符的SVR模型F=8 465，再次证实本研究室前期发展的高维特征非线性选择新方法BMSF与低维特征非线性选择新方法WDEM是有效的。分子描述符与生物活性间往往存在复杂的非线性关系。与Bhonsle等[15]报道的含30个分子描述符的多元线性回归模型(R2=0.989 4)相比，本研究所建的8个保留分子描述符的非线性SVR模型精度更高(R2=0.999 6)；图 2也显示多个保留分子描述符与驱避活性的单因子效应呈抛物线变化。因此，在QSAR研究中，应优先选用基于结构风险最小、非线性、适于小样本、能有效避免过拟合的SVR为基本建模工具。
本研究所获的8个保留分子描述符分别属于化合物分子特性、径向分布描述子、边缘邻接指数、信息指数和三维原子对5组，表明驱避活性与化合物多种性质间存在复杂关联。其中，拓扑极性分子表面积TPSA(Tot)最为重要，其值越小，活性越高，支持文献报道的“良好的驱避剂化合物需要合适的分子表面积来跟气味感受器发生交互作用，达到驱避效果，相对偏小的分子表面积有利于跟气味受体作用”的推论[5, 19]。RDF035e、RDF075e和RDF055u均属于表示电子相关性的电负性径向分布函数，在本研究中作为保留描述符得到Katritzky等[13]认为的“酰胺类化合物的电荷分布对其驱避活性有很大影响”的支撑。单因子效应分析进一步显示，RDF055u和RDF075e与驱避活性呈正相关(图 2)，与王宗德等[5]“负电性有助于蚊虫触角上嗅觉感受器感知气味而起到驱避作用，即负电性越强，驱避活性越高”的结论一致。佐证了本研究所得的保留分子描述符与单因子效应分析方法的合理性。
本研究建立了8个保留分子描述符及40个样本的SVR-QSAR模型，为设计新的、高效的酰胺类驱蚊剂分子奠定了基础。将来拟进一步虚拟构建多个酰胺类化合物，经PCLIENT量子化学计算在线获取虚拟化合物的8个保留分子描述符，代入模型预测，取预测活性最高且大于6 h(原数据集7c化合物的保护时间最高为 6 h)的少量虚拟化合物，通过化学合成及生物试验验证其驱蚊活性。
| [1] | 廖圣良, 姜志宽, 宋杰, 等. 蚊虫驱避剂的QSAR研究[J]. 中华卫生杀虫药械, 2011, 17(6): 469-471. LIAO Shengliang, JIANG Zhikuan, SONG Jie, et al. The study of mosquito repellent QSAR[J]. Chinese J Hygienic Insecticides Equip, 2011, 17(6): 469-471.(in Chinese) |
| [2] | 陶波, 张大伟. 蚊虫驱避剂的研究进展[J]. 东北农业大学学报, 2014, 45(2): 123-128. TAO Bo, ZHANG Dawei. Research advance of mosquito repellents[J]. J Northeast Agric Univ, 2014, 45(2): 123-128. (in Chinese) |
| [3] | 韩晓峰, 刘莹, 高莹, 等. 非肽类凝血酶抑制剂的比较分子力场分析[J]. 化学学报, 2003, 61(7): 1136-1139. HAN Xiaofeng, LIU Ying, GAO Ying, et al. Comparative molecular field analysis of non-peptidic inhibitors of thrombin[J]. Acta Chimica Sinica, 2003, 61(7): 1136-1139. (in Chinese) |
| [4] | 梅虎, 周原, 孙立力, 等. 一种新的氨基酸描述子及其在肽QSAR中的应用[J]. 物理化学学报, 2004, 20(8): 821-825. MEI Hu, ZHOU Yuan, SUN Lili, et al. A new descriptor of amino acids and its application in peptide QSAR[J].Acta Phys Chim Sin, 2004, 20(8): 821-825. (in Chinese) |
| [5] | 王宗德, 宋杰, 姜志宽, 等. 驱避剂的构效关系和驱避机理的研究[J]. 中华卫生杀虫药械, 2008, 14(6): 472-476. WANG Zongde, SONG Jie, JIANG Zhikuan, et al. Study of the structure-activity relationship and mechanism of repellent[J]. Chin J Hygienic Insecticides Equip, 2008, 14(6): 472-476. (in Chinese) |
| [6] | NATATAJAN R, BASAK S C, MILLS D, et al. Quantitative structure-activity relationship modeling of mosquito repellents using calculated descriptors[J]. Croatica Chemica Acta, 2008, 81(2): 333-340. |
| [7] | 代志军, 周玮, 袁哲明. 基于支持向量机的高维特征非线性快速筛选与肽QSAR建模[J]. 物理化学学报, 2011, 27(7): 1654-1660. DAI Zhijun, ZHOU Wei, YUAN Zheming. A novel method of nonlinear rapid feature selection for high dimensional data and its application in peptide QSAR modeling based on Support Vector Machine[J]. Acta Phys Chim Sin, 2011, 27(7): 1654-1660. (in Chinese) |
| [8] | ZHANG Hongyan, WANG Haiyan, DAI Zhijun, et al. Improving accuracy for cancer classification with a new algorithm for genes selection[J]. BMC Bioinformatics, 2012, 13(1): 298. |
| [9] | VAPNIK V N. The Nature of Statistical Learning Theory[M]. New York: Springer Verlag Press, 1995: 87-189. |
| [10] | 李俊, 谭显胜, 谭泗桥, 等. 改进支持向量机在棉铃虫人工饲料配方优化中的应用[J]. 昆虫学报, 2010, 53(4): 420-426. LI Jun, TAN Xiansheng, TAN Siqiao, et al. Application of improved Support Vector Machine in the optimization of artificial diet for the cotton bollworm, Helicoverpa armigera (Lepidoptera Noctuidae)[J]. Acta Entomologica Sinica, 2010, 53(4): 420-426. (in Chinese) |
| [11] | 周世豪, 李俊, 姚润贤, 等. 基于均匀设计与支持向量回归的棉铃虫幼虫全纯人工饲料配方优化[J]. 昆虫学报, 2012, 55(1): 124-132. ZHOU Shihao, LI Jun, YAO Runxian, et al. Optimization of chemically defined diet for larvae of the cotton bollworm(Helicoverpa armigera) based on uniform design and Support Vector Regression[J]. Acta Entomologica Sinica, 2012, 55(1): 124-132. (in Chinese) |
| [12] | 戴长庚, 李凯龙, 王立峰, 等. 基于均匀设计优化的大螟实用饲料配方及继代饲养[J]. 中国水稻科学, 2013, 27(4): 434-439. DAI Changgeng, LI Kailong, WANG Lifeng, et al. An oligidic diet for Sesamia inferens optimized by uniform design and successive rearing[J]. Chin J Rice Sci, 2013, 27(4): 434-439. (in Chinese) |
| [13] | KATRITZK A R, DOBCHEV D A, TULP I, et al. QSAR study of mosquito repellents using Codessa Pro[J]. Bioorg Med Chem Lett, 2006, 16(8): 2306-2311. |
| [14] | OLIFERENKO P V, OLIFERENKO A A, PODA G I, et al. Promising Aedes aegypti repellent chemotypes identified through integrated QSAR, virtual screening, synthesis, and bioassay[J]. PLoS One, 2013, 8(9): e64547. |
| [15] | BHONSLE J B, BHATTACHARJEE A L, GUPTA R K. Novel semi-automated methodology for developing highly predictive QSAR models: application for development of QSAR models for insect repellent amides [J]. J Mol Model, 2007, 13(1): 179-208. |
| [16] | SURYANARAYANA M V S, PANDEY K S, PRAKASH S, et al. Structure-activity relationship studies with mosquito repellent amides[J]. J Pharm Sci, 1991, 80(11): 1055-1057. |
| [17] | 谭显胜, 王志明, 谭泗桥, 等. 支持向量回归可解释性体系的建立[J]. 系统仿真学报, 2009, 21(24): 7795-7797. TAN Xiansheng, WANG Zhiming, TAN Siqiao, et al. Establishing interpretability for Support Vector Regression[J]. Journal of System Simulation, 2009, 21(24): 7795-7797. (in Chinese) |
| [18] | CHANG Chih-Chung, LIN Chih-Jen. LIBSVM: a library for Support Vector Machines[J]. ACM Trans Intell Syst Technol, 2011, 2(3): Article No. 27. |
| [19] | 廖圣良, 姜志宽, 宋杰, 等. 蚊虫驱避剂的驱避机理研究[J]. 中华卫生杀虫药械, 2012, 18(4): 280-283. LIAO Shengliang, JIANG Zhikuan, SONG Jie, et al. Repelling mechanism of mosquitoes repellent[J]. Chin J Hygienic Insecticides Equip, 2012, 18(4): 280-283.(in Chinese) |