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

摘要

本发明公开了一种有机污染物定量结构活性相关建模的变量筛选方法。首先计算全部的单变量和不同双变量组合的线性模型，二者均保留一定数量的最优模型。然后从保留的双变量线性模型中依次取出一个模型，将其中的两个变量与余下的变量中每一个变量分别组成三变量模型；直到所有保留的双变量模型处理完成，然后比较这些三变量模型的优劣，保留一定数量的最优三变量模型。依此类推，计算到构成模型的变量数满足要求时，停止计算。其中模型的优劣以Leave-One-Out Cross Validation(LOOCV)或者Leave-Multiple-Out Cross Validation(LMOCV)计算的q

说明书

技术领域

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

背景技术

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

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

通常，解决大量描述变量建模方法还有主成分回归(Principle Component  Regression，PCR)和偏最小二乘法(Partical Least Squares，PLS)，这两种方法都是 从描述变量中分解出有原变量线性组合构成的几个重要的正交变量，然后建立 模型的方法，可以解决部分线性相关及大变量集的问题(刘树深；易忠胜，基础化 学计量学.科学出版社：北京，1999.王惠文；吴载斌；孟洁，偏最小二乘回归的 线性与非线性方法.国防工业出版社：北京，2006.)。实际上，这两种方法所得模 型质量不会超过通过适宜的变量筛选方法得到模型质量(Kubinyi，H.，Variable  Selection 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 using  nonstandard 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 Selection  Methods 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 novel  variable 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 Selection  Methods in QSAR：An Overview.Curr.Top.Med.Chem.2008，8(18)，1606-1627.)。 遗传算法(Kubinyi，H.，Variable Selection in QSAR Studies.I.An Evolutionary  Algorithm.QSAR Comb.Sci.1994，13(3)，285-294.)、模拟退火(simulated  annealing，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 multivariate  pattern 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 Based  on 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),$$\beta =\left(\begin{array}{c}{\beta }_1\\ {\beta }_2\\ ·\\ ·\\ ·\\ {\beta }_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),$$ϵ=\left(\begin{array}{c}{ϵ}_1\\ {ϵ}_2\\ ·\\ ·\\ ·\\ {ϵ}_n\end{array}\right)$

则式(1)可写成

$\left(\begin{array}{c}Y=\mathrm{X\beta }+ϵ\\ ϵ~N_n(0,{\sigma }^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 Cross  Validation(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-activity  relationships of antifilarial antimycin analogs：a multivariate pattern recognition  study.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。最终得到如下表所示 的结果。