基于支持向量机方法的HERG钾离子通道抑制剂分类模型

基于支持向量机方法的HERG 钾离子通道抑制剂分类模型李平2谈宁馨2饶含兵1李泽荣1,*陈宇综3(1四川大学化学学院,成都610065;2四川大学化学工程学院,成都610065;3Department of Pharmacy,National University of Singapore,Singapore 117543)摘要：对human ether -a -g ò-g òrelated genes(HERG)钾离子通道(钾通道)抑制剂,计算了表征分子组成、电荷分布、拓扑、几何结构及物理化学性质等特征的1559个分子描述符,采用Fischer Score(F -Score)排序过滤和Monte Carlo 模拟退火法相结合从中筛选与HERG 钾通道抑制剂分类相关的分子描述符.采用支持向量机(SVM)方法,分别以IC 50=1.0、10.0μmol ·L -1为分类标准,建立了三个分类预测模型.对367个训练集分子,用五重交叉验证,得到正、负样本的平均预测精度分别为84.8%-96.6%、80.7%-97.7%,其总的平均预测精度为87.1%-97.2%,优于其它文献报道结果.对97个外部测试集分子,所建三个模型的总样本预测精度在67.0%-90.1%之间,接近或优于其它文献报道结果.关键词：支持向量机;HERG 钾通道抑制剂;Monte Carlo 模拟退火法中图分类号：O641Classification Models for HERG Potassium Channel InhibitorsBased on the Support Vector Machine ApproachLI Ping 2TAN Ning -Xin 2RAO Han -Bing 1LI Ze -Rong 1,*CHEN Yu -Zong 3(1College of Chemistry,Sichuan University,Chengdu 610065,P.R.China ;2College of Chemical Engineering,Sichuan University,Chengdu 610065,P.R.China ;3Department of Pharmacy,National University of Singapore,Singapore 117543)Abstract ：We calculated 1559molecular descriptors including constitutional,charge distribution,topological,geometrical,and physicochemical descriptors to characterize the molecular structure of human ether -a -g ò-g òrelated genes (HERG)potassium channel inhibitors.A hybrid filter/wrapper approach combing the Fischer Score (F -Score)and Monte Carlo simulated annealing was used to select molecular descriptors relevant to the discrimination of HERG potassium channel inhibitors.Three classification models with threshold values of IC 50=1.0,10.0μmol ·L -1,respectively,were built using the support vector machine (SVM)approach.Models developed from 367training set molecules were validated through 5-fold cross -validation (CV)and the average prediction accuracies were 84.8%-96.6%,80.7%-97.7%,and 87.1%-97.2%for the positive,negative,and overall samples,respectively,which showed better performance than models previously reported in literature.Overall prediction accuracies for the three models using an external test set of 97molecules were between 67.0%and 90.1%,which were close to or better than the results reported in literature.Key Words ：Support vector machine;HEGR potassium channel inhibitor;Monte Carlo simulated annealing[Article]物理化学学报(Wuli Huaxue Xuebao )Acta Phys.-Chim.Sin .,2009,25(8)：1581-1586Received:February 20,2009;Revised:April 20,2009;Published on Web:June 1,2009.*Corresponding author.Email:lizrscu@;Tel:+86-28-85406139鬁Editorial office of Acta Physico -Chimica SinicaHERG 基因编码的钾离子通道存在于人类心室和心房肌细胞中,其HERG 通道编码心脏快速延迟整流钾电流I Kr ,为心肌细胞动作电位三相复极期的主要外向钾电流,在心脏动作电位复极化过程中发挥着重要作用[1].此外,快速延迟整流钾通道(I Kr /HERG 通道)是III 类抗心律失常药物的作用靶点,被某些药物作用时,会引起长QT 间期综合征(LQTS),进而导致尖端扭转性室性心动加速(TdP,Torsade deAugust 1581Acta Phys.-Chim.Sin.,2009Vol.25Pointes),甚至猝死[2].因此,药物通过阻滞HERG通道产生心律失常的不良反应已成为药物开发中一个必须考虑的问题.药物从研发到上市的过程非常漫长,其中很多化合物因为药物吸收、分布、代谢、排泄及毒性(ADMET)性质不佳而被中止开发,给制药企业造成巨大的经济损失.目前很多公司运用高通量筛选的方法和虚拟ADMET预测的方法以降低药物研发后期的风险性.具有HERG阻滞活性的部分药物,其化学结构具有多样性,如能借助计算机模拟技术对新化学实体是否具有致心律失常的不良反应进行有效的预测和筛选,这将为新药的设计和研究中避免心律失常的不良反应产生指导作用.近年来,有一些关于药物对HERG钾通道抑制活性的虚拟预测方法的报道,例如定量构效关系(QSAR)[3-7]和分类方法[1,8],包括朴素贝叶斯(naive bayes)[9]、决策树(decision tree)[10]、随机森林(random forests)[11]、偏最小二乘(PLS)[11,12]和SVM[12-15]等方法.然而,这些研究大多使用结构差异性小的数据组构建的分类模型,其应用范围有限;有的模型虽用了结构差异性大的数据组,结果却不很理想.因此,为了提高模型的预测精度,拓宽模型的应用范围,有必要详细研究如何使用结构差异性大的数据组来构建分类预测模型.1原理与方法1.1数据收集所研究的化合物包含两个部分:训练集和外部测试集.训练集含367个结构多样性的化合物,其实验半数抑制浓度IC50数据来自文献[3],列于附表S1(见Supporting Information,在物理化学学报网站上可免费下载).该文献共提供了371个HERG钾通道抑制剂,其中有两个(附表S1中No.277(RP58866)和No.290(terikalant)[3])化合物为旋光异构体,由于本文的分子描述符中没有包含区别旋光异构体的描述符,所以本研究仅保留了其中一个,即terikalant.另外,附表S1中145号、238号和246号化合物[3] (分别为dofetilide N-acetyl,N-demethylerythromycin, noracetylmethadone),因不清楚其结构,所以本文数据集中也没有包含这三个化合物.在分类模型的建立中,需根据化合物IC50大小将数据集中化合物分类为有活性化合物和无活性化合物.本文采用与文献[3]相同的分类标准,即IC50=1.0、10.0μmol·L-1为阈值,构造了三个分类模型:模型I,按IC50=1.0μmol·L-1为标准,有118个化合物IC50<1.0μmol·L-1,被划分为活性化合物,有249个化合物IC50≥1.0μmol·L-1,被划分为非活性化合物;模型II,按IC50= 10.0μmol·L-1为标准,有229个化合物IC50<10.0μmol·L-1,被划分为活性化合物,有135个化合物IC50≥10.0μmol·L-1,被划分为非活性化合物;模型III,同时采用IC50=1.0μmol·L-1和IC50=10.0μmol·L-1为分类标准,将IC50<1.0μmol·L-1的化合物划分为活性化合物,而将IC50>10.0μmol·L-1的化合物划分为非活性化合物,这样,数据集中有118个化合物被划分为活性化合物,132个化合物作为被划分为非活性化合物.对每一个模型,分别使用五重交叉方法进行模型验证,即将所有数据随机地分成样本数目相等(或近似相等)的五组,每次取其中的四组用于训练模型,称为训练集,剩余的一组用于测试模型的预测精度,称为测试集,如此交替反复五次,使每一组依此用作测试集.外部测试集为从文献[16-19]中收集的97个新化合物(列于附表S2,见Supporting Information),用于衡量由训练集建立模型的预测能力.所有分子的稳定结构均用Hyperchem7.0[20]中的MM+分子力场优化得到.1.2分子描述符分子描述符用本小组自行编制的MODEL[21]软件计算.此软件可以计算3778个描述符,共有六类,分别为组成描述符、电荷相关描述符、物理化学性质描述符、拓扑描述符、几何描述符和量子化学描述符.为了节约计算时间和成本,本文未包含量子化学描述符.因而,对每个分子计算了1559个描述符,包括21个组成描述符、19个电荷相关描述符、2个物理化学性质描述符、471个拓扑描述符和1046个几何描述符.1.3描述符筛选方法由于所计算的1559个描述符并非都与分子的所研究的活性相关,且描述符之间可能存在线性相关,因此,有必要从中选择恰当的描述符子集以使所建立的模型具有最优的预测能力.本文采用一套混合描述符筛选方法来寻找最佳的描述符子集.其步骤如下.(1)预处理.首先,如果有90%样本的同一个描述符值是相同的,而这些样本的活性又具有显著区别,显然,这个描述符对区别分子的活性无效,将它删除.其次,如果对任意一个描述符,遍及所有的样1582No.8李平等：基于支持向量机方法的HERG钾离子通道抑制剂分类模型本计算该描述符的相对标准偏差,如果小于0.05,则删除该描述符.最后,如果两个描述符的Pearson相关系数[22]大于0.9,则说明这两个描述符高度相关,删除其中之一描述符.(2)描述符排序.经过预处理后的描述符,按照Fischer Score(F-Score)[23]值递减的顺序排序.F-Score 值(F(i))是单个描述符的分辨能力的简单度量,定义为,F(i)=(μ+i-μ-i)2(σ+i+σ-i)2(1)其中μ+i和μ-i分别是第i个描述符的正、负样本的平均值,σ+i和σ-i分别是它们的标准偏差.某个描述符的F-Score值越大,它的分辨能力就越强.筛选步骤如下:①计算每一个描述符的F-Score值;②将所有描述符按照F-Score值从大到小的顺序排序;③用五重交叉方法得到SVM模型平均推广误差,并优化模型参数σ,使SVM模型平均推广误差最小;④删除最后10个描述符,回到第③步,直至推广误差最小.由此,得到按F-Score值过滤的优化描述符子集及其对应的优化σ值.(3)Monte Carlo模拟退火算法.模拟退火算法是一种物理模拟过程,在此,运用蒙特卡洛模拟退火法(MC-SA)结合SVM,文献[24]有类似的描述.可以简要概括步骤如下:①给定高斯核函数的初值σ;②给定模拟初始温度T;③随机选择一组描述符子集作为初始解;④用五重交叉验证计算SVM模型的平均总预测精度作为适应度函数Q old;⑤对描述符子集给一随机微扰,得到一新描述符子集,按照上述方法获得新的适应度函数Q new;⑥如果Q old<Q new,则接受新的描述符子集替代原来的描述符子集;如果Q old>Q new,则以概率P=e-(Q old-Q new)/T接受为新描述符子集;⑦回到步骤⑤,直到达到此温度下的最大平衡步数;⑧降低模拟退火温度T,回到步骤③,直到达到最低模拟退火温度;⑨系统地调整σ值,返回步骤②,直至得到最大平均总的预测精度.由此,可得最佳描述符子集和对应的σ值,给出最终的SVM-MC模型.1.4SVM方法SVM是Vapnik等在20世纪90年代提出的一种学习算法,文献[25]中有详细介绍,在此仅简述.对于两类分类问题,给定一训练集{(x i,y i),i=1,2,…, n},其中y i∈{-1,+1},表示任一样本x i的分类标识,如果训练集是线性可分的,SVM就是寻求超平面,f(x)=w·x+b=0(2)其中,“·”是点积,w是n维向量,b为偏移量.使正样本(y i=+1)和负样本(y i=-1)可分,且使其边界上的点到该超平面的距离最大.这可以转化为以下条件w·x i+b≥+1(y i=+1)和w·x i+b≤-1(y i=-1)(3)限制下求函数ψ(w,b)=12||w||2的最小值.Lagrange 乘数法可得解w=iΣαi y i x i(4)αi是Lagrange乘子,(4)式满足限制条件iΣαi y i=0,并有最优分类决策函数f(x)=sgn(w·x+b)=sgniΣαi y i(x i·x)+ΣΣb(5)式(5)中sgn为符号函数.由于很多两类情形并非线性可分,为此,SVM将样本点x通过函数准(x)投影到高维空间以使其线性可分.但SVM并非直接引入准(x),而是通过核函数K(x i,x)方法间接引入的: K(x i,x)=准(x i)·准(x)(6)其分类决策函数变为:f(x)=sgniΣαi y i K(x i,x)+ΣΣb(7)核函数K(x i,x)可以有多种形式,本文使用Gaussian 核函数,即K(x i,x)=exp-||x i-x||2/2σΣΣ2(8) SVM模型参数σ的选择是通过最小化推广误差进行的.1.5结果评价采用广泛使用的预测正确率[26]来评估分类模型的预测能力,模型的敏感性SE=TP/(TP+FN)用于评估正样本的预测正确率,专一性SP=TN/(TN+FP)用于评估负样本的预测正确率,Q用于评估总样本的预测正确率,表达式如下:Q=TP+TN(9)其中,TP、TN、FP、FN分别代表预测正确的正样本数、预测正确的负样本数、预测错误的正样本数以及预测错误的负样本数.2结果与讨论2.1SVM模型的预测结果根据不同的分类标准,分别建立了模型I、模型1583Acta Phys.-Chim.Sin.,2009Vol.25II 和模型III 三个模型.首先对每个分子计算了1559个描述符以表征其结构.然后,经预处理、过滤筛选和模拟退火筛选,分别建立分类模型.变量筛选对模型预测精度的影响列于表1,最终筛选出的描述符列于附表S3(见Supporting Information).从表1可见,经预处理后,三个模型的描述符的数目分别从1559个减小到633、631和617个.结果表明:在最初计算的1559个描述符中,大量描述符信息含量低或与其它描述符高度相关,经预处理即可去掉.经过预处理,用五重交叉验证得到的三个模型,正样本的平均预测正确率分别为63.7%、84.7%、84.8%,负样本的平均预测正确率分别为88.3%、63.0%、84.0%,总的预测正确率分别为80.4%、76.6%、84.4%.经过描述符按F -Score 排序并过滤后,模型I 、II 、III 的描述符数目分别由633个减少到118个,631个减少到95个,617个减少到166个.正样本的预测正确率分别为72.9%、83.0%、89.9%;负样本的预测正确率分别为88.8%、72.6%、92.4%,总的预测正确率为83.7%、79.1%和91.2%.可见无论正、负样本预测正确率,还是总样本预测正确率都有所提高.结果表明,排序过程在描述符的选择中起了重要的作用,能大大减少描述符的数目,有利于提高模型建立的效率.然而,描述符的数目仍然很大,因此,采用Monte Carlo 模拟退火法进一步筛选描述符,最终模型I 、II 、III 的描述符数目分别分为42、25、40个,这三个模型的正样本的预测正确率分别为84.8%、90.8%、96.6%,负样本的预测正确率分别为92.8%、80.7%、97.7%,总样本的预测正确率分别是90.2%、87.1%、97.2%.研究结果表明,变量选择在简化模型的同时明显地提高了各模型的预测精度,最后的SVM -MC 模型具有良好的预测能力.对模型I 、II 、III 其五重交叉验证平均预测精度对正样本在84.8%-96.6%范围,负样本在80.7%-97.7%范围,总样本在87.1%-97.2%范围.而且,值得注意的是,模型III 具有优良的预测能力,其正、负样本的预测精度分别为96.6%和97.7%,优于其它文献报道的最好结果[3].同时,从这三个模型的所有化合物及所筛选的分子描述符出发,并固定以上交叉验证优化得到的模型参数σ,分别建立模型,并计算含97个化合物的外部测试集的预测精度,得到三个训练模型总的预测精度分别为67.0%、85.5%、90.1%.模型I 结果接近文献[14]结果,而模型II 、III 结果明显优于文献[14]结果.2.2与文献报道结果[3,14]比较表1变量筛选对SVM 性能的影响Table 1Effect of feature selection on the performance of SVMa)step 1:preprocessing,step 2:ranking and backward selection,step 3:Monte Carlo simulated annealing;b)optimized values for the exponent of the Gaussian kernel of SVM;c)CV set:cross -validation set;d)SE:sensitivity;e)SP:specificity;f)Q :overall prediction accuracyStep a Number of descriptor(σb )CV set cPrediction accuracySE d (%)SP e (%)Q f (%)IIIIIIIII III I II III I II III 1633(19.0)631(13.0)617(20.0)175.087.091.792.074.188.986.582.290.2258.389.179.294.063.096.382.479.588.2345.887.075.092.066.780.877.079.578.0469.674.087.080.063.073.176.769.979.6569.686.791.383.748.280.879.272.285.7average63.784.784.888.363.084.080.476.684.42118(8.5)95(8.5)166(10.0)179.289.1100.094.074.192.689.283.696.1275.080.483.386.070.4100.082.476.792.2362.584.875.090.085.288.581.184.982.0469.678.395.780.077.892.376.778.193.9578.382.295.793.955.688.588.972.291.8average72.983.089.988.872.692.483.779.191.2342(6.0)25(3.5)40(9.0)187.593.595.896.085.2100.093.290.498.0287.587.095.890.081.596.389.284.996.1370.893.595.896.081.5100.087.889.098.0482.689.195.786.085.2100.084.987.798.0595.791.1100.095.970.492.395.883.395.9average 84.890.896.692.880.797.790.287.197.21584No.8李平等：基于支持向量机方法的HERG 钾离子通道抑制剂分类模型将本文建立的SVM -MC 模型,与文献报道的采用其它机器学习方法建立的模型所得的预测结果进行比较,有助于检测本研究的模型精度是否与其它模型的精度处于相似的水平.但是,值得注意的是,由于使用不同的样本、分子描述符、分类方法和参数以及验证的方法,将不同模型的预测结果进行直接比较是不恰当的.无论如何,可以进行尝试性比较,以便对HERG 钾通道抑制剂预测精度的大致水平,进行粗略的评价.本研究结果与文献结果[3,14]比较列于表2.由表2可知,本研究的SVM -MC 模型的预测结果无论是交叉验证结果还是外部测试结果,都接近或优于文献[3,14]报道结果,尤其是模型III,明显优于文献[3,14]报道的结果.因此,本研究的SVM -MC 模型有助于HERG 钾通道抑制剂先导化合物的筛选.3结论在一个差异性较大的HERG 钾通道抑制剂数据集的基础上,从一个大的分子描述符集合出发,经描述符的预处理、按F -Score 排序过滤和Monte Carlo 模拟退火法的筛选,建立了按不同IC 50为分类标准的分类预测模型.结果表明,在去除多余的描述符和提高SVM 模型的预测能力方面,本文所提出的描述符选择方法是很有效的,在简化模型的同时模型的预测能力可显著提高.通过本研究建立的SVM -MC 模型与文献报道结果比较表明,SVM -MC 模型用于预测化合物的HERG 钾通道抑制活性有较好的预测能力,从而可用于从化合物数据库进行虚拟筛选,有助于发现潜在的HERG 钾通道抑制剂.总之,采用SVM 构建了较高精度的HERG 钾通道抑制剂鉴别模型.Supporting Information Available:Table S1listed 367HERG potassium channel inhibitors in the training set with their experimental HERG IC 50value (μmol ·L -1),therapeutic area,and cell line.Table S2listed 97compounds in the external testing set with their experimental HERG IC 50value (μmol ·L -1).The selected molecular descriptors by the three models were listed in Table S3.Table S4listed 1559molecular descriptors.This information is available free of charge via the internet at .References1Roche,O.;Trube,G.;Zuegge,J.;Pflimlin,P.;Alanine,A.;Schneider,G.ChemBioChem ,2002,3:4552Trudeau,M.C.;Warmke,J.W.;Ganetzky,B.;Robertson,G.A.Science ,1995,269:923Thai,K.M.;Ecker,G.F.Bioorg.Med.Chem .,2008,16:41074Cianchetta,G.;Li,Y.;Kang,J.;Rampe,D.;Fravolini,A.;Cruciani,G.;Vaz,R.J.Bioorg.Med.Chem.Lett .,2005,15:36375Coi,A.;Massarelli,I.;Murgia,L.;Saraceno,M.;Calderone,V.;Bianucci,A.M.Bioorg.Med.Chem .,2006,14:31536Aptula,A.O.;Cronin,M.T.SAR QSAR Environ.Res .,2004,15:3997Keseru,G.M.Bioorg.Med.Chem.Lett .,2003,13:27738Chen,X.;Li,H.;Yap,C.W.;Ung,C.Y.;Jiang,L.;Cao,Z.W.;Li,Y.X.;Chen,Y.Z.Cardiovasc.Hematol.Agents Med.Chem .,2007,5:119Sun,H.M.ChemMedChem ,2006,1:31510Gepp,M.M.;Hutter,M.C.Bioorg.Med.Chem .,2006,14:532511Gavaghan,C.L.;Arnby,C.H.;Blomberg,N.;Strandlund,G.;Boyer,put.Aided Mol.Des .,2007,21:18912Song,M.;Clark,M.J.Chem.Inf.Model .,2006,46:39213Yap,C.W.;Cai,C.Z.;Xue,Y.;Chen,Y.Z.Toxicol.Sci .,2004,79:17014Tobita,M.;Nishikawa,T.;Nagashima,R.Bioorg.Med.Chem.a)10-fold cross -validation;b)LOO cross -validation;c)5-fold cross -validation;d)TdP -and non -TdP -causing agent,unavailable IC 50values;e)result for compounds with IC 50<1.0μmol ·L -1and IC 50>10.0μmol ·L -1only表2文献[3,14]和本文中HERG 抑制剂分类模型的预测精度比较Table 2Comparison of the prediction accuracies of different classification models for HERGinhibitors from literature [3,14]and this workStudy Method Threshold (μmol ·L -1)Training setExternal testing Compd.SE(%)SP(%)Q (%)Compd.Q (%)Ref.[14]SVM 40.07397.087.095.0a 827d 67.0-78.01.07386.093.090.0a827d67.0-78.0Ref.[3]Bayes1.022370.093.085.0b5884.010.022390.068.083.0b 5878.0active(<1.0),inactive(>10.0)15095.080.087.0b,e5886.0this work SVM -MC 1.036784.892.890.2c9667.010.036490.880.787.1c9785.5active(<1.0),inactive(>10.0)25096.697.797.2c,e7190.11585Acta Phys.-Chim.Sin.,2009Vol.25Lett.,2005,15:288615Leong,M.K.Chem.Res.Toxicol.,2007,20:21716Duncton,M.A.J.;Chekler,E.L.P.;Katoch-Rouse,R.;Sherman,D.;Wong,W.C.;Smith,L.M.;Kawakami,J.K.;Kiselyov,A.S.;Milligan,D.L.;Balagtas,C.;Hadari,Y.R.;Wang,Y.;Patel,S.N.;Rolster,R.L.;Tonra,J.R.;Surguladze,D.;Mitelman,S.;Kussie,P.;Bohlen,P.;Doody,J.F.Bioorg.Med.Chem.,2009,17:73117Kawai,M.;Ando,K.;Matsumoto,Y.;Sakurada,I.;Hirota,M.;Nakamura,H.;Ohta,A.;Sudo,M.;Hattori,K.;Takashima,T.;Hizue,M.;Watanabe,S.;Fujita,I.;Mizutani,M.;Kawamura,M.Bioorg.Med.Chem.Lett.,2007,17:555818Berglund,S.;Egner,B.J.;Gradén,H.;Gradén,J.;Morgan,D.G.A.;Inghardt,T.;Giordanetto,F.Bioorg.Med.Chem.Lett.,2008,18:485919Guay,D.;Boulet,L.;Friesen,R.W.;Girard,M.;Hamel,P.;Huang, Z.;Laliberté,F.;Laliberté,S.;Mancini,J.A.;Muise,E.;Pon,D.;Styhler,A.Bioorg.Med.Chem.Lett.,2008,18:555420Hyperchem release7.0(Beta1.04),software available at http:// 21Li,Z.R.;Han,L.Y.;Xue,Y.;Yap,C.W.;Li,H.;Jiang,L.;Chen, Y.Z.Biotechnol.Bioeng.,2007,97:38922Yang,C.W.;Su,J.Y.;Tsou,A.P.;Chau,G.Y.;Liu,H.L.;Chen,C.H.;Chien,C.Y.;Chou,mun.,2005,330:48923Fisher,R.A.Annals of Eugenics,1936,7:17924Ajmani,S.;Jadhav,K.;Kulkarni,S.A.J.Chem.Inf.Model.,2006, 46:2425Czerminski,R.;Yasri,A.;Hartsourgh,D.Quant.Struct.-Act.Relat.,2001,20:22726Bhasin,M.;Zhang,H.;Reinherz,E.L.;Reche,P.A.FEBS Lett., 2005,579:43021586。