CN102682209B - 有机污染物定量结构活性相关建模的变量筛选方法 - Google Patents

Abstract

本发明公开了一种有机污染物定量结构活性相关建模的变量筛选方法。首先计算全部的单变量和不同双变量组合的线性模型，二者均保留一定数量的最优模型。然后从保留的双变量线性模型中依次取出一个模型，将其中的两个变量与余下的变量中每一个变量分别组成三变量模型；直到所有保留的双变量模型处理完成，然后比较这些三变量模型的优劣，保留一定数量的最优三变量模型。依此类推，计算到构成模型的变量数满足要求时，停止计算。其中模型的优劣以Leave-One-Out Cross Validation(LOOCV)或者Leave-Multiple-Out Cross Validation(LMOCV)计算的q²或者均方根差RMSEV为终止标准。本发明原理简单，易于理解，容易编程实现，方法快速有效，保障了变量筛选的合理性以及模型的稳定性和预测能力。

Description

有机污染物定量结构活性相关建模的变量筛选方法

技术领域

本发明涉及一种有机污染物定量结构活性相关建模的变量筛选方法，具体地说是从大量的分子结构描述符变量中，选取一定数量相互作用较大的n-变量组合；然后以此为基础，每次增加一个变量，与所有选取出n-变量组合构成(n+1)-变量组合，并筛选出一定数量的(n+1)-变量组合，依次类推，直到满足要求为止，从而得到不同变量数的最优线性模型的一种变量筛选方法。

背景技术

有机污染物的定量结构活性相关(Quantitative Structure and ActivityRelationship，QSAR)研究方法作为一种计算机建模技术，能够深入挖掘有机污染物结构与其对人体及生态环境造成危害之间的量变规律和因果关系，为从分子水平上将污染物结构与其相对宏观的多样性环境行为和毒性效应进行关联提供了一种强有力的体外模拟技术。QSAR模型可以预测并弥补有机污染物环境行为和生态毒理数据的缺失、节省人力物力，有助于减少和替代实验(尤其是动物实验)测试，因此QSAR目前已经成为污染物环境生态风险评价和人体健康风险评价的重要工具(王连生，韩朔睽，有机污染物的定量结构-活性相关，1993，北京：中国环境科学出版社；国家自然科学基金委员会化学科学部组编，叶常明，王春霞，金龙珠主编，21世纪的环境化学，2004，科学出版社：北京。)

随着化学信息学的高速发展，QSAR建模研究的日益深入，表征分子结构的描述变量越来越多(目前可提取5000以上)。然而，因为结构表征描述符的提取方式、原理和途径等不一致，包含信息量的不一致，导致描述符变量之间存在大量的信息重叠或冗余，为QSAR模型的建立带来了困难。同时根据奥卡姆剃刀原则，要求模型尽可能简单，以避免复杂模型带来的过拟合问题。因此，筛选最佳的分子结构描述符变量子集，建立QSAR模型是成了有机污染物毒性效应建模与预测中急需解决的问题。

通常，解决大量描述变量建模方法还有主成分回归(Principle ComponentRegression，PCR)和偏最小二乘法(Partical Least Squares，PLS)，这两种方法都是从描述变量中分解出有原变量线性组合构成的几个重要的正交变量，然后建立模型的方法，可以解决部分线性相关及大变量集的问题(刘树深；易忠胜，基础化学计量学.科学出版社：北京，1999.王惠文；吴载斌；孟洁，偏最小二乘回归的线性与非线性方法.国防工业出版社：北京，2006.)。实际上，这两种方法所得模型质量不会超过通过适宜的变量筛选方法得到模型质量(Kubinyi，H.，VariableSelection in QSAR Studies.I.An Evolutionary Algorithm.QSAR Comb.Sci.1994，13(3)，285-294.)；其次，这两种方法，最后建立的模型是一个隐模型，也就是不能给出一个明确函数关系，虽然可以通过拟线性化得到线性模型(刘树深，有机物分子电性距离矢量表征及其应用.高等教育出版社：北京，2005.)，但变量数比较多的时候，非常不方便；因主成分为原变量的线性组合，不能直接给出分子结构描述符变量与有机污染物毒性效应之间的影响关系。因此，现在QSAR建模预测有机污染物毒性效应时，通常从众多的分子结构描述中筛选最佳子集，然后采用线性方法进行建模和预测。

目前QSAR建模时可用的变量筛选方法可以分作两类，一是基于全回归的方法，二是随机型的筛选方法。但不管是哪种方法都没有真正彻底地解决变量筛选的方法(Kapetanios，G.，Variable selection in regression models usingnonstandard optimisation of information criteria.Comput.Stat.Data Anal.2007，52(1)，4-15.González，M.P.；Terán，C.；Saíz-Urra，L.；Teijeir，M.，Variable SelectionMethods in QSAR：An Overview.Curr.Top.Med.Chem.2008，8(18)，1606-1627.)。如果变量的数量比较小，如变量数N＝30(所有的变量组合有2³⁰＝1073741824)，可以采用一些控制手段在短时间内通过比对全部的可能变量组合，获得最优模型，如VSMP方法(Liu，S.S.；Liu，H.L.；Yin，C.S.；Wang，L.S.，VSMP：A novelvariable selection and modeling method based on the prediction.J.Chem.Inf.Comput.Sci.2003，43(3)，964-969.)，并获得全局最优变量组合；当变量数比较多的时候，如变量数N＝100时，VSMP方法也很难到达要求，更不用说变量达到1000以上的数量了。其它的方法如前进法、后退法和逐步回归法，也是不宜处理大变量集(González，M.P.；Terán，C.；Saíz-Urra，L.；Teijeir，M.，Variable SelectionMethods in QSAR：An Overview.Curr.Top.Med.Chem.2008，8(18)，1606-1627.)。遗传算法(Kubinyi，H.，Variable Selection in QSAR Studies.I.An EvolutionaryAlgorithm.QSAR Comb.Sci.1994，13(3)，285-294.)、模拟退火(simulatedannealing，SA)、蚁群算法(ant colony algorithm，ACA)和微粒群(particle swarms，PS)等随机型的算法，在理论上可以解决大变量集筛选的问题，但实际操作过程中没有办法证明所得的结果为最佳(González，M.P.；Terán，C.；Saíz-Urra，L.；Teijeir，M.，Variable Selection Methods in QSAR：An Overview.Curr.Top.Med.Chem.2008，8(18)，1606-1627.)。因此变量筛选问题一直QSAR建模研究中的一个重要的方面。

通常分子结构描述符变量并不是独立地影响有机污染物毒性效应，同时各变量之间也不会相互独立、互不干扰，描述符变量之间总是存在相互作用的效应。如某个描述符变量与有机污染物毒性效应的相关系数非常小，几乎接近0，但是与另外的变量组合，可能是最佳的模型。以QSAR变量筛选的“标准”数据集——Selwood数据集(Selwood，D.L.；Livingstone，D.J.；Comley，J.C.W.；O′Dowd，A.B.；Hudson，A.T.；Jackson，P.；Jandu，K.S.；Rose，V.S.；Stables，J.N.，Structure-activity relationships of antifilarial antimycin analogs：a multivariatepattern recognition study.J.Med.Chem.1990，33(1)，136-142.)为例，第36、38、39和40个变量与因变量的相关系数几乎等于0，但是当与第50个变量组合时，却成了二变量模型中最好的四个。其原因正是分子结构描述变量并不是单独作用的，变量之间总会或多或少地存在相互作用(交互效应)，也就是说两个以上变量的组合会引入其它的信息。因此，根据变量的相互作用开发新变量筛选方法对采用QSAR技术对有机污染物毒性效应进行建模和预测有着重要的现实意义。

文献检索结果表明，采用基于变量相互作用进行有机污染物毒性效应与分子结构描述符建立QSAR模型的变量筛选方法未见报道。

发明内容

1、本发明要解决的技术问题

目前QSAR研究中通常采用的变量筛选方法通常存在两个基本的问题，一是全回归类型的变量筛选方法不能对大尺寸变量集进行有效快速地筛选，二是随机型的变量筛选方法无法验证是否最优及无法保证不同的筛选过程得到相同的结果。本发明的目的是针对以上技术难题，提供一种多元线性回归建模时，以描述变量的相互作用为基础，逐步从大变量集中筛选出变量数从1到指定数目变量的最佳变量子集，建立高稳定性和预测能力的多元线性模型，并科学地对有机污染物毒性效应进行预测。

2、技术方案

本发明的原理是将一定数量的最优n-变量子集系统地从余下变量的选取一个组成(n+1)-变量子集，并对新的变量子集进行线性模型质量评价，筛选出一定数量的最优(n+1)-变量子集，依次循环，直到变量数达到要求或者变量数大于样本数的5倍。一般采用留一法或者留多法交互验证的相关系数评价模型质量。

为便于基于变量相互作用的变量筛选方法(Variable Selection Method Basedon Variable Interaction，VSMVI)的技术方案的叙述，首先介绍多元线性回归(multiple linear regression，MLR)数学模型。

在许多实际问题中，和某随机变量y有关系的变量往往不是一个，而是多个，如有p个变量：x₁，x₂，…，x_n。研究因变量y与描述变量x₁，x₂，…，x_n之间的线性相关关系问题称之为多元线性回归问题。为寻找y与x₁，x₂，…，x_n之间的相关关系，假设收集到n组独立观测值：

(x_i1，x_i2，…，x_ip，y_i)，i＝1，2，…，n

并假定它们之间有如下关系式：

其中β₀，β₁，…，β_p是p+1个未知参数，x₁，x₂，…，x_p是描述变量，其值是可以精确测量或被控制的，y是可观测其值的因变量，ε是不可观测的随机误差。式(1)即为p元线性回归的数学模型。为方便起见，常采用矩阵表达式，并通过矩阵进行研究，记

$Y=\left(\begin{array}{c}{y}_1\\ y_2\\ \·\\ \·\\ \·\\ y_n\end{array}\right),$ $\mathrm{\β}=\left(\begin{array}{c}{\mathrm{\β}}_1\\ {\mathrm{\β}}_2\\ \·\\ \·\\ \·\\ {\mathrm{\β}}_n\end{array}\right),$ $X=\left(\begin{array}{cccc}1& x_{11}& \·\·\·& x_{1p}\\ 1& x_{21}& \·\·\·& x_{2p}\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ 1& x_{n1}& \·\·\·& x_{\mathrm{np}}\end{array}\right),$ $\mathrm{\ϵ}=\left(\begin{array}{c}{\mathrm{\ϵ}}_1\\ {\mathrm{\ϵ}}_2\\ \·\\ \·\\ \·\\ {\mathrm{\ϵ}}_n\end{array}\right)$

则式(1)可写成

$\left\{\begin{array}{c}Y=\mathrm{X\β}+\mathrm{\ϵ}\\ \mathrm{\ϵ}~N_n(0,{\mathrm{\σ}}^2{I}_n)\end{array}\right.$

数学上已经证明了如果漏掉了对y有显著影响的变量会使模型的估计与预测失去无偏性，若在回归方程中引入了对y无显著影响的变量后会降低估计和预测的精度。因此，当可用于建模的变量数f＞＞p时，需要通过变量筛选的方法，选择p个对y具有显著影响的变量进行回归建模。目前常用的变量筛选方法有后退法、前进法、逐步回归法、全回归法以及各种随机型变量筛选方法等。进行变量筛选时，需要一个筛选的终止标准。这些标准所采用统计参数大致有以下的几种：修正复相关系数R_adj、均方根差RMSE、预测偏差的方差JJ_p、平均预测均方差S_p，C_p统计量、预测平方和PRESS、AIC准则以及常用的交叉验证相关系数q²或者均方根差RMSEV等。这些统计参数的计算方法，可以参看相关的书籍与文献。本发明采用的交叉验证相关系数q²和RMSEV作为终止标准。

3、本发明采用的技术方案如下：

首先计算表征有机污染物分子结构特征的描述变量的单变量和不同双变量组合的线性模型，二者均保留一定数量的最优模型。然后从保留的双变量线性模型中依次取出一个模型，将其中的两个变量与余下的变量中每一个变量分别组成三变量模型；直到所有保留的双变量模型处理完成，然后比较这些三变量模型的优劣，保留一定数量的最优三变量模型。依此类推，计算到构成模型的变量数满足要求时，停止计算。其中模型的优劣以Leave-One-Out CrossValidation(LOOCV)或者Leave-Multiple-Out Cross Validation(LMOCV)计算的q²或者均方根差RMSEV为终止标准。LOOCV和LMOCV的计算方法可参见相关文献。

有机污染物定量结构活性相关建模的变量筛选方法具体步骤为：

(1)、单变量模型计算过程如下：

按照变量顺序，依次取一个变量，即模型变量数V_n＝1，建立一元线性模型，计算其相关系数r²，如果相关系数大于给定的临界值r_cri，则进行LOOCV或者LMOCV计算，否则选取下一个变量建立一元线性模型，并保留预先设定的Ns个最优的模型(包括模型的变量，r²、q²、RMSEE和RMSEV等统计参数)。直到所有变量计算完成。

(2)、双变量模型计算过程如下：

从所有两个变量的组合中依次选取一个，即模型变量数V_n＝2，建立二元线性模型，计算其复相关系数r²，如果相关系数大于给定的临界值r_cri，则进行LOOCV或者LMOCV计算，否则选取下一个双变量组合建立二元线性模型，并保留预先设定的Ns个最优的模型(包括模型的变量，r²、q²、RMSEE和RMSEV等统计参数)。直到所有变量计算完成。

(3)、三个及三个以上变量模型计算过程如下：

三个及三个以上变量模型的变量数V_n≥3，从保留的Ns个变量数为V_n-1的模型中选取一个模型，以该模型中的变量为基础，从余下的变量中每次选择一个变量与原来的V_n-1个变量构成变量数为V_n的变量组合，建立V_n元的线性模型，计算其复相关系数r²，如果相关系数大于给定的临界值r_cri，则进行LOOCV或者LMOCV计算，否则选取下一个变量组合，建立V_n元线性模型，直到所有变量计算完成；接着选取下一个变量数为V_n-1的模型，重复以上过程，直到保留的Ns变量数为V_n-1的模型全部处理完成。然后进行变量数为V_n+1的变量筛选，直到变量数满足要求为止。并保留预先设定的Ns个最优的模型(包括模型的变量，r²、q²、RMSEE和RMSEV等统计参数)。

3、有益效果

本发明充分利用了表征有机污染物分子结构特征的描述变量之间是相互作用的而不是独立影响因变量的特点，提供了一种快速有效的MLR建模过程中的变量筛选方法，将没有相互作用的单变量模型单独筛选，然后筛选出具有相互作用的双变量模型，然后以为基础，逐步增加变量，直到变量数满足要求，即可终止筛选计算。本方法是基于全回归方法，原理简单，没有涉及MLR计算以外的理论，易于理解，容易编程实现；方法快速有效，通过控制参数将质量不高的变量子集剔除，避免进行过多消耗机时的交互验证的计算，计算量与目前常用的方法相比大大地减少了；方法系统地对比指定变量数的所有子集，因而避免了变量子集的遗漏，保障了变量筛选的合理性。同时，以交互验证的相关系数为变量筛选终止标准和模型预测能力的判断指标，保证了模型的稳定性和预测能力。

附图说明

附图为本发明的整个计算流程图。

具体实施方式

以下通过实施例结合附图的计算流程进一步说明本发明。

实施例1：

选用变量筛选所谓“标准”测试集——Selwood数据集进行测试。该数据集最先出自文献(Selwood，D.L.；Livingstone，D.J.；Comley，J.C.W.；O′Dowd，A.B.；Hudson，A.T.；Jackson，P.；Jandu，K.S.；Rose，V.S.；Stables，J.N.，Structure-activityrelationships of antifilarial antimycin analogs：a multivariate pattern recognitionstudy.J.Med.Chem.1990，33(1)，136-142.)。该数据集包含31个样本，53个描述子。筛选过程中设置参数如下：保留模型数量Ns＝100、变量间相关系数r_int＝0.9、决定是否进行LOOCV或LMOCV计算的相关系数临界值的初始值r_cri＝0.1(但该值随着变量数的增加，也做相应的调节)，经过计算后得到如下表所示的结果，该数据集从未见文献报道变量数大于6的模型。

实施例2：

根据文献(易翔，郭宗儒，噻唑烷二酮和芳酮酸类PPAR-γ激动剂三维定量构效关系研究.药学学报2001，36(4)，262-268.)给出的58个PPAR-γ激动剂的结构和生物效应，采用虚拟计算化学实验室(VCCLAB)提供的E-Dragon软件，计算得到1664个分子结构描述子，经过预筛选后得到814个描述子。然后利用VSMVI方法筛选，筛选参数同实施例1。最终得到如下表所示的结果。

实施例3：

采用Dr.Igor V.Tetko提供的“Environmental Toxicity Prediction Challenge”训练集进行变量筛选测试。该训练集包括644个有机物，结构采用虚拟计算化学实验室(VCCLAB)的E-Dragon软件计算的1664个描述子表征，可以在http://www.cadaster.eu./node/65上得到。该数据及经过变量预筛选后得到827个描述子，VSMVI的参数同实施例1。最终得到如下的结果。

实施例4：

采用Dr.Igor V.Tetko提供的“Environmental Toxicity Prediction Challenge”训练集进行变量筛选测试。该训练集包括644个有机物，结构采用60个E-states指数表征，同样可以在http://www.cadaster.eu./node/65上得到。经过预筛选后得到50个描述子，VSMVI的参数同实施例1。最终得到如下的结果。

实施例5：

根据文献(Cronin，M.；Netzeva，T.；Dearden，J.；Edwards，R.；Worgan，A.，Assessment and modeling of the toxicity of organic chemicals to Chlorella vulgaris：development of a novel database.Chem.Res.Toxicol 2004，17(4)，545-554.)给出的91个有机物的结构和生物效应，采用虚拟计算化学实验室(VCCLAB)提供的E-Dragon软件，计算得到1664个分子结构描述子，经过预筛选后得到835个描述子。然后利用VSMVI方法筛选，筛选参数同实施例1。最终得到如下表所示的结果。

Claims (5)

1.一种有机污染物定量结构活性相关建模的变量筛选方法，其特征在于具体步骤为：

(1)、单变量模型计算过程如下：

按照变量顺序，依次取一个变量，即模型变量数V_n＝1，建立一元线性模型，计算其相关系数r²，如果相关系数大于给定的临界值r_cri，则进行LOOCV或者LMOCV计算，否则选取下一个变量建立一元线性模型，并保留预先设定的Ns个最优的模型，包括模型的变量，r²、q²、RMSEE和RMSEV统计参数；直到所有变量计算完成；

(2)、双变量模型计算过程如下：

从所有两个变量的组合中依次选取一个，即模型变量数V_n＝2，建立二元线性模型，计算其复相关系数r²，如果复相关系数大于给定的临界值r_cri，则进行LOOCV或者LMOCV计算，否则选取下一个双变量组合建立二元线性模型，并保留预先设定的Ns个最优的模型，包括模型的变量，r²、q²、RMSEE和RMSEV统计参数；直到所有变量计算完成；

(3)、三个及三个以上变量模型计算过程如下：

三个及三个以上变量模型的变量数V_n≥3，从保留的Ns个变量数为V_n-1的模型中选取一个模型，以该模型中的变量为基础，从余下的变量中每次选择一个变量与原来的V_n-1个变量构成变量数为V_n的变量组合，建立V_n元的线性模型，计算其复相关系数r²，如果复相关系数大于给定的临界值r_cri，则进行LOOCV或者LMOCV计算，否则选取下一个变量组合，建立V_n元线性模型，直到所有变量计算完成；接着选取下一个变量数为V_n-1的模型，重复以上过程，直到保留的Ns变量数为V_n-1的模型全部处理完成；然后进行变量数为V_n+1的变量筛选，直到变量数满足要求或者变量数大于样本的5倍为止；并保留预先设定的Ns个最优的模型，包括模型的变量，r²、q²、RMSEE和RMSEV统计参数。

2.根据权利要求1中所述有机污染物定量结构活性相关建模的变量筛选方法，其特征在于通过全回归的方法计算所有单变量即V_n＝1所构成一元线性模型，并根据模型的相关系数大小决定是否计算单变量模型的交叉验证统计参数，以加快计算速度，并保留Ns个最优的模型。

3.根据权利要求1中所述的有机污染物定量结构活性相关建模的变量筛选方法，其特征在于通过全回归的方法计算所有双变量组合即V_n＝2所构成的二元线性模型，并根据模型的复相关系数大小决定是否计算双变量模型的交叉验证统计参数，以加快计算速度，并保留Ns个最优的模型。

4.根据权利要求1中所述的有机污染物定量结构活性相关建模的变量筛选方法，其特征在于对三个及三个以上的变量组合即V_n>2，依次选取保留的Ns个V_n-1最优模型，每次增加一个变量，不包含每次选取的V_n-1最优模型中的变量，构成V_n个变量的组合，并根据模型的复相关系数大小决定是否计算V_n个变量模型的交叉验证统计参数，以加快计算速度，直到保留的Ns个V_n-1最优模型全部处理完成，并保留Ns个最优的模型，循环计算直到V_n达到预定的要求为止。

5.根据权利要求1～4中任一项所述的有机污染物定量结构活性相关建模的变量筛选方法，其特征在于将每一变量数下的最好的模型作为当前的最佳模型，然后通过不同变量数下最佳模型确定全局最优模型。