CN101131391B

CN101131391B - 基于分子亲电矢量和扩展支持向量机的基因毒性概率预测方法

Info

Publication number: CN101131391B
Application number: CN2006100303584A
Authority: CN
Inventors: 蒋华良; 罗小民; 朱维良; 陈凯先; 郑明月; 刘治国; 薛春霞
Original assignee: Shanghai Institute of Materia Medica of CAS
Current assignee: Shanghai Institute of Materia Medica of CAS
Priority date: 2006-08-24
Filing date: 2006-08-24
Publication date: 2011-07-20
Anticipated expiration: 2026-08-24
Also published as: CN101131391A

Abstract

本发明涉及一种基于分子亲电矢量(MEV)和扩展支持向量机(SVM)的基因毒性概率预测方法，适用于根据有机化合物分子结构信息对该化合物进行虚拟毒性评价和筛选。首先，对输入分子结构使用基于文本的化学结构提问语言SMARTS和可编程原子分类PATTY回溯算法根据预定义规则进行原子分类；然后，根据半经验Hückel方法计算分子中每种原子类型的原子描述符(前线轨道电子密度、电子超离域度、原子π-电荷)建立描述分子的亲电性质的特征矢量MEV；最后，使用扩展的支持向量机方法(SVM)对分子的基因毒性数据和及其特征MEV进行统计学习，得到该分子基因毒性的后验概率估计。

Description

基于分子亲电矢量和扩展支持向量机的基因毒性概率预测方法

技术领域

本发明涉及一种基于分子亲电矢量(MEV)和扩展支持向量机(SVM)的基因毒性概率预测方法，适用于根据有机化合物分子结构信息对该化合物进行虚拟毒性评价和筛选。

背景技术

药物开发依赖于发现具有目标活性，同时毒副作用低的化合物。在过去的几十年中，组合化学和高通量筛选(HTS)等药物发现技术在早期确定先导化合物方面已经取得了实质性进展。然而，毒性问题仍是后期药物失败的一个重要的因素(Caldwell，et al.，2001)。当前，为了评价药物的安全性，需要进行一系列毒性测试试验，除了巨大的经济和人力开销，这类测试普遍受到低通量筛选性能的限制。因此，既不能用于早期发现阶段的药物毒性评价，也不能在大规模临床应用之前检测药物的毒副作用(Johnson and Wolfgang，2000)。近来，作为体内、体外试验有力补充的虚拟毒性预测模型开发领域开始引起越来越多的关注。长期的毒性试验已经积累了大量的结构活性关系(SAR)信息，为建立基于分子结构和反应性的计算模型奠定了基础；众多能够支持产生分子碎片模式和计算拓扑、物化、电子描述符等的软件以及统计学工具，也使得预测毒性的模建成为可能(Benfenati and Gini，1997；Benigni，2005；Fielden，et al.，2002；Greene，2002；Helma，2005；Johnson and Wolfgang，2000)。

基因毒性是指化合物造成DNA突变的能力，是一类与药物开发紧密关联的毒性指标。一方面，标准的基因毒性检测试验-Ames试验，目前已经成为药物审批必需的一个环节。另一方面，越来越多试验结果表明，基因毒性与致癌性高度相关(77％-90％)，可以说Ames试验是所有体外生物实验中对致癌性预测最准确的方法。因此，成功的虚拟基因毒性预测模型甚至可以作为啮齿类动物致癌性的评价工具(Kimand Margolin，1999；Zeiger，et al.，1990)。此外，开发准确的毒性预测模型要求数据集化合物的毒性来源主要集中于一种机制。从这个角度来讲，因为其相对单一的发生机理，基因毒性也应该比其他类型的毒性更容易预测(Snyder and Smith，2005)。

到目前为止，已经开发了一些计算工具可以对基因毒性进行评价，主要可以分为基于知识的和统计学驱动的模型。前者可以进一步划分为基于规则的和基于碎片的模型。基于规则的程序如DEREK(deductiveestimate of risk from existing knowledge)，可以根据预定义规则(专家从大量试验数据中总结出来的)对某一毒性反应是否发生给出标记值。基于碎片的方法以MCASE(Multiple Computer AutomatedStructure Evaluation)为代表，将每个输入分子拆分成为原子和基团碎片，然后从统计学角度评价这些碎片和数据库中相似碎片(基因毒性片断)之间的关联度，并给出相关性打分。统计学驱动的程序TOPKAT(Toxicity Prediction by Komputer Assisted Technology)使用电子和拓扑描述符，根据化合物结构类别得到单独的预测模型，是一种经典QSAR分析方法的扩展。对这些计算系统更详细的评价可见参考文献(Greene，2002)。尽管乐观的观点认为这些工具也许在将来可能会减轻或取代试验毒性测试的需要，但事实上与这一目标还有相当的距离。近来的研究表明，上述的所有商业程序对基因毒性的预测能力都十分有限(表1)，尤其是以整体敏感度(真阳性率)作为性能衡量指标时，仅对43～52％的真阳性化合物可以给出正确的预测(Snyder and Smith，2005)。在对一些专利药品的毒性预测中，White等人也得到类似的结果，说明这些模型在预测性能方面，尤其是基因毒性预测的敏感度方面亟待提高(White，et al.，2003)。

表1 现有商业程序对已上市药物的细菌基因毒性评价中的预测性能比较。

	敏感度	特异性	整体正确率
				MCASE	13/27(48％)	307/330(93％)	320/357(90％)
DEREK	14/27(52％)	260/346(75％)	274/373(74％)
				topKAT	10/23(43％)	267/316(85％)	277/339(82％)

除上述商业软件包之外，文献中也报道了一些基因毒性预测的方法。总体来说基于知识的方法更倾向于提供毒性发生的机制信息，但这些预先定义的碎片或规则是现有知识的归纳表达，而非新知识的演绎推理。统计驱动的方法正相反，具有对现有知识进行扩展的能力，不会对某一种毒性作用机制的预测形成偏倚；但这些模型的性能通常会受到诸如分子描述符的质量，训练集和测试集的分子结构多样性，以及统计学习方法的有效性等因素的影响。最近，Helma等人开发了一套分子特征发掘算法(MOLFEA)用于基因毒性的预测，结果表明基于碎片的分子特征描述符建立的模型比基于整体性质的描述符的模型可以给出更为精确的预测(10～15％以上)(Helma，et al.，2004)。该研究还比较了不同的统计学习方法的建模能力，支持向量机(SupportVector Machine)给出最高的预测精度：对10倍交叉验证(CrossValidation)整体正确率为78％。

除了性能方面的局限性，当前的虚拟毒性预测方法很少涉及毒性概率预测，如基因毒性的预测通常被认为是一个二元分类问题。但实际应用中简单的化合物是否有毒的表示方法通常是不能满足需要的。比如在化合物定购或已知有机合成过程中，一个错误的阳性预测通常意味着淘汰掉一个有前景的药物候选物。概率预测不仅能帮助研究人员指定更为灵活的取舍标准，而且可为他们的决定提供定量的可靠性评估。

发明内容

本发明的目的主要是提供一种判别能力强，同时又提供相关概率信息的基因毒性评价方法。

该目的、主要是采用分子表征技术和统计学习方法，通过对现有分子结构和相关基因毒性数据进行训练而达到的：

1)结构活性(毒性)关系模型的建立要求选择与感兴趣的性质最相关的结构或化学特征。对于基因毒性，描述符的选择应该充分考虑到可能引发DNA突变的机制。小分子基因毒性的产生机制主要包括如下几个方面：

a.形成DNA加合物或者造成碱基缺失；

b.非反应性分子代谢后形成活性中间体、与DNA反应；

c.通过非共价相互作用插入DNA双螺旋，造成DNA空间结构扭曲。

这些效应干扰正常的DNA修复和复制，从而导致突变，最终造成基因毒性。当前的研究表明典型的基因毒性主要来自于由分子亲电性引起的直接的化合物-DNA相互作用(Snyder and Smith，2005)。David等人的QSAR研究(David F.V.Lewis，2003)显示，化合物的亲电性对间接致突变活性(机制b)的产生也有重要作用。基于上述原因，我们设计了一套基于原子指标的分子亲电矢量(MEV)来描述化合物的亲电性特征。

2)成功的结构活性模型还依赖于采用的统计学方法的有效性。在这里我们选择具有强泛化能力支持向量机(SVM)方法。近期，大量的模式识别研究证明支持向量机在不同程度上可以比其他监督的统计学习方法给出更为精确的预测(Bock and Gough，2001；Doniger，et al.，2002；Li，et al.，2005；Lo，et al.，2005)；另外，Wu等人提出一种扩展的支持向量机方法可以用于分类的概率性预测(Wu and Lin，2004)。在本发明中我们采用Chang等人实现的这种扩展支持向量机(Chang and Lin，2001)，建立最后的毒性预测模型。

3)训练集来自于文献(Kazius，et al.，2005)，其中包括4337化合物结构和药物审批所需要的指示是否具有对Salmonella Typhimurium菌的致突变毒性信息。我们使用同文献(Kazius，et al.，2005)相同的分类方法，得到的数据集中54％是毒性性化合物(2401个致突变化合物，1936个非致突变化合物)。验证数据集从各种公共来源得到，如美国环保局EPA(http://www.epa.gov)和国立卫生研究院NIH(http://www.nih.gov)。去除与训练集相同的化合物，混合物和共振异构结构，共得到592个数据样本，其中54.90％具有致突变性。数据集中的化合物涵盖广泛的化学多样性空间，为建立健壮、有效的预测模型奠定了基础。

表2中列出了本发明的SVM+MEV模型的整体性能参数。由于使用了相同训练集数据，我们还同时给出了Kazius等人的毒效团模型(Kazius，et al.，2005)的统计结果以作比较。对于训练集，我们的SVM+MEV方法显示了卓越的数据拟合能力，整体预测正确率达到了91.86％。敏感度和专一性、分别为93.63％和89.67％，超过Kazius等人的毒效团模型10％以上。对于外部测试集，我们的模型可以给出84.80％的预测正确率，接近Salmonella试验测定基因毒性的方法重现性(～85％)(Benigni and Giuliani，1988)。如Kazius等人指出，鉴于当前的试验条件和结构活性关系方法的内在局限性，虚拟毒性预测模型理论上不可能给出超过85％整体正确率(Kazius，et al.，2005)。在本次研究中，我们还使用Accelry公司的TOPKAT对外部测试集进行评价。TOPKAT中的基因毒性预测模块是基于“电子-拓扑”描述符建立的，与我们的MEV模型同属于统计学驱动模型，这种对于同一数据集的对照测试有助于比较我们的模型与商业程序的性能差别。如表3所示的统计学结果，TOPKAT程序具有正确识别阴性样本的能力，专一性(真阴性率)可以达到85.10％；然而，该程序对阳性样本并没有同样有效的预测能力，敏感度(真阳性率)和整体预测正确率仅分别为77.32％和80.81％，显著的低于我们的MEV+SVM模型。其他著名商业程序诸如MCASE和DEREK也存在这种对阳性化合物预测能力差的问题(较差的敏感度)在中也普遍存在(Snyder and Smith，2005)。我们的SVM+MEV模型的预测敏感度和专一性都具有比较高的水平，显示了该模型在识别致突变和非致突变化合物方面具有同样非常平衡的能力。

表2.MEV+SVM模型与Kazius’(2005)毒效团模型和商业程序TOPKAT预测性能指标的比较。

作为进一步的分析，我们还对训练集数据进行了交叉验证，结果如表3所。同样是采用10倍交叉验证，与Helma等人提出的MOLFEA模型(Helma，et al.，2004)和Kazius等人的完备化学表征方法(Kazius，etal.，2006)相比，我们的SVM+MEV模型的预测精度高出10～12％。由于SVM+MEV模型可以得到足够高的交叉验证正确率，而且在对外部测试集的预测中基本保持这样的性能，因此，可以认为该方法能够有效的从测试集数据中抽取基因毒性相关的信息。

表3.MEV+SVM模型与其他分子结构信息挖掘方法的性能(10倍交叉验证)比较

模型	样本数	敏感度(％)	特异性(％)	整体正确率(％)	相关系数
						MEV+SVM	4337	92.09	87.71	90.13	0.8003
MOLFEA(Helma，et al.，2004)	684	77.50	79.40	78.51	0.5692
						ECR(Kazius，et al.，2006)	4069	83.00	74.04	80.63	0.5735

除了二元分类之外，实际应用中人们通常对分类的预测可信度更感兴趣；而实际上，试验毒性评价也应该是一个宏观和统计的结果。例如，一个化合物的预测毒性概率(P)为90％的化合物应该比P为60％的化合物更可能具有基因毒性，这对虚拟毒性评价模型来说应该是更合理的输出结果。图1显示了模型对所有训练集样本的整体正确率(Q)随预测毒性概率(P)之间的变化关系：当P达到两个端点0％和100％时，模型的预测性能达到最大值，当P接近中值50％时，预测正确率显著下降。因为毒性概率P在0％附近时实际意味着非毒性概率接近100％，这一点在某种程度上容易造成混乱。为清楚起见、本发明使用P’同时代表致毒性和非毒性可能性，其值可用一分段函数表示为：

\{\begin{matrix} P^{'} = P & (P &GreaterEqual; 50 %) \\ P^{'} = 100 - P & (P < 50 %) \end{matrix}.

经过这一转换，我们可以清楚的发现P′和预测精度高度相关：P′值越高，模型的预测结果越准确，该化合物引发毒性的倾向也就越确定。对于外部测试集，预测P′值高于75％的化合物(405 compounds of test set)可以以高达91.85％的Q值进行估计，显著的(也是局部的)提高了预测性能。尽管平均预测精度没有改变，我们仍然可以通过对化合物的P′值对应的性能水平制定更为合理的决策。

本发明中发明者基于一种新型的分子表征方法MEV和扩展的SVM，发展了一套用于基因毒性概率预测的虚拟药物分子毒性评价模型。一方面，在MEV的设计中我们吸收和整合现有知识，充分考虑了与基因毒性引发机制高度相关的分子亲电和拓扑性质；另一方面，我们利用SVM强大的统计学习能力扩展现有知识，避免预测结果对某种特定毒性机制发生偏倚。我们的模型对大规模的数据集也可保持良好的预测性能，说明这种MEV和SVM组合的方法对药物分子基因毒性预测模建是合理和有效的。根据最终的MEV+SVM+RFE模型开发的C++应用程序ZGTOX，可以快速(～200个化合物/秒，P4 2.4G)、准确的对化合物的毒性概率进行预测，在高通量虚拟筛选和虚拟毒性评价中有良好应用前景。

附图说明

图1为整体预测正确率随预测基因毒性概率(P，％)的变化关系(左下坐标系)，取样数(N)对P的变化关系(右下坐标系)。

图2为以苯胺为例，将化合物结构信息转化为该化合物的特征分子亲电性矢量。

具体实施方式

基于MEV和扩展SVM的分子基因毒性评价模型的建立主要涉及五个步骤：

1)对数据集中的化合物进行原子分类：

鉴于实现的可移植性和便捷性，我们使用基于文本的化学结构提问语言SMARTS(SMiles ARbitrary Target Specification)描述所有原子分类(表4)。对每个原子的类型由其自身的化学性质和反应其化学环境的邻接原子和键类型所决定。然后我们使用OpenBabel(http://openbabel.sourceforge.net)C++库中的可编程原子分类PATTY(Programmable atom typer)回溯算法(Bush and Sheridan，1993)完成原子类型指定。通过使用SMARTS和PATTY，我们可方便的从化学和毒理学角度出发对原子进行灵活有效的分类。

表4.52个原子分类定义规则和3个毒效团子结构定义，使用SMARTS语法描述。

类型

编号

SMARTS字串

描述

C.3G	1	[C＾3]	generic sp3 carbon
				C.3X	2	[$([C＾3][Cl，Br，I，$(C#N)，$(N(～O)～O)，$(OS(＝O)＝O)])，$(ON＝O)]	sp3 carbon connected to strong electronwithdrawing group
C.3H	3	[$([C＾3]([#6＾3，#1])([#6＾3，#1])([#6＾3，#1])[#6＾3，#1])]	sp3 carbon connected to sp3 carbon orhydrogen
				C.2G	4	[C＾2；！c]	generic sp2 carbon
C.2C	5	[$(C＝C)]	alkene carbon
				C.2N	6	[$(C＝[N，P])]	carbon in imine or C＝P group
C.2NP	7	[$(C＝[N+，P+])]	carbon in iminium group
				C.2O	8	[$(C＝[O，S])&！$(C(～[OX1]～[OX1])]	carbon doubly bonded to O or S
C.2U	9	[$(Cl(＝O)O--1)，$(C(＝O)-＝)，$(C(＝O)[F，Cl，Br，I])]	carbon in β-propiolactone，α，β-unsaturated aldehyde or carboxylicacid halide
				C.1G	10	[C＾1]	general sp1 carbon
C.1G2	11	[$(C(＝)＝)]	central carbon in diene
				C.3R	12	[$([#6]1～～[#6]1)，$([#6]1～～*～[#6]1)]	carbon in 3-or 4-membered ring
N.3G	13	[N＾3，$([PX3](-)(-)-*)]	generic sp3 nitrogen
				N.3H	14	[$([N＾3]([#6＾3，#1])([#6＾3，#1])-*)]	amino nitrogen
N.3R	15	[$(N1--1)]	nitrogen in 3-membered ring
				N.2G	16	[N＾2]	generic sp2 nitrogen
N.2C	17	[$(N＝C)]	nitrogen in imine group
				N.1G	18	[N＾1]	general sp1 nitrogen
N.1U	19	[$(N#[*+])，$(N(＝N)＝[C，N])，$(N(＝N)N)]	nitrogen in azide，diazo，and triazenegroup
				N.2OA	20	[$([NX2](a)＝O)]	nitrogen in nitroso group attached toaromatic ring
N.PL	21	[N＾2+0；！$(NC＝O)；！$(N＝*)；！$(N(＝N)N)]	trigonal planar nitrogen
				N.HA	22	[$([NX3]([#6＾3，#1])([#6＾3，#1])a)]	amine attached to aromatic ring
N.AM	23	[$([NX3]([#6＾3，#1])([#6＾3，#1])*＝O)]	nitrogen in amide carbonyl group
				N.PG	24	[N+]	positively charged nitrogen
N.O2	25	[$(N(～O)～O)]	nitrogen in nitro group

N.O2A	26	[$(N(a)(～O)～O)]	nitrogen in aromatic nitro group
				O.3G	27	[O＾3]	generic sp3 oxygen
O.3A	28	[$(Oa)]	oxygen attached to aromatic ring
				O.3H	29	[$([OX2]([#6＾3，#1])[#6＾3，#1])]	ether oxygen
O.3HA	30	[$([OX2]([#6＾3，#1])a)]	oxygen attached to aromatic ring and sp3carbon or hydrogen
				O.3R	31	[$(O1--1)]	oxygen in 3-membered ring
O.2G	32	[O＾2]	generic sp2 oxygen
				O.2C	33	[$(O＝[#6])；！$([OX1]～[SX4](～[OX1]) ～[OX1])；！$([OX1]～[PX4](～[OX1])～[OX1])；！$([OX1]～[CX3]～[OX1])；！$([OX1]～[PX5](～[OX1])(～[OX1])(～[OX1]))]	generic carbonyl carbon
O.2N	34	[$(O＝[#7，P])；！$([OX1]～[SX4](～[OX1])～[OX1])；！$([OX1]～[PX4](～[OX1])～[OX1])；！$([OX1]～[CX3]～[OX1])；！$([OX1]～[PX5](～[OX1])(～[OX1])(～[OX1]))]	nitroso oxygen
				O.2S	35	[$(O＝S)；！$([OX1]～[SX4](～[OX1])～[OX1])；！$([OX1]～[PX4](～[OX1])～[OX1])；！$([OX1]～[CX3]～[OX1])；！$([OX1]～[PX5](～[OX1])(～[OX1])(～[OX1]))]	doubly bonded sulfoxide oxygen

O.O2	36	[$([OX1]～[SX4](～[OX1])～[OX1])，$([OX1]～[PX4](～[OX1])～[OX1])，$([OX1]～[CX3]～[OX1])，$([OX1]～[PX5](～[OX1])(～[OX1])～[OX1])]	oxygen in carboxylate or sulfate group
				O.AG	37	[o]	generic aromatic oxygen
S.3G	38	[$([SX2](-)-)]	generic sp2 sulfur
				S.3H	39	[$([SX2]([#6＾3，#1])[#6＾3，#1])]	sulfur in thiol or thioether
S.2C	40	[$([SX1]＝C)]	sulfur doubly bonded to carbon
				S.2O	41	[$(S(＝O)＝[O，N])，$([SX3](＝[O，N])(-)-)，$([SX3](～[OX1])～[OX1，SX1])，$([SX2](＝O)＝*)]	sulfur in sulfone，sulfonate or sulfategroup
S.AG	42	[s]	generic aromatic sulfur
				P.4O	43	[$([PX4]＝O)，$([PX5]＝O)]	phosphate phosphorous
X	44	[#17，Br，I]	chlorine bromine or iodine

X.A	45	[$([#17，Br，I]a)]	chlorine bromine or iodine attached toaromatic ring
				C.AG	46	[$(c(:a):a)]	generic aromatic carbon
C.AGH	47	[$(c([#6＾3，#1])(:a):a)]	aromatic carbon connected to sp3 carbonor hydrogen
				N.A5	48	[$(n1aaaa1)]	generic aromatic 5-membered ringnitrogen
N.A5L	49	[$([nX3]1aaaa1)]	aromatic 5-membered ring nitrogenwithσlone pair
				N.A6	50	[$(n1aaaaa1)]	generic aromatic 6-membered ringnitrogen
N.A6P	51	[$([n+]1aaaaa1)，$([N+]1＝-＝-＝1)，$([n+]1aaaa1)，$([N+]1＝-＝-*1)]	positive charged nitrogen in aromatic ring
				A.FR	52	[$(a(:a)(:a):a)]	bridge atom in fused cyclic aromatic ringsystem
BAY^*	53		bay-region in polycyclic aromatichydrocarbons
				K^*	54		K-region in polycyclic aromatichydrocarbons
POLY^*	55		polycyclic aromatic system

2)计算分子中每种原子类型的原子描述符：

在MEV中我们共定义三类电子描述符以表征分子亲电性：

a.前线轨道电子密度：原子上的前线轨道电子密度计算可以作为一种有效的电子供体-受体相互作用表征方法，是一种最直接的亲电性指标(Prabhakar，1991；Tuppurainen，et al.，1991)。由于大部分对DNA的进攻是通过亲电反应完成，并且电子转移通常是来自富电子的DNA碱基对，在本方法中我们只计算化学试剂的亲核电子密度(f_N)。

b.电子超离域度：亲电性描述的另一个重要指标是电子超离域度(Fukui，1975；Fukui，et al.，1957)，反映了分子接受或贡献电子密度的能力，可以用作描述占据和非占据轨道的反应性指标。该参数在识别分子软相互作用(Brown and Simas，1982；Kikuchi，1987)和比较不同分子的相应原子的反应性方面(Kikuchi，1987)有重要意义。本方法中我们计算亲核超离域度(S_N)，描述了小分子与另外的反应物(DNA结构)中的亲核中心之间的相互作用(Brown and Simas，1982)。

c.原子净电荷：原子净电荷显然是分子间静电相互作用的驱动力。象电子超离域度这类的轨道描述符代表了动态反应性指标(Franke，1984)，而原子电荷描述了孤立分子基态时的电子性质，可以作为静态反应性指标。当前我们计算三种不同原子净电荷，包括σ-电荷(q_σ)、π-电荷(q_π)和全电荷(q)。这些描述符用于衡量小分子与DNA结构之间的弱分子间相互作用。

我们使用OpenBabel中的OBGastChrg模块根据Gasteiger-Marsili模型对分子分配q_σ电荷(Gasteiger and Marsili，1980)。另外我们设计了一个自定义的C++模块OBHMO根据半经验Hückel方法(Hückel，1931)实现对q_π、fN和SN的计算。对给定的原子r，q_π，r由下式计算得到：

q_{π, r} = k_{r} - Σ_{i = 1}^{m} n_{i} {| C_{r, i} |}^{2}

(式1)

其中m为分子轨道数目；k_r和n_i分别为r原子提供的和处于第i个分子轨道的π-电子的数目；C_r，i是r分子的第i个π分子轨道的原子轨道线性组合(LCAO)的系数。f_N，r由下式计算得到：

(式2)

又被称为最低空分子轨道(LUMO)前沿电子轨道密度；S_N，r由下式给出：

(式3)

其中occ代表已占据分子轨道数目，E_j是第j个非占据分子轨道能级的能量。简单Hückel分子轨道计算只能处理具有π体系的结构而不能考虑分子中的σ键。由于该研究中许多化合物结构中都包括σ键，我们对其进行近似处理当作连接于π体系的甲基杂原子(Streitweiser，1961)。不同原子类型和键类型的库仑积分和共振积分参数值取自文献(Purcell and Singer，1967)。

除了对DNA直接和间接的反应活性，基因毒性也可能由具有平面芳香并环的结构插入到DNA中，引起DNA结构扭曲而造成(Garrett andGrisham，1995)。这种机制引发的基因毒性可用毒效团方法进行准确的预测，Kazius等人使用三个预定义的子结构碎片可以对这种机制下的基因毒性达到93～95％的预测精度(Kazius，et al.，2005)。因此，在MEV表征方法中我们也预先保留三位(如表4所示，即BAY、K和POLY)来保存这些预定义子结构的出现信息。这样，已经验证的毒效团可以便捷的整合到我们的模型中，包含了分子碎片方法易于集成专家知识的优点。

3)根据预定义规则，利用原子描述符和分子碎片模式构建特征向量：

得到对每个原子计算的亲电性描述符和毒效团检索信息之后，下一步骤是构建分子的MEV矢量。具体可以分为以下几步：

a.给定输入分子M，先根据预定规则对其产生一个长度为N×5+3(52×5+3＝263)的浮点数数组V_M，其中N是所有亲电性相关的原子类型数目；

b.将V_M的所有位都初始化为零，并且按每五位合并成为一个子集对应某一特定原子类型；

c.M分子的所有原子根据指定原子类型进行排序，具有同种类型的原子将被映射到相同的数组子集中；

d.对原子类型A所对应的数组子集，每一位分配一种类型的原子描述变量(F_A)，其值为M分子中所有A类型原子的F_A之和。对于预定义毒效团B对应的位，如果分子M中具有B碎片则其值F_B为1，否则为0。

结果得到的数组V_M即为分子M的MEV向量，其中特定位的值及其占据情况是编码了分子亲电性和子结构信息的描述符变量。MEV产生的整个过程可以化合物苯胺为例用(图2)表示。

4)使用基于F-score的递归特征消除(RFE)方法优化MEV向量：

在本方法中我们使用一种基于F-score的递归特征消除(RFE)方法(Guyon，et al.，2002)排列并选择对基因毒性有较高贡献的特征描述变量。对训练集MEV x_i(i＝1，...N)，如果阳性和阴性样本数分别为N₊和N_-，则第j个特征的F-score可由下面公式计算得到：

F (j) = \frac{{({\overset{&OverBar;}{x}}_{j}^{(+)} - {\overset{&OverBar;}{x}}_{j})}^{2} + {({\overset{&OverBar;}{x}}_{j}^{(-)} - {\overset{&OverBar;}{x}}_{j})}^{2}}{\frac{1}{N_{+} - 1} Σ_{i = 1}^{N_{+}} {(x_{i, j}^{(+)} - {\overset{&OverBar;}{x}}_{j}^{(+)})}^{2} + \frac{1}{N_{-} - 1} Σ_{i = 1}^{N_{-}} {(x_{i, j}^{(-)} - {\overset{&OverBar;}{x}}_{j}^{(-)})}^{2}}, - - - (5.5)

这里

和

分别为第j个特征对整体、阳性和阴性样本的平均值。x_i，j ⁽⁺⁾表示第i个阳性样本的第j个特征值。对每个特征变量计算其F-score，值越大表明该特征的判别能力越强。在递归特征消除过程中逐步剔除低F-score打分的特征，同时进行5倍交叉验证跟踪验证正确率的变化，最后选择给出最高的模型预测正确率的最优特征变量子集用于最后的模建。

5)使用支持向量机对训练集分子的基因毒性数据和及其特征MEV进行统计学习：

训练集分子产生的MEV向量作为输入数据进入支持向量机通路以建立并确定最终模型。有关支持向量机的详细理论介绍参见文献(Burges，1998；Vapnik，1995)。本方法中的实施细节描述如下：对给定具有标签值y_i(y_i＝+1 or -1，代表将要被分类的两种类型，在本方法中分别为致突变性和非致突变性化合物)的数据集x_i∈Rⁿ(i＝1，...N)，支持向量机给出下面的决策函数(类型判别器)：

f (x) = sgn (Σ_{i = 1}^{N} γ_{i} α_{i} K (x, x_{i}) + b), - - - (5.1)

这里α_i是待定系数K是内核函数。参数α_i通过最大化下面的拉格朗日表达式训练得到：

\max_{α_{i}} [Σ_{i = 1}^{N} α_{i} - \frac{1}{2} Σ_{i, j = 1}^{N} α_{i} α_{j} γ_{i} γ_{j} K (x_{i}, x_{j})],

其中0≤α_i≤C(i＝1，...N)和

Σ_{i = 1}^{N} α_{i} γ_{i} = 0

(5.2)

我们使用Platt等人的方法得到观测值x_i的类成员f(x_i)的后验概率(Platt，1999)。对所有预测的g(x_i)拟合sigmoid函数以得到具有下面形式的概率表达：

P (y = 1 | g (x_{i})) = p_{i} = \frac{1}{1 + e^{Ag (x_{i}) + B}}, - - - (5.3)

此处A和B由最小化训练集数据的最大似然函数估计得到：

\min_{A, B} [- Σ_{i = 1}^{n} (\frac{y_{i} + 1}{2} \log (p_{i}) + (1 - \frac{y_{i} + 1}{2}) \log (1 - p_{i}))] . - - - (5.4)

标签值和决策值(预测的g(x_i)值)需要彼此独立，所以我们进行5倍交叉验证得到决策值。

本发明中我们使用LibSVM工具包(version 2.81)(Chang and Lin，2001)进行SVM模建，对训练集采用格点搜索策略调整惩罚参数C和径向基函数(RBF)内核参数γ以得到SVM分类器。

6)根据得到的分子基因毒性的概率预测模型，使用C++开发模块化的分子基因毒性预测程序；

即、根据最终的MEV+SVM+RFE模型开发的C++应用程序ZGTOX，可以快速(～200个化合物/秒，P4 2.4G)、准确的对化合物的毒性概率进行预测。

Claims

1.一种基于分子亲电矢量(MEV)和扩展支持向量机(SVM)的基因毒性概率预测方法，适用于根据化合物分子结构信息对该化合物进行虚拟毒性评价和筛选，其特征在于，该方法步骤包括：

a.对数据集中的化合物进行原子分类的步骤；

b.计算分子中每种上述步骤所分类的原子类型的原子描述符；

c.根据预定义规则，利用在上述步骤中获取的原子描述符和分子碎片模式构建特征MEV向量；

d.使用基于F-score的递归特征消除(RFE)方法优化上述的MEV向量；

e.使用支持向量机方法对训练集分子的基因毒性数据和及其特征MEV进行统计学习，得到观测值类成员的后验概率进而得到基因毒性的概率预测模型；

f.根据得到的分子基因毒性的概率预测模型，使用C++开发模块化的分子基因毒性预测程序，实现对化合物的基因毒性的预测。

2.根据权利要求1所述的基因毒性概率预测方法，其特征在于上述步骤a包括：使用基于文本的化学结构提问语言SMARTS描述原子分类；使用SMARTS语法描述的52个原子分类定义规则和3个毒效团子结构定义；使用可编程原子分类PATTY回溯算法完成原子类型指定。

3.根据权利要求1所述的基因毒性概率预测方法，其特征在于上述步骤b中，使用基于根据半经验Hückel方法计算得到原子描述符，用于描述分子的亲电性信息并用于分子基因毒性模建。

4.根据权利要求1或3所述的基因毒性概率预测方法，其特征是上述描述符包括，前线轨道电子密度、电子超离域度、原子π-电荷。

5.根据权利要求1所述的基因毒性概率预测方法，其特征是上述步骤c中、由以下步骤构建特征MEV向量：

(i)给定输入分子M，先根据预定规则对其产生一个长度为N×5+3的浮点数数组V_M，其中N是所有亲电性相关的原子类型数目；

(ii)将V_M的所有位都初始化为零，并且按每五位合并成为一个子集对应某一特定原子类型；

(iii)M分子的所有原子根据指定原子类型进行排序，具有同种类型的原子将被映射到相同的数组子集中；

(iv)对原子类型A所对应的数组子集，每一位分配一种类型的原子描述变量F_A，其值为M分子中所有A类型原子的F_A之和。对于预定义毒效团B对应的位，如果分子M中具有B碎片则其值F_B为1，否则为0。

6.根据权利要求1所述的基因毒性概率预测方法，其特征是上述步骤d中，使用基于F-score的递归特征消除方法优化得到描述分子基因毒性信息的最优描述变量子集。