技术领域
[0001] 本
发明涉及一种CNS药物关键特征识别方法，属于计算机辅助药物设计领域。
背景技术
[0002] 目前，全世界数亿人受到中枢神经系统(Central Nervous System，CNS)
疾病的影响。由于大脑环境的特殊性，相关药物的研发具有成功率低、成本高、周期长等弊端，开发新的CNS药物迫在眉睫。设计可靠的CNS候选药物可以大大缩减新药研发的周期和成本，并显著提高成功率。而了解CNS药物与non-CNS药物的特征差异是设计有效CNS候选药物的前提。因此，发现CNS药物中的关键特征有助于我们了解CNS药物的特殊性，并指导CNS药物设计。
[0003] 对于如何从CNS药物的大量特征中筛选出关键特征，Shahid M(《SVM Based Descriptor Selection and Classification of Neurodegenerative Disease Drugs for Pharmacological Modeling》.Molecular Informatics,2013,32(3):241-249.)等人利用
支持向量机构建模型，通过每个特征的系数计算特征评分来排序特征，也可以用来进行特征选择。但根据单个特征的评分来删除不重要的特征，忽略了特征间组合的效果，有些特征单独作用不大，但是两个不重要的特征组合起来可能作用很大。Lu J(《Analysis of a drug target-based classification system using molecular descriptors》.Combinatorial chemistry&high throughput screening,2016,19(2):129-135.)等人从0开始逐个增加特征；但这种方法在最初单个特征时，含有的信息量很少，很可能会出现SEN＝0％，SPE＝100％或者SEN＝100％，SPE＝0％的情况，那这种情况下该选择留下哪个特征无法度量，而最开始选择的特征对后续的特征组合的预测性能有很大影响；若是从多个特征开始使用IFS
算法，那最开始的多个特征又需要采用其他方法确定。
[0004] 因此准确找到CNS药物与non-CNS药物小分子间的关键特征对帮助人们设计CNS药物分子，开发新的CNS药物有着很大的作用。
发明内容
[0005] 为了找到CNS药物与non-CNS药物小分子间的关键特征，从而达到可以指导CNS药物设计的目的，本发明提供了一种针对CNS药物关键特征的识别方法，所述方法将支持向量机和贪心算法相结合，利用贪心算法逐步删除对提升预测结果作用最小的特征，进而准确筛选出区分CNS药物与non-CNS药物小分子的关键特征。
[0006] 可选的，所述方法包括：
[0007] 步骤一、从CNS药物与non-CNS药物小分子的所有特征中初步筛选出对区分二者有作用的特征；
[0008] 步骤二、利用步骤一初步筛选出的对区分CNS药物与non-CNS药物有作用的特征构建支持向量机模型并优化参数c和g，得到优化后的支持向量机模型；
[0009] 步骤三、利用贪心算法逐步删除步骤一初步筛选出的对区分CNS药物与non-CNS药物有作用的特征，从删除过程中筛选出区分CNS药物与non-CNS药物的关键特征。
[0010] 可选的，假设步骤一初步筛选出的对区分CNS药物与non-CNS药物有作用的特征的个数为n；则所述步骤三包括：
[0011] 3.1逐个删除每一个特征，得到n个不同的特征组合：{a2，a3，a4，…an}，{a1，a3，a4，…an}，{a1，a2，a4，…an}，…{a1，a2，a3，a4，…an-1}；
[0012] 3.2将所述n个不同的特征组合作为步骤二得到的优化后的支持向量机模型的输入向量，得到n个不同的特征组合分别对应的预测性能，并保留预测性能最好的一个特征组合；
[0013] 3.3以3.2得到的预测性能最好的一个特征组合中的n-1个特征执行3.1至3.2，以此循环直至n个特征被删除完毕；
[0014] 3.4从上述3.1至3.3执行过程中选择预测性能最好的一个特征组合，该特征组合中的特征即为区分CNS药物与non-CNS药物的关键特征。
[0015] 可选的，所述预测性能包括灵敏度SEN和特异性SPE；SEN表示CNS药物的预测率，SPE表示non-CNS药物的预测率。
[0016] 可选的，所述3.2中保留预测性能最好的一个特征组合，包括：
[0017] 分别比较各个特征组合对应的SEN值和SPE值，选出最高的SEN值和SPE值；
[0018] 若最高的SEN和SPE属于同一个特征组合，则保留这个特征组合；
[0019] 若最高的SEN和SPE属于两个不同的特征组合，则根据两个不同的特征组合各自的SEN和SPE综合确定所要保留的特征组合。
[0020] 可选的，假设最高的SEN和SPE分别属于两个不同的特征组合A和B，则所述根据两个不同的特征组合各自的SEN和SPE综合确定所要保留的特征组合，包括：
[0021] 比较特征组合A的SPE和特征组合B的SEN；
[0022] 若特征组合A的SPE大于特征组合B的SEN，则选择保留特征组合A；
[0023] 若特征组合A的SPE小于特征组合B的SEN，则选择保留特征组合B；
[0024] 若特征组合A的SPE等于特征组合B的SEN，则比较特征组合A的SEN和特征组合B的SPE的大小，选择较大者对应的特征组合。
[0025] 可选的，若两个特征组合的SPE相等、SEN也相等，则随机保留特征组合A或者特征组合B。
[0026] 可选的，所述步骤一从CNS药物与non-CNS药物小分子的所有特征中初步筛选出对区分二者有作用的特征，采用
随机森林算法，使用信息增益率作为属性划分评价函数，进行初步特征选择。
[0027] 可选的，所述步骤二采用穷举法得到优化后的支持向量机模型。
[0028] 本发明还提供一种CNS药物分子设计方法，所述设计方法采用上述方法识别CNS药物的关键特征。
[0029] 本发明有益效果是：
[0030] 通过将支持向量机和贪心算法相结合，利用贪心的思想逐步删除掉对提升预测结果作用最小的特征，进而准确筛选出区分CNS药物与non-CNS药物小分子的关键特征。本发明方法首次将支持向量机与贪心算法结合起来应用于CNS药物的关键特征识别，以逐步删除的方式筛选出关键特征考虑了各个特征之间组合的效果，且避免了增加特征的方法所带来的最初特征选择的难题，使得所筛选出的关键特征更能有效区分CNS药物与non-CNS药物小分子，这为从根本上设计CNS药物候选小分子提供了一种重要的指导方法。
具体实施方式
[0031] 为使本发明的目的、技术方案和优点更加清楚，下面将对本发明实施方式作进一步地详细描述。
[0033] 本实施例提供一种基于支持向量机和贪心算法的CNS药物关键特征的识别方法，所述方法将支持向量机和贪心算法相结合，利用贪心的思想逐步删除掉对提升预测结果作用最小的特征，进而准确筛选出区分CNS药物与non-CNS药物小分子的关键特征，包括：
[0034] 步骤(1)采用随机森林算法进行初步特征选择：
[0035] 构建随机森林模型，使用信息增益率作为属性划分评价函数，进行初步特征选择；
[0036] 具体的，构建包含100棵
决策树的随机森林模型，可参考“姚登举,杨静,詹晓娟.基于随机森林的特征选择算法[J].吉林大学学报(工学版)(01).”，使用信息增益率作为属性划分评价函数。为了提高选择效率及性能，每次随机选择2/3的样本及1/2特征来构建决策树；为防止过拟合，当未分配样本数小于5时
节点中止分裂。
[0037] 统计所有出现在树上的特征，即为初步选择出的特征。
[0038] 步骤(2)采用穷举法优化支持向量机模型参数c和g，c为惩罚系数，g为核参数；
[0039] 使用LIBSVM包中带有径向基核函数的支持向量机算法来构建识别CNS药物与non-CNS药物小分子的分类器。
[0040] 在[2-4，24]范围内穷举所有的c和g的组合，在每种组合下进行5折交叉验证，找出最优的c和g的组合。
[0041] 支持向量机的目标优化问题就是找到一个可以将CNS样本与non-CNS样本尽可能区分的
超平面，公式如下：
[0042] f(x)＝wTx+b
[0043] 其中，x为输入
特征向量，w为超平面的垂直法向量，b为偏置。
[0044] 根据拉格朗日乘子，求得w；采用映射函数Φ将特征向量xi和x映射至高维空间，具体如下式所示：
[0045]
[0046]
[0047] 其中，λi是拉格朗日乘子，yi是与特征向量xi相关的样本标签，m是样本个数，1≤i≤m。
[0048] 如果没有核函数，高维空间的计算很可能会导致维度爆炸，为了避免这个问题，径向基核函数K(xi,x)被用来替代显式映射ΦT(xi)Φ(x)，如下所示：
[0049]
[0050] K(xi,x)＝exp(-g*|xi-x|2)
[0051] 其中，g参数是非常重要的核参数，对模型的训练具有很大影响。另一个重要的参数是惩罚系数c，它影响着分类平面的平滑程度。
[0052] 将步骤(1)初步选择出的特征作为输入向量计算所有的c和g的组合对应的支持向量机模型的预测性能，选出预测性能最好的一组对应的c和g，以该组c和g作为优化后的支持向量机模型参数c和g。
[0053] 步骤(3)利用贪心算法识别出关键特征
[0054] 将步骤(1)初步选择出的特征的不同组合分别作为步骤(2)优化后的支持向量机模型的输入向量，根据对应的预测性能筛选出关键特征。
[0055] 具体的，包括以下步骤：
[0056] S1假设上述步骤(1)初步选择出n个特征：{a1，a2，a3，a4，…an}；
[0057] S2逐个删除每一个特征，得到n个不同的特征组合：{a2，a3，a4，…an}，{a1，a3，a4，…an}，{a1，a2，a4，…an}，…{a1，a2，a3，a4，…an-1}；
[0058] S3将步骤S2中每个特征组合作为支持向量机模型的输入向量，保留预测性能最好的特征组合，并记录预测性能pj，1≤j≤n；
[0059] S4使用步骤S3中保留下来的一组特征组合，以该组中的n-1个特征执行S2和S3，以此循环执行S2至S4直到所有特征被删完；
[0060] S5上述S2至S4过程中所有预测性能pj中最好的pj所对应的特征组合，即筛选出的关键特征。
[0061] 上述步骤(2)和步骤(3)中，预测性能pj包括灵敏度SEN和特异性SPE：
[0062] 灵敏度SEN，即正样本预测率(CNS药物的预测率)；
[0063] 特异性SPE，即负样本预测率(non-CNS药物的预测率)；
[0064] 通过灵敏度SEN和特异性SPE共同评估构建的CNS药物识别模型，值越大表示模型的性能越好。
[0065] 具体的，在确定预测性能pj中最好的pj时，若最高的SEN和SPE属于同一个特征组合，则保留这个特征组合。
[0066] 若最高的SEN和SPE属于不同的特征组合，则根据其各自对应的SPE和SEN的大小选择所要保留的特征组合；比如，最高的SEN和SPE分别属于特征组合A和B，即特征组合A的SEN最高，而特征组合B的SPE最高，则比较特征组合A的SPE和特征组合B的SEN：
[0067] 若特征组合A的SPE大于特征组合B的SEN，则选择保留特征组合A；
[0068] 若特征组合A的SPE小于特征组合B的SEN，则选择保留特征组合B；
[0069] 若特征组合A的SPE等于特征组合B的SEN，则比较特征组合A的SEN和特征组合B的SPE的大小，选择较大者对应的特征组合。
[0070] 若两个特征组合的SPE相等、SEN也相等，则随机保留特征组合A或者特征组合B。
[0071] 为验证本
申请提出的方法的能够有效识别出CNS药物的关键特征，本申请以现存的CNS药物和non-CNS药物小分子为实验对象，数据来源于ZINC15(http://zinc15.docking.org/)和DrugBank(https://www.drugbank.ca/)
数据库。
[0072]
发明人从上述两个数据库下载SDF格式的药物数据(对应于879个non-CNS药物小分子和273个CNS药物小分子)，包括每个药物分子中所有
原子的初始坐标和原子间的键类型信息；将该药物数据作为PaDEL
软件的输入，得出每个药物分子的所有特征值，本实施例中计算出的每个药物分子的所有特征值个数为1875个。
[0073] 步骤(1)采用随机森林算法进行初步特征选择：
[0074] 构建随机森林模型，选择对CNS药物和non-CNS药物小分子的区分有贡献的特征，从1875个特征中选择了941个有用的特征。
[0075] 步骤(2)采用穷举法优化支持向量机模型参数c和g，得到优化后的支持向量机模型参数c和g：
[0076] 在879个non-CNS药物小分子中随机选择355个作为负样本，以273个CNS药物小分子作为正样本，将支持向量机的参数c和g置于[2-4，24]的范围内，在所有c和g组合中搜索最佳的5折交叉验证结果，并使用测试集测试模型泛化性能。重复上述过程5次，结果如下表1所示，5个样本中正样本CNS药物相同，负样本分别为从879个non-CNS药物小分子中随机选择出的355个，也即5个样本中正样本CNS药物相同，负样本对应不同的355个non-CNS药物。
[0077] 将每个药物小分子的941个特征作为不同c和g组合对应的支持向量机模型的输入向量对其进行性能预测。
[0078] 表1.支持向量机模型在外部测试集上的参数及性能
[0079]
[0080] 根据表1可知，样本3对应的参数c和g所对应的支持向量机模型的预测性能最优，因此以此该组c和g对应的支持向量机模型进行关键特征的筛选，并且以该组对应的随机选出的355个non-CNS药物小分子以及273个CNS药物小分子作为筛选关键特征的样本。
[0081] 步骤(3)利用贪心算法识别出关键特征：
[0082] 将上述步骤(2)确定的样本3对应的随机选出的355个non-CNS药物小分子以及273个CNS药物小分子对应的941个特征作为优化后的支持向量机模型的输入向量进行关键特征的选择，具体的：
[0083] S1：941个特征为：{a1，a2，a3，a4，…an}，n＝941；
[0084] S2：逐个删除每一个特征，得到941个不同的特征组合：{a2，a3，a4，…an}，{a1，a3，a4，…an}，{a1，a2，a4，…an}，…{a1，a2，a3，a4，…an-1}；
[0085] S3：将步骤S2中每个特征组合作为支持向量机模型的输入向量，保留预测性能最好的特征组合，并记录各组特征组合的灵敏度SEN和特异性SPE的值；
[0086] 将941个特征组合各自的灵敏度SEN和特异性SPE值分别进行比较，若最高的SEN和SPE属于同一个特征组合，则保留该特征组合。
[0087] 若最高的SEN和SPE属于不同的特征组合，则根据其各自对应的SPE和SEN的大小选择所要保留的特征组合；比如，最高的SPE属于第800个特征组合，而最高的SEN对应于第900个特征组合，则将第800个特征组合的SEN与第900个特征组合SPE进行比较：
[0088] 若第800个特征组合的SEN大于第900个特征组合SPE，则保留第800个特征组合；
[0089] 若第800个特征组合的SEN小于第900个特征组合SPE，则保留第900个特征组合；
[0090] 若第800个特征组合的SEN等于第900个特征组合SPE，则比较第800个特征组合的SPE与第900个特征组合SEN的大小：
[0091] ——若第800个特征组合的SPE大于第900个特征组合SEN，则保留第800个特征组合；
[0092] ——若第800个特征组合的SPE小于第900个特征组合SEN，则保留第900个特征组合；
[0093] ——若第800个特征组合的SPE等于第900个特征组合SEN，则随机保留第800个特征组合或第900个特征组合。
[0094] S4：使用步骤S3中保留下来的一组特征组合，以该组中的940个特征执行S2和S3，以此循环执行S2至S4直到所有特征被删完；
[0095] S5：上述S2至S4过程中所有预测性能pj中最好的pj所对应的特征组合，即筛选出的关键特征。
[0096] 通过上述循环过程，本实施例最终筛选了40个关键特征，使用筛选出的40个关键特征作为模型的输入变量，测试结果的SEN和SPE都达到了94％以上。所筛选出的40个关键特征如下表2所示：
[0097] 表2.筛选出的40个关键特征及其说明
[0098]
[0099]
[0100]
[0101] 本发明实施例中的部分步骤，可以利用软件实现，相应的软件程序可以存储在可读取的存储介质中，如光盘或
硬盘等。
[0102] 以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何
修改、等同替换、改进等，均应包含在本发明的保护范围之内。
