CN101587510A - 基于复杂抽样和改进决策森林算法的化合物致癌毒性预测方法 - Google Patents

Abstract

本发明涉及一种基于复杂抽样和改进决策森林算法的化合物致癌毒性预测方法，适用于根据有机小分子结构信息对该化合物进行计算致癌毒性评价和虚拟筛选。首先对分子结构分子采用相关力场进行优化并进行电荷计算，对初始训练集中的化合物进行复杂抽样用于训练子集的产生，并根据复杂抽样算法结果固定描述符的组成计算分子中的各种相关描述符；然后使用基于相关矩阵分析和因子分析方法优化描述符池；最后使用改进的决策森林方法对训练集分子的致癌毒性数据和及其相应化学表征进行数据挖掘，得到分类的预测可信区间，致癌毒性预测模型和判断规则。该方法在高通量虚拟筛选和计算致癌毒性评价中具有良好应用前景。

Description

基于复杂抽样和改进决策森林算法的化合物致癌毒性预测方法

技术领域

本发明涉及一种基于复杂抽样和改进决策森林算法的化合物致癌毒性预测计算方法，适用于根据有机化合物分子结构信息对该化合物进行虚拟致癌毒性评价和筛选。

背景技术

毒性问题是后期药物研发失败的一个重要的因素。化合物的致癌毒性指化合物诱发人体内滋生恶性肿瘤或良性肿瘤的一种远期性作用。啮齿类生物测试是目前采用的主要致癌毒性实验方法。然而，该方法存在以下问题：(1)测试费用高(平均测试费用超过两百万美元)；(2)耗时(3至5年)；(3)伦理考虑和公众压力——在研发和测试中减少或不再使用动物。由于计算毒性预测具有低损耗、高通量、虚拟现实等特点，该领域近年已经成为化学，毒理学，生物信息学，化学信息学，药物学，统计学等学科交叉中的一个研究热点(Van deWaterbeemd and Gifford，2003)。与此同时，长期以来生物测试累积的大量实验数据也为小分子致癌毒性预测打下了基础(Yu and Adedoyin，2003)。

致癌化合物根据其分子致癌机制主要分为两大类(Choy，2001)：基因毒性致癌物(Genotoxic carcinogens)和非基因毒性致癌物(Non-genotoxic carcinogens)。前者是通过自身或其活性代谢物与DNA的直接相互作用引起正常细胞转变成癌细胞，经常是跨种属，跨性别，多靶标组织的广泛致癌物。后者则是在癌变初期通过各种不同的机制，如与特定的细胞内受体结合，引起下游细胞分裂异常，经常是在一个种属，一种性别，特定的组织中表现其致癌性。基因致癌物由于其机制相对简单，体外可以采用Salmonella typhimurium细菌实验(即Ames检验)间接帮助致癌性预测，因而对其进行计算毒性预测也相对容易；非基因毒性致癌物则由于牵涉到复杂的生物化学机制网络，如转运，代谢，DNA诱导等等，这是目前计算毒性预测方面的一个重点。

近年来已经有许多不同的计算毒性预测工具问世，依其原理主要可以分为基于统计和基于专家系统的模型(Butina，et al，2002)。前者通过计算各种理化描述符，应用多种统计和模式识别方法，对这些表征小分子化学结构的描述符与毒理属性进行拟合，即在一系列化合物的结构和毒性效应之间推导建立数学意义上的结构与毒性关系(Greene，2002)，代表软件为TOPKAT，MCASE，ToxScope，ToxFilter等等。后者利用根据已知毒性机制建立的规则，把一个分子的某个特殊片段或结构与一个特殊的毒理事件联系在一起，当对新分子进行预测的时候，从存储规则的规则库中调出与该分子相关的规则并依此进行毒理属性的判断，代表软件为DEREK，Ocologic，HazardExpert等等。尽管通过计算进行毒性预测已经成为一个非常流行的方法，当前的商业软件由于其预测准确度尚显不足还不能在无监督情况下使用——近来美国NTP专门对上述的商业程序在致癌毒性的预测能力方面先后进行了两次专项测试(Toivonen，et al，2003)，所选化合物都严格按照NTP啮齿类动物致癌实验协议在被预测后进行实验测定，结果发现，上述软件在致癌毒性预测能力方面均十分有限，以整体准确度作为性能指标衡量时，仅对49～59％的化合物给出正确的预测，说明这些模型在致癌预测性能方面亟待提高。

除商业软件包外，文献中也报道了一些致癌毒性预测的方法(Vendrame，et al，2001)。总体来说，这些模型基本上只采用小样本的同系列化合物，即毒性机制相同或相近的化合物，进行模型的建立，以获得较高的预测正确率，实际上对于药物发现早期对结构多样且机制未知的合成前化合物的致癌性预测，缺乏普遍意义，这表明传统构效关系方法对于当活性和分子描述符之间存在线性相关或化合物为单一作用机制时候时比较有效，但对于模拟多种机制的非线性相关属性效果不明显。目前致癌性预测面临的一个主要挑战就是如何发展更新更好的方法从结构多样的化合物数据集中发现对正确预测有价值的信息(Votano，et al，2004)。

数据挖掘是一个很有潜力的解决上述问题的一个方法。数据挖掘是“从数据中提取出隐含的特别的，过去未知的，潜在有价值的信息”(Dunham，2003)，也是“一门从大量数据或者数据库中提取有用信息的科学”，是一个受到统计学，机器学习，算法语言，数据库系统和信息科学等多个学科影响的交叉领域(Han and Micheline，2000)。其最主要的目的有两个，根据历史预测未来和了解数据中潜在的规律。在药物开发过程中积累的符合标准毒理学实验测试结果的化合物数据库为采用数据挖掘方法进行致癌性预测提供了良好的数据挖掘基础。

除了性能方面的局限，当前的计算毒性预测方法经常会碰到随机预测的问题，即改变哪怕很小的训练数据集的组成都可能作出截然相反的预测结果。因此，一个能够产生具有统计意义的预测结果和提炼出稳定判断规则的计算方法对于指导研究者作出可靠决策具有重要意义。

发明内容

本发明的目的是提出一种能够避免随机预测又同时提供稳定致癌性判断规则的基于决策森林的致癌毒性评价方法。

该目的主要是采用分子描述，复杂抽样和改进的决策森林方法，通过对现有分子结构和相关致癌毒性数据进行数据挖掘达到：

1)结构与毒性关系模型的建立要求选择与属性最相关的结构片断或化学特征。对于致癌毒性，描述符的选择需要充分考虑其分子致癌机制。如前所述，小分子致癌毒性的产生机制可以概括为如下两个方面：a)自身或其代谢物与DNA的直接相互作用引起正常细胞转变成癌细胞；b)通过复杂的生物化学机制网络，如转运，代谢，DNA诱导等，小分子与特定的细胞内受体结合，引起下游细胞分裂异常。这些效应影响正常的细胞分裂周期，最终导致致癌毒性。当前的研究表明，致癌毒性可能与分子的大小、形状、电子信息、极化度等多种化合物内在属性相关。Votano等人的QSAR研究(Votano，et al，2004)显示，化合物中杂原子的电负性和共价半径对致癌毒性的产生也有重要作用。基于上述原因，我们从五类约49个符合上述特征的描述符中分别经过相关矩阵分析，因子分析(Factor Analysis，FA)等方法的挑选获得18个描述符用于建模和规则推导。

2)决策森林是一系列标准决策单树的集合，由于每一株树都能对任意一个观测值给出预测值，大多数单树的预测结果将是稳定的预测结果，而最终预测值是所有决策单树即决策森林预测结果的平均值。这样的策略能有效降低决策单树出现随机预测的缺点。理论上，构成决策森林的每株决策单树可以通过a)改变算法建立；b)固定训练集，改变变量建立；c)固定变量改变训练集的组成建立(Breiman，2001)。本研究采用第三种方式建立决策森林。这是因为决策单树依赖于训练集的组成，改变训练集的组成(对于决策森林而言就是训练子集)意味着改变树形，而采用合理方式生成的训练子集有助于决策单树之间彼此合理交叉，使训练子集的化学空间能够有效重叠在一起，在充分代表原有总训练集化学空间的同时，能够对其中组成重复进行预测，降低决策单树中出现随机预测的不利影响，并依据每一株决策单树中的描述符测试值提取出稳定的判断依据(规则)。基于这个原因，我们设计编写了一个采用复合抽样算法(Chaudhuri and Stenger，1992)的程序用于优化训练子集的产生。相对决策单树而言，决策森林在药物发现中的应用还不多(Zhang，et al，2003)，尤其是用于计算毒性预测，尚无报道。

3)数据集由预测毒理学官方网站获得(http://www.predictive-toxicology.org/ptc/，2000-2001年预测毒理学竞赛数据集)(Helma，2001)，分成训练集(下称内部训练集)和测试集(下称外部训练集)两部分。内部数据集的活性数据来自美国国家毒理计划(US National Toxicology Progarm，US NTP)测定结果，包括502个由“简化分子输入行录入系统”(simplified molecular input lineentry system，SMILES)格式表示的经四种啮齿类动物(雄性大鼠，雌性大鼠，雄性小鼠，雌性小鼠)致癌实验测定相应致癌性的化合物。外部训练集包括282个来自美国食品和药品管理局(US Food and DrugAdministration，FDA)的药物的经同样四种啮齿类动物致癌性实验结果。根据分子特征挖掘软件(Molecular Feature Miner，MFM)的计算结果(Toivonen，et al，2003)，NTP数据集中包含50.60％的为MFM完全已知的结构片段，在FDA数据集中该数字是8.07％，因此表明两个数据集合之间有足够的结构多样和差异性。由于雄性大鼠实验中有更均衡的致癌物与非致癌物，可以避免训练时发生计算偏斜，因此我们选择该组数据建立模型，其中，带有EE(Equivocal Evidence，证据不足)，IS(Inadequated Study，研究不足)和E(Equivocal，不确定)标记的化合物从数据集中剔除，带有CE(Clear Evidence，明确证据)，SE(SomeEvidence，部分证据)和P(positive，阳性)标记的化合物赋值为+1，统一标记为致癌物，带有NE(No Evidence，无证据)和N(Negative，阴性)标记的化合物赋值为-1，统一标记为非致癌物。在进一步除去不正确的SMILES字串，最后内部验证集中包括342个分子(192个非致癌物，150个致癌物)，外部验证集中包含217个分子(159个非致癌物，58个致癌物)。

具体地，本发明提供一种基于复杂抽样和改进决策森林算法的化合物致癌毒性预测方法，适用于根据有机化合物结构信息对该分子进行致癌毒性计算评价和筛选，所述方法包括如下步骤：

a.使用复杂抽样算法完成训练子集产生，使用如图5所示的采用三层循环基于系统时间的随机数发生器进行化合物的随机挑选程序流程，并使用C程序完成训练子集抽样；

b.根据复杂抽样算法结果固定描述符的组成计算分子中的各种相关描述符，用于分子的表征并用于分子致癌毒性建模，所述描述符包括电子描述符、空间描述符、结构描述符、热力学描述符和拓扑描述符；

c.使用基于相关矩阵分析和因子分析方法优化描述符池；

d.使用改进的决策森林方法对训练集分子的致癌毒性数据及其相应化学表征进行数据挖掘，得到分类的预测可信区间和致癌毒性预测模型；

e.根据改进的决策森林方法使用C语言开发模块化的有机化合物致癌毒性预测程序；

f.根据得到的致癌毒性预测模型，统计所有决策单树中描述符的出现频率，发现与致癌性最相关的描述符与相关规则：Jurs-DPSA-2点电荷比表面指数(Jurs-DPSA-2)，最高占用分子轨道(HOMO)，油水分配系数(AlogP)，原子极化度和(Apol)和Balaban-JX分子连接指数(Balaban-JX)与致癌物最为相关；Jurs-DPSA-2点电荷比表面指数(Jurs-DPSA-2)，最高占用分子轨道(HOMO)和分子折射率(Molref)与非致癌物最为相关。

有益效果

本发明的有益效果体现在，一方面，我们整合了统计知识-相关矩阵分析和因子分析，并考虑了与致癌毒性引发机制高度相关的分子属性进行描述符的挑选；另一方面，我们改进现有的决策森林方法，固定描述符的组成并通过复杂抽样算法进行训练子集的产生，将该方法用于致癌毒性数据挖掘，展现了决策森林方法的三个新颖的优点：即其独特的树型结构产生均衡的预测能力，化合物平均得分值的高低可以作为结果可信程度的判断依据，以及可以回溯判断依据，形成预测规则模型。计算结果表明我们的方法对化合物致癌毒性预测明显优于现有方法，测试集中化合物的改变不会明显波动我们建立的模型的预测性能，说明这种改进的决策森林方法对药物分子致癌毒性预测建模是稳定有效的。根据模型开发的C应用程序DDDC-CTOX可快速的对化合物致癌毒性进行预测，在高通量虚拟筛选和计算致癌毒性评价中具有良好应用前景。

附图说明

图1：DF和SVM采用FA挑选的描述符对来自外部验证集的所有子集进行预测的TP(a)和TN(b)对子集中化合物数的散点图；

图2：DF采用FA挑选的描述符对内部验证集中342个化合物的正确预测分布，其中正确预测率是平均预测分类值的函数，化合物平均得分小于0被认为是非致癌物，大于0被认为是致癌物；

图3：决策森林中所有描述的出现频率，其中(a)为致癌物预测描述符；(b)为非致癌物预测描述符；

图4：SSDT采用6个由DF选出的描述符对内部验证集进行预测建立的树型结构；

图5：训练子集产生流程图，其中三层循环用于控制最后训练子集的产生：第一层循环，数据类型time_t确保抽样的真正随机性；第二层循环，变量i记录子集的数量；第三层循环，变量j实时记录所有子集中每个元素的出现次数。

具体实施方式

下文将对本发明进行详细描述。

其中，基于复杂抽样和改进决策森林算法的化合物致癌毒性预测模型的建立主要涉及五个步骤：

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

由于致癌毒性可能与分子的大小、形状、电子信息、极化度、化合物中杂原子的电负性和共价半径对等多种化合物内在属性相关，我们首先用如下五类约49个符合上述特征的描述符进行化合物属性表征：a)电子描述符，所述电子描述符为原子极化度和(Apol)，偶极矩(Dipole)，最高占用分子轨道(HOMO)，最低未占用分子轨道(LUMO)，超离域能(Sr)；b)空间描述符，所述空间描述符为转动半径(RadOfGyration)，分子表面积(Area)，分子体积(Vm)，分子密度(Density)，惯性矩(PMI)，Jurs点电荷比表面指数(Jurs indices)，分子形状投影(Shadow indices)；c)结构描述符，所述结构描述符为旋转键(Rotlbonds)，氢键受体(Hbondacceptor)，氢键供体(Hbond donor)，分子量(MW)，手性中心(chiralcenters)；d)热力学描述符，所述热力学描述符为油水分配系数(Alogp)，分子折射率(MolRef)，水中去溶剂化自由能(Fh2o)，油中去溶剂化自由能(Foct)，生成热(Hf)；e)拓扑描述符，所述拓扑描述符为Balaban分子连接指数(Balaban JX和JY)，Wiener化学键指数(Wiener，CHI分子连接指数(CHI indices)，Kappa形状指数(Kappa indices)，Hosoya分子图指数(log Z)，Zagreb分子图指数(Zagreb)。所有分子采用UNIVERSERAL1.02力场进行优化(UNIVERSERAL ForceField，UFF)：分子的电荷采用电荷平衡公式计算，最大迭代次数和收敛标准分别设为500和5.0E-04；最大优化循环次数设为50,000。分子描述符采用Cerius 2(http://www.accelrys.com)软件计算。

2)使用相关矩阵分析和因子分析优化描述符集：

计算整个描述符集的相关矩阵，在此过程中，彼此之间相关系数大于0.85的描述符被舍弃以保持描述符的独立性，产生的25个独立描述符被进一步用因子分析进行数据降维。因子分析属于无监督的数据降维方法(Malinowski and Howery，1984)，即无须目标属性为连续性数值，非常适合我们研究中致癌属性为+1/-1时的情况，用于挑选那些起主要贡献的描述符。下面以对这25个独立描述符进行数据降维为例，说明因子分析具体的步骤和计算原理：

a)先计算内部验证集中所有342个化合物由25个描述符代表的协方差矩阵。该矩阵由式(1.1)计算：

$A_{\mathrm{ij}}=\frac{\underset{k=i}{\overset{N}{\mathrm{\Σ}}}{X}_{\mathrm{ki}}{X}_{\mathrm{kj}}-\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{X}_{\mathrm{ki}}\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{X}_{\mathrm{kj}}/N}{N-1}---\left(1.1\right)$

i，j＝1，2，...，P

其中，A表示一个包括N行化合物和P列描述符的矩阵，i，j是描述符下标，k是化合物下标。

b)用雅可比方法(Jacobi method)(Causey，1958)计算A的实对称矩阵的特征值λ与相应特征向量U，从而得综合向量F：

F＝U′X＝U′[X₁，X₂，...，X_p]    (1.2)

c)根据特征值λ计算所有主因子方差的累计方差百分比：

$(\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_i/\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_i)\&times;100\%\&GreaterEqual;M(m<p)---\left(1.3\right)$

M是第m个主因子的贡献率，该值决定主因子的排序。

d)由最大正交旋转法(Varimax algorithm)(Harman，1976)计算每个描述符对主因子的贡献，根据每个描述符的因子载荷(＞0.7)挑选描述符：

$\max \left(V_j\right)=\underset{i=1,m}{\mathrm{\&Sigma;}}[\frac{1}{p}\underset{j-1,p}{\mathrm{\&Sigma;}}{\left(b_{\mathrm{ij}}^2\right)}^2-\frac{1}{p^2}{\left(\underset{j=1,p}{\mathrm{\&Sigma;}}{b}_{\mathrm{ij}}^2\right)}^2]---\left(2.4\right)$

其中，b_ij是载荷矩阵的因子载荷，max(V_j)是相应因子载荷平方的最大方差，m，p为载荷矩阵的行与列，分别对应综合向量数与描述符数。经过因子分析挑选描述符，最终确定18个描述符用于后续的决策森林计算。

3)训练子集产生：

如何合理生成训练子集对于各决策单树之间产生有效重叠的化学空间至关重要。我们编写了一个采用复合抽样算法(Chaudhuri andStenger，1992)的程序用于优化训练子集的产生。图5是该程序的流程图。程序分为三个模块分别实现三个功能，(1)在训练子集中重现原训练集的致癌物与非致癌物的比例。程序用一个二维数组代表所有342个化合物以及相应致癌属性(数组的行元素为1-342之间的自然数，列元素为+1/-1，分别代表致癌物和非致癌物)，统计数组的列元素以确定致癌物与非致癌物的比例，在此限制下，采用基于系统时间的随机数发生器进行自然数的随机挑选；(2)确保同一个训练子集中不出现重复的化合物。为实现这一个目的，每个训练子集按第一元素升序排序，如果任意相临元素相等，则舍弃整个子集；(3)保证每个化合物的出现次数有统计意义。在符合前两个条件的训练子集产生后，对所有子集中的自然数进行加和，当所有自然数的出现次数大于等于50的时候，则认为产生有效的训练子集用于决策森林的建立。在本次研究中，致癌物与非致癌物的比例为0.78∶1(150/192)，每个子集中的化合物数人为设定为280(约占总训练集中342个化合物的80％)。总共生成80个子集即对应80株决策单树。

4)使用改进的决策森林算法对训练集分子的致癌毒性数据和及其表征描述符进行数据挖掘：

我们采用循环化的决策树C4.5递归区分(Recursive Partitioning，RP)算法处理训练子集及相应固定的描述符数组。RP是决策树的核心，是一个根据划分类别的信息熵进行二叉区分的有监督算法，从样本中抽提出同质元素，并建立基于树型的分类规则(Breiman，et al，1984)。建树是一个递归过程，可用下面的语句表示。其中，在第五步根据最佳区分限对节点进行分割。

1)  BuildingTree(Training Set T)

2)      Partition(T)

1)  Partition(Data S)

2)     If(all points in S are in the same class)then return；

3)     Else

4)        Evaluate splits for each attribute A；

5)        Use best split found to partition S into S1 and S2；

6)        Splitting(S1)；

7)   Splitting(S2)；

设有包括N个对象的训练集为S，所有对象都隶属于K类，用C1，C2，...Ck表示；每个对象均由带有V个不同值的属性A{a1，a2，...av}描述；P(Ci，S)为属于Ci类的对象数占总对象数的比例，那么

某一类对象的信息量(Information)和系统总熵(Entropy)如下：

Information＝-log₂(P(C_i，S))      (2.5)

$\mathrm{Entropy}=-\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}P(C_i,S){\log }_2\left(P(C_i,S)\right)---\left(2.6\right)$

属性A可以把S分成V个子集{S1，S2，...，Sv}，其中Sj＝{x|x∈S&A＝aj}。这些子集表示从集合S出发的所有树枝。令Sij为在Sj中类为Ci的记录个数，则按A的每个属性值进行分割后的信息量，也就是系统区分熵Entropy(T)为：

$\mathrm{Entropy}\left(T\right)=\underset{j=1}{\overset{v}{\mathrm{\&Sigma;}}}((S_{1j}+S_{2j}+...+S_{\mathrm{mj}})/S)\&times;I(S_{1j}+S_{2j}+...+S_{\mathrm{mj}})---\left(2.7\right)$

其中，((S_1j+S_2j+...+S_mj)/S)表示第j个子集的权重，s＝|S|；子集的总熵I为：

$I(S_{1j}+S_{2j}+...+S_{\mathrm{mj}})=-\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{P}_{\mathrm{ij}}{\log }_2\left(P_{\mathrm{ij}}\right)---\left(2.8\right)$

又因为熵越大，表示系统越混乱。而分类的目的是提取系统信息，使系统向更加有序、有规则组织的方向发展，所以，当满足下式时，达到最佳的分割方案，递归区分停止，实现熵减少量Gains(A)最大：

$\mathrm{Gains}\left(A\right)=\frac{{\log }_2(N-1)}{N}+\frac{\mathrm{\&Delta;}(A,T;S)}{N}---\left(2.9\right)$

这里，Gains(A)＝Entropy-Entropy(T)，

Δ(A，T；S)＝log₂(3^K-2)-[K·Entropy(S)-K₁·Entropy(S₁)-K₂·Entropy(S₂)-..-K_m·Entropy(S_m)]

N为对象数，Ki为相应子集的对象数。

5)致癌规则的提出：

由于决策树能够回溯判断依据，即显式提出进行相应预测的依据形成规则。我们统计了所有80株决策单树中描述符的出现频率，提出下述6个描述符是本方法中与致癌毒性预测最相关的表征：Jurs-DPSA-2，HOMO，AlogP，Apol和Balaban-JX是预测致癌物的关键描述符，Jurs-DPSA-2，HOMO和Molref是预测非致癌物的关键描述符。这6个描述符的内在属性表述如下：a)Jurs-DPSA-2是空间描述符，它综合分子的形状和电子信息表示其特征，数值上，它等于总电荷除以正电荷的溶剂可及表面减去总电荷除以负电荷的溶剂可及表面。MolRef是热力学描述符，它反映的是分子大小和极化度。也就是说，此两个描述符体现的是分子结构特征，从毒效动力学(Toxicodynamics)的角度，能够用来定义可能的分子间的相互作用，如毒物与靶标分子如DNA的相互作用。b)HOMO属于电子描述符，表示最高占用分子轨道能。该描述符是化合物电负性的量度。在基因毒性致癌机制上，化合物要与亲核性的DNA以共价键形式形成加合物需要有一个亲电中心，而对于具有越负的HOMO能量的化合物，其亲核性很强，因此，这样的化合物有较大的可能性是非致癌物。Apol是另一种电子描述符，是衡量原子极化度的指标，其数值越小，表示存在杂原子的程度越高，因此存在具有致癌性的极化官能团的可能性越高(Kier.and Hall，1999)。C)Balaban-Jx是一个拓扑描述符，着重反映杂原子的电负性和共价半径，这对应于很多含杂原子官能团与致癌的化学结构相关(Votano，et al，2004)，如亚硝胺，三卤甲烷，某些芳香胺等等。AlogP属于热力学描述符，即油水分配系数，与分子的疏水性直接相关，而分子的疏水性是影响分子顺利运输到作用靶点的重要因素，所以，从毒物代谢动力学(Toxicokinetics)的角度，在毒性预测时候考虑这一描述符很有必要。实际上，LogP已经在很多文献中用于构毒关系研究(Lipnick，1995)。

本发明中我们基于复合抽样和改进决策森林算法进行建模，对训练集采用固定描述符池的，设定裁剪因子为中度裁剪，最大树深为5，其余参数为可调以优化决策单树的方式以得到致癌毒性预测分类模型。单株决策单树的挑选依据为：采用10倍交叉验证时必须对训练子集有70％以上的整体预测正确率，内部验证集的真阳性率大于0.5，假阳性率低于0.5。

本发明中我们首先对分子结构分子采用相关力场进行优化并进行电荷计算，对初始训练集中的化合物进行复杂抽样用于训练子集的产生，并根据复杂抽样算法结果固定描述符的组成计算分子中的各种相关描述符；然后使用基于相关矩阵分析和因子分析方法优化描述符池；最后使用改进的决策森林方法对训练集分子的致癌毒性数据和及其相应化学表征进行数据挖掘，得到分类的预测可信区间，致癌毒性预测模型和判断规则。

为进一步客观评价模型的预测能力，我们将建立的决策森林(Decision Forest，DF)，支持向量机(Support Vector Machine，SVM)和TOPKAT模型对外部验证集(如前述，由来自美国食品和药品管理局的282个药物组成)作为实例进行预测，所获得的结果与来自预测毒理学竞赛中最好的Gonzales模型的结果(Toivonen，et al，2003)放在一起进行比较，见表1。

表1我们的决策森林模型与Gonzales模型和商业程序TOPKAT对外部验证集的预测结果比较

^aGonzales模型是预测毒理学竞赛中表现最好的预测模型，该模型采用的外部验证集包括133个非致癌物和52个致癌物31。

通过表1的结果表明，DF是其中最好的预测模型，其致癌物的真阳性预测正确率超过62％，优于Gonzales模型30％以上，SVM居次，其致癌物的真阳性预测正确率比随机预测略高，为51.7％。我们还使用TOPKAT对外部验证集进行评价，如表1所示的结果，该程序对致癌化合物与Gonzales模型一样并没有有效的致癌物预测能力，真阳性率仅为32.0％，显著的低于我们的决策森林模型。从表1还可以发现，TOPKAT和Gonzales模型在识别非致癌物方面有明显的倾向性。由于假阴性结果有可能接受进一步的实验研究，这样实际上等于并没有经过有效评价，因此在评价一个模型预测能力的时候，更应关注其真阳性率的高低。我们的决策森林模型满足这一要求，同时显示了该模型在识别致癌物和非致癌物方面具有同样均衡的预测能力。

由于训练集与测试集之间的结构差异对模型预测表现会有影响，但这种影响总是希望越小越好。为进一步的分析我们模型的预测稳定性，我们用前面介绍过的复杂抽样程序对外部测试集进行抽样产生了5个测试子集，然后提交DF和SVM模型进行预测。这五个测试子集中的化合物分别是40，80，120，160和200，其中的致癌物与非致癌物的比例都固定为原外部验证中这两类化合物的比例0.36∶1(58/159)。DF和SVM对这五个测试子集和原外部验证集的预测结果列于表2，TP/TN对化合物数的散点图见图1。

表2DF和SVM对来自外部验证集的所有子集中的致癌物与非致癌物的预测表现a

a所有模型的预测表现都用真阳性率(True Positive，TP)，真阴性率(TrueNegative，TN)，假阳性率(False Positive，FP)，假阴性率(False Negative，FN)四个参数表示。“+FA”代表方法使用的是FA挑选出的18个描述符。

通过表2和图1中的结果可以发现，对任一测试子集，这DF和SVM的四个衡量参数都变化不大，DF在致癌物预测上稳定的优于SVM，而SVM比DF有稍好一些的非致癌物预测能力。因此，我们认为这两个数据挖掘方法都能够有效进行化合物致癌毒性预测，而DF以其更均衡的表现更加值得注意。此外，因为测试集中化合物的改变并没有明显引起模型预测表现的波动，这也说明训练集与测试集的结构差异不是影响模型表现的主要因素(Toivonen，et al，2003)。

除了分类预测外，人们在实际应用时通常对分类的预测可信度和判断规则更感兴趣。对于预测的可信度而言，由于DF最终预测值是所有决策单树预测结果的平均值，故平均得分值的高低可以作为结果可信程度的判断依据。图2提供了DF对内部验证集中342个化合物正确预测的分布。对非致癌物预测，[-0.4，-1.0]区域有约为88.9％(120/135)的正确预测率，对致癌物预测，[0.4，1.0]区域有约为77.1％(121/157)的正确预测率。主要的错误或者模糊预测区域在[-0.4，0.4]之间，总共只有～46％(23/50)的预测正确率。所以，上述预测区域可以作为预测可信度高低的判据：在横坐标两端可信度最高，中间最低，整体表现为一个近似U形的曲线。

对判断规则，相对于其他算法，DF能够回溯判断依据，即能够显式提出进行相应预测的依据，形成规则。我们统计了所有决策单树中描述符的出现频率，以发现与致癌性最相关的描述符。图3给出了相应统计结果。

从图3a可以发现，Jurs-DPSA-2，HOMO，AlogP，Apol和Balaban-JX是预测致癌物的关键描述符，图3b则表明Jurs-DPSA-2，HOMO和Molref是预测非致癌物的关键描述符。因此，这6个描述符是本研究中与致癌毒性预测最相关的描述符。

另外为测试其表现，我们用基于这6个描述符的SSDT和SVM模型对内部和外部验证集进行预测，发现它们同样具有较好的预测能力，图4提供了相应建立的决策单树。

综上所述，本发明中我们基于电子，空间，结构，热力学和拓扑等五类描述符，设计编写了一个采用复合抽样算法的决策森林程序，发展了一套用于药物分子致癌毒性预测模型。一方面，我们整合了统计知识-相关矩阵分析和因子分析，并考虑了与致癌毒性引发机制高度相关的分子属性进行描述符的挑选；另一方面，我们改进现有的决策森林方法，固定描述符的组成并通过复杂抽样算法进行训练子集的产生，将该方法用于致癌毒性数据挖掘，展现了决策森林方法的三个新颖的优点：即其独特的树型结构产生均衡的预测能力，化合物平均得分值的高低可以作为结果可信程度的判断依据，以及可以回溯判断依据，形成预测规则模型。计算结果表明我们的方法对化合物致癌毒性预测明显优于现有方法，测试集中化合物的改变不会明显波动我们建立的模型的预测性能，说明这种改进的决策森林方法对药物分子致癌毒性预测建模是稳定有效的。根据模型开发的C应用程序DDDC-CTOX可快速的对化合物致癌毒性进行预测，在高通量虚拟筛选和计算致癌毒性评价中具有良好应用前景。

参考文献

Breiman，L.Random forests.Mach.Learn.，2001a，45：5-32.

Breiman，L.，Friedman，J.，Stone，C.，and Olshen，R.Classification andRegression Trees，1984，Wadsworth International Group，Belmont.Van deWaterbeemd，H.and Gifford，E.ADMET in silico modeling：towardsprediction paradise？Nature Rev.Drug Discov.，2003，2：192-204.

Butina，D.，Segall，M.D.and Frankcombe，K.Predicting ADMEproperties in silico：methods and models.Drug Discov.Today.，2002，7(11Suppl)：S83-88.

Causey，R.L.Computing eigenvalues of non-Hermitian matrices bymethods of Jacobi type.J.Soc.Indust.Appl.Math.，1958，6：172-181.

Chaudhuri，A.and Stenger，H.，Survey sampling，1992，Dekker，N.y.

Choy，W.N.Genotoxic and non-genotoxic mechanisms ofcarcinogenesis.In：Genetic toxicology and cancer risk assessment(ed.Choy，W.N.)，2001，PP.47-72.New York，NY，USA：Marcel Dekker.

Dunham，M.H.Data mining.Introductory and advanced topics.2003，Prentice Hall.

Greene，N.Computer systems for the prediction of toxicity：an update.Adv.Drug Deliver Rev.，2002，54：417-431.

Han，J.and Micheline，K.Data mining：concepts and techniques.2000，Morgan Kaufmann.

Harman，H.H.Modern Factor Analysis(3rd Ed.).1976，University ofChicago Press.

Helma.，C.，King，R.D.，Kramer，S.，and Srinivasan，A.The PredictiveToxicology Challenge 2000-2001.Bioinformatics，2001，17：107-108.

Kier，L.B.and Hall，L.H.Molecular Structure Description；TheElectrotopological State，1999，Academic Press，SanDiego，CA.

Lipnick，R.L.Structure-activity relationships.In：Fundamentals ofAquatic Toxicology(Rand，G.M.，second Ed.)pp 609-655，1995，Taylor&Francis，London.

Malinowski，E.R.，and Howery，D.G.Factor Analysis in Chemistry，1980，Wiley-Interscience，New York.

Toivonen，H.，Srinivasan，A.，King，R.D.，Kramer，S.，and Helma，C.Statistical evaluation of the Predictive Toxicology Challenge 2000-2001.Bioinformatics，2003，19：1183-1193.

Vendrame，R.Braga，R.S.，Takahata，Y.and

，D.S.Structure-careinogenic activity relationship studies of polycyclic aromatichydrocarbons(PAHs)with pattern-recognition methods.J.Mol.Struct.(Theochem)，2001，539：253-265.

Votano，J.R.，Parham，M.，Hall，L.H.，Kier，L.B.，Oloff，S.，Tropsha，A.，Xie，Q.and Tong，W.Three new consensus QSAR models for theprediction of Ames genotoxicity.Mutagenesis，2004，19：365-377.

Zhang，H.，Yu，C.and Singer，B.Cell and tumor classification usinggene expression data：construction of forests.Proc.Natl.Acad.Sci.USA.，2003，100：168-4172.

Yu，H.and Adedoyin，A.ADME-Tox in drug discovery：integration ofexperimental and computational technologies.Drug Discov.Today，2003，8：852-861.

Claims (4)

1、一种基于复杂抽样和改进决策森林算法的化合物致癌毒性预测方法，其特征在于，所述方法包括如下步骤：

a.使用复杂抽样算法完成训练子集产生，采用三层循环基于系统时间的随机数发生器进行化合物的随机挑选，并使用C程序完成训练子集抽样；

b.根据复杂抽样算法结果固定描述符的组成计算分子中的各种相关描述符，用于分子的表征并用于分子致癌毒性建模，所述描述符包括电子描述符、空间描述符、结构描述符、热力学描述符和拓扑描述符；

c.使用基于相关矩阵分析和因子分析方法优化描述符池；

d.使用改进的决策森林方法对训练集分子的致癌毒性数据及其相应化学表征进行数据挖掘，得到分类的预测可信区间和致癌毒性预测模型；

e.根据改进的决策森林方法使用C语言开发模块化的有机化合物致癌毒性预测程序；

f.根据得到的致癌毒性预测模型，统计所有决策单树中描述符的出现频率，发现与致癌性最相关的描述符与相关规则。

2、根据权利要求1所述的基于复杂抽样和改进决策森林算法的化合物致癌毒性预测方法，其中，所述步骤a中的三层循环包括如下三个模块：(1)在训练子集中重现原训练集的致癌物与非致癌物的比例的模块；(2)确保同一个训练子集中不出现重复的化合物的模块；(3)保证每个化合物的出现次数有统计意义的模块。

3、根据权利要求1所述的基于复杂抽样和改进决策森林算法的化合物致癌毒性预测方法，其中，在所述步骤b中，所述电子描述符为原子极化度和、偶极矩、最高占用分子轨道、最低未占用分子轨道、超离域能；所述空间描述符为转动半径、分子表面积、分子体积、分子密度、惯性矩、Jurs点电荷比表面指数、分子形状投影；所述结构描述符为旋转键、氢键受体、氢键供体、分子量、手性中心；所述热力学描述符为油水分配系数、分子折射率、水中去溶剂化自由能、油中去溶剂化自由能、生成热；所述拓扑描述符为Balaban分子连接指数、Wiener化学键指数、CHI分子连接指数、Kappa形状指数、Hosoya分子图指数、Zagreb分子图指数。

4、根据权利要求1所述的基于复杂抽样和改进决策森林算法的化合物致癌毒性预测方法，其中，在所述步骤f中，与致癌物最为相关的描述符为Jurs-DPSA-2点电荷比表面指数、最高占用分子轨道、油水分配系数、原子极化度和和Balaban-JX分子连接指数，与非致癌物最为相关的描述符为Jurs-DPSA-2点电荷比表面指数、最高占用分子轨道和分子折射率，Jurs-DPSA-2点电荷比表面指数综合分子的形状和电子信息表示其特征，数值上，它等于总电荷除以正电荷的溶剂可及表面减去总电荷除以负电荷的溶剂可及表面；最高占用分子轨道能是化合物电负性的量度；油水分配系数与分子的疏水性直接相关；原子极化度和是衡量原子极化度的指标，其数值越小，表示存在杂原子的程度越高，因此存在具有致癌性的极化官能团的可能性越高；Balaban-Jx分子连接指数反映杂原子的电负性和共价半径，这对应于很多含杂原子官能团与致癌的化学结构相关；分子折射率反映分子大小和极化度。

Images