一种基于作用模式的人类发育毒性预测的方法
技术领域
本发明涉及化学品的人类发育毒性虚拟筛选与活性预测领域，更具体的说，涉及一种基于作用模式的人类发育毒性预测的方法。
背景技术
大量动物实验及流行病学研究结果发现，多数环境小分子化合物具有干扰核受体介导的核酸翻译和表达功能因而影响人类个体生长发育过程，称为发育毒性(Developmental Toxicity)。特别地，以美国环保署(U.S.Environmental ProtectionAgency,U.S.EPA)和欧洲化学品管理局(European Chemicals Agency,ECHA)在内的多个官方和非官方组织发现许多环境相关污染物会潜在影响有害结局路径(Adverse OutcomePathway,AOP)并导致人类生殖相关器官(男性和女性)的发育异常。例如。对于男性而言，污染物会通过影响雄激素受体介导有害结局路径并导致男性生殖相关器官的发育异常，包括睾丸、前列腺的生长迟缓、功能不全或异常等有害效应，进而产生个体水平的男性发育毒性(Male Developmental Toxicity,MDT)。具有MDT的环境污染物不仅会导致男性生殖相关器官的发育异常，还会潜在导致男性生殖健康异常，包括增加睾丸生殖细胞肿瘤、低精液质量、隐睾和尿道下裂的发生率，最终产生生殖毒性(Reproductive Toxicity)。对于女性而言，污染物会通过影响雌激素受体介导有害结局路径并导致女性生殖相关器官的发育异常，包括卵巢、输卵管、子宫、胎盘和乳腺等，进而产生个体水平的女性发育毒性(FemaleDevelopmental Toxicity,FDT)。女性进行成功的怀孕生殖需要这些女性生殖相关器官的正常发育与运作，若功能失调会导致难以怀孕或无法怀孕，无法顺利地怀孕到足月，或难以喂养婴儿。
原则上，基于作用模式，环境污染物通过干扰有害结局路径，从分子水平、细胞水平、及器官组织产生人类发育毒性干扰，导致这些人类生殖相关器官发育不完全，或者发育畸形导致人类无法进行正常发育与运作，最终会导致人类生殖相关的功能障碍。目前，对发育相关的危害识别主要基于大量的动物实验以及流行病学研究结果，和少量的干扰机制研究。商业上有成千上万种化学品，但很少有进行过个体相关发育毒性测试，更少有针对人类生殖相关器官的发育异常为终点的流行病学研究。因此，近二十年来开展了大量的传统动物实验对化学品的人类发育毒性进行毒性测试。然而，自2004起，欧盟因伦理问题将传统动物测试禁止之后，基于体外测试和虚拟筛查的非动物替代测试(non-animal testing)方法应运而生。体外测试耗时常，成本高，无法完全测试数以万计的已登记在册的化学品，因此虚拟筛查技术显得尤为重要。
科学家们发展了基于计算机的虚拟筛选方法来进行化学品相关毒性终点的活性预测。定量结构效应关系(Quantitative Structure-Activity Relationship，QSAR)可以利用分子描述符提取并描绘化合物生物活性与结构特征之间的关系。QSAR作为一种成熟的方法已在多种毒性预测中广泛使用。例如发明创造名称：一种人运甲状腺素蛋白干扰物虚拟筛选方法(专利公开号：CN106407665A，公开日：2017-02-15)曾利用过QSAR技术构建了一种干扰物筛选方法。还有发明创造名称：一种基于金属定量构效关系的淡水急性基准预测方法(专利公开号：CN104820873A，公开日：2015-08-05)和发明创造名称：基于金属定量构效关系的海水急性基准预测方法(专利公开号：CN105447248A，公开日：2015-11-24)也利用QSAR技术来对海水及淡水急性基准进行预测。值得注意的是，存在许多非QSAR技术在发育毒性预测中进行了运用。例如发明创造名称：一种利用黑腹果蝇评价三唑类农药的生长发育毒性的方法(专利公开号：CN110150236A，公开日：2021-04-06)。然而，该方法不仅适用范围低(局限于三唑类农药)，而且预测发育终点为黑腹果蝇，不是高关注的人类健康毒性终点。还有发明创造名称：一种化学品发育毒性预测方法、预测模型及其构建方法与应用(专利公开号：CN112063681A，公开日：2020-12-11)。虽然该方法运用了人类心肌细胞数据进行预测，但是预测终点为人类心肌细胞中α-Actinin和SOX17蛋白活性，局限于细胞水平的活性预测，对器官、个体水平的发育毒性无法进行预测。
有害结局路径(Adverse outcome pathway,AOP)是用以描述已有的关于一个直接的分子起始事件(molecular initiating event,MIE)(如：配体-受体结合)与在生物不同组织结构层次(如:细胞、器官、机体、群体)所出现的与危险度评定相关的“有害结局”之间的相互联系。一条AOP的建立不仅确定了毒性发生过程中各个独立的毒性事件，以及毒性事件发生间的前后关系，从而将这些毒性事件模块化。AOP最终将细胞水平、器官水平以及个体水平的毒性事件相结合，综合评估化学品基于有害结局路径的毒性效应，并可预测化学品的毒性作用机制(如：干扰某个关键事件)。但通过对现有技术进行分析，现有技术中，缺少基于作用模式的人类发育毒性高通量预测的方法。
发明内容
技术问题：本发明的目的在于克服现有技术中，不能有效地基于作用模式进行人类发育毒性高通量预测的不足，提供一种基于作用模式的人类发育毒性预测的方法。该方法可以对潜在作用于有害结局路径而进一步产生发育毒性的化学品进行高通量筛查，从而可以准确快速的判断出化合物是否存在发育毒性，且该毒性是否是通过干扰有害结局通路而产生。
技术方案：本发明提供一种基于作用模式的人类发育毒性预测的方法，包括：
构建化合物活性数据集，包括基于作用模式，选择有害结局路径，确定第一数量的研究事件，所述第一数量的研究事件能够描述第二数量的毒性终点；收集并整理每个研究事件的相关试验及其或化合物活性数据；
根据化合物活性数据集，结合第三数量的分子描述符库，针对每一毒性终点，通过第四数量的不同种机器学习算法，分别训练出第四数量的对应的QSAR模型，并筛选出每一毒性终点对应的预测效果最佳的QSAR模型，获得的第二数量的QSAR模型的集合形成第一预测模型；
利用若干种具有体内实验数据的化合物，以及利用所述第一预测模型对具有体内实验数据的化合物进行预测的预测结果，利用朴素贝叶斯算法训练得到第二预测模型；
将待测化合物输入到所述第一预测模型中进行定性预测，并将定性预测结果输入所述第二预测模型，完成对化学品的人类发育毒性进行预测。
优选地，所述基于作用模式，选择有害结局路径，确定第一数量的研究事件的方式为：
基于有害结局路径的人类发育毒性机制，从分子水平、细胞水平、个体或器官水平收集包括分子启动事件和关键事件的第一数量的研究事件；
所述收集每个研究事件的相关试验的方式为：除基于动物实验的个体/器官水平研究事件外，从分子启动事件开始的所有研究事件数据均为高通量化学信息、体外试验；并以个体或器官水平研究事件作为研究有害结局或最后一个研究事件，所使用的动物试验需满足美国环保署或经济合作与发展组织的测试指南；
所述整理每个研究事件的活性数据的方式为：
首先，去除无结构信息、高分子型、离子型、和混合型化合物；
然后，采用如下公式将化合物活性进行标准化：
公式中，Activity value代表活性强度值，Ki代表抑制常数，Kd代表离解常数，AC50代表半数活性浓度，IC50代表半数抑制浓度，EC50代表半数效应浓度，uM表示微摩尔量；
最后，利用细胞毒性实验作为活性过滤器，去除潜在由于细胞毒性而导致的假阳性化合物。
优选地，活性数据整理后，针对每个研究事件，对收集的化合物进行活性分类，包括：
活性：针对某一研究事件，至少有一个实验存在活性，同时该活性若为抗性体外实验数据，其活性强度需大于同实验下的细胞毒性活性强度；
非活性：针对某一研究事件，所有实验测定结果都为无活性，或仅有抗性体外实验数据，且其抗性活性强度小于等于同实验下的细胞毒性活性强度；
拟性：针对某一研究事件，化合物存在活性前提下，存在拟性活性数据；
抗性：针对某一研究事件，化合物存在活性前提下，存在抗性活性数据，且抗性活性强度大于细胞毒性活性强度。
优选地，所述根据化合物活性数据集，结合第三数量的分子描述符库，针对每一毒性终点，通过第四数量的不同种机器学习算法，分别训练出第四数量的对应的QSAR模型，并筛选出每一毒性终点对应的预测效果最佳的QSAR模型，获得第一预测模型的方法包括：
将所构建的化合物活性数据集以4:1的比例分为训练集和测试集，其中训练集用作模型构建和内部验证，测试集用作外部验证；
选择第三数量的分子描述符库进行化合物的结构信息数据计算；
对于每个毒性终点，通过第四数量的不同种机器学习算法，分别训练得到第四数量的QSAR模型；
采用第五数量的不同指标，对QSAR模型的预测效果进行验证，对每个毒性终点，选择具有最佳预测能力的QSAR模型；
筛选出的第二数量的QSAR模型的集合，构成所述第一预测模型。
优选地，第四数量的所述机器学习算法包括K近邻算法、朴素贝叶斯算法、随机森林、支持向量机和决策树。
优选地，第三数量的所述分子描述符库分别为：OEState、Mold2和Dragon v.7。
优选地，第五数量所述的指标包括真阳性、假阳性、真阴性、假阴性、敏感度、特异度、精确度和曲线下面积。
优选地，在构建QSAR模型时，利用5折交叉验证对构建的QSAR模型的训练集进行内部验证，测试数据的稳定性。
优选地，所述体内实验数据为动物实验数据，部分所述动物实验数据能够检测化合物在器官/个体水平上表现出的多种毒性效应，使得所述第二预测模型包括多个独立的预测模型，每个模型单独预测一种毒性效应，当待测化合物经过多个第二预测模型中多个独立的预测模型预测后，若至少在一个模型预测中呈现阳性，则化合物存在人类发育毒性。
进一步地，还包括构建应用域，所述构建应用域的方法为：
利用训练集中化合物的26种物理化学性质来描述相关QSAR模型的应用域，26种所述物理化学性质包括1D、2D、3D，分别为nAcid，ALogP，AMR，apol，naAromAtom，nAromBond，nAtom，nHeavyAtom，nBonds，nBondsD，nBondsT，nBondsQ，bpol，ETA_Alpha，FMF，nHBAcc，nHBDon，TopoPSA，VABC，MW，AMW，XLogP，TPSA，nRing，nRotB和nRotBt；
采用欧式距离计算出训练集的距离矩阵作为每个QSAR模型的应用域。
有益效果：本发明与现有技术相比，具有以下优点：
本发明的实施例中所提供的方法，基于作用模式，选择有害结局路径，确定研究事件，收集并整理每个研究事件的相关试验及其化合物活性数据，进行化合物活性数据集的构建；然后，利于构建的化合物活性数据集，并结合分子描述符库以及机器学习算法，训练得到第一预测模型；然后利用若干种具有体内实验数据的化合物，以及利用所述第一预测模型对具有体内实验数据的化合物进行预测的预测结果，利用朴素贝叶斯算法训练得到第二预测模型；将待测化合物通过第一预测模型和第二预测模型两阶段的预测验证，从而准确快速地判断出化合物是否存在发育毒性，以及该毒性是否是通过通过干扰有害结局通路而产生。
利用本发明所提供的方法，可以对潜在作用与有害结局路径而进一步产生发育毒性的化学品进行高通量筛查，从而准确快速地判断出化合物是否存在发育毒性，以及该毒性是否是通过通过干扰有害结局通路而产生，因而弥补了现有技术中缺少基于作用模式进行人类发育毒性高通量预测的不足。
在本发明的实施例中，所构建的模型是基于环境污染物干扰有害结局路径信号通路的分子机制构建而成，突破性的将分子水平、细胞水平和个体器官水平上化合物产生的人类发育毒性活性机制产生直接联系，为实现基于计算毒理的人类发育毒性外推提供关键的技术方法和理论基础。
附图说明
图1为本发明的实施例中基于作用模式的人力发育毒性预测的方法的流程图；
图2为基于雄激素受体介导有害结局路径的男性发育毒性预测模型化学品预测流程图；
图3为基于雄激素受体介导有害结局路径的男性发育毒性预测模型建模流程图；
图4为基于雄激素受体介导有害结局路径的男性发育毒性预测模型的第一预测模型的构建流程预测指标图；
图5为基于雄激素受体介导有害结局路径的男性发育毒性预测模型第一预测模型外部验证评估结果图；
图6为基于雄激素受体介导有害结局路径的男性发育毒性预测模型中第一和第二预测模型的预测结果图；
图7为基于雄激素受体介导有害结局路径的男性发育毒性预测模型的预测流程图；
图8基于雄激素受体介导有害结局路径的男性发育毒性预测结果图。
具体实施方式
实施例1
本实施例选取雄激素受体介导有害结局路径的男性发育毒性预测方法为例，对本发明的具体实施过程进行详细说明。图1示出了本发明的实施例中基于作用模式的人力发育毒性预测的方法的流程图，图2示出了基于雄激素受体介导有害结局路径的男性发育毒性预测模型化学品预测流程图，结合图1和图2所示，本实施例中的方法包括：
步骤S100：构建化合物活性数据集，包括基于作用模式，选择有害结局路径，确定第一数量的研究事件，所述第一数量的研究事件能够描述第二数量的毒性终点；收集并整理每个研究事件的相关试验及其或化合物活性数据。
其中，第一数量是指所确定的研究事件的数量，在该实施例中，基于雄激素受体介导有害结局路径的男性发育毒性机制,从分子水平、细胞水平和个体或器官水平收集到了包括分子启动事件(MIE)和关键事件(KE)的七个研究事件的实验数据。七个研究事件分别是配体-受体结合(MIE)、共因子招募(KE1)、DNA结合(KE2)、蛋白转录活性异常(KE3)、转录异常(KE4)、细胞增殖(KE5)和器官发育异常(KE6)，详细地，如图3所示。
除基于动物实验的器官发育异常事件(KE6)外，前六个研究事件的高通量体外(invitro)测试数据来源于由四个美国官方组织开发的ToxCASTTM/Tox21高通量筛查(High-Throughput Screening,HTS)项目，包括美国国家毒理学计划(National ToxicologyProgram,NTP)，美国国家转化科学发展中心(National Center for AdvancingTranslational Sciences，NCATS)，美国食品药品监管局(U.S.Food and DrugAdministration，FDA)和属于U.S.EPA的国家计算毒理学中心(National Center forComputational Toxicology，NCCT)。基于2019年2月26日更新的最新ToxCASTTM/Tox21数据库总共收集到了描述前六个研究事件的20个试验，如表1所示。
表1基于雄激素受体介导有害结局路径的男性发育毒性预测模型所用实验的总结
配体-受体结合(Receptor Binding)作为MIE，是雄激素受体介导有害结局路径活性产生的关键第一步，存在三个分子实验(NVS_NR_cAR,NVS_NR_rAR,NVS_NR_hAR)。
共因子招募(COA Recruitment)作为KE1，存在(OT_AR_ARSRC1_0480，OT_AR_ARSRC1_0960)两个试验，旨在测定雄激素受体介导有害结局路径中招募共激活调节因子(Coactivator,COA)的作用。因此，DNA结合(Chromatin Binding)作为KE2存在一个实验(TOX21_ARE_BLA_agonist_ratio)测定AR是否与DNA存在结合。值得注意的是，基于体外的实验中测定的蛋白生成抑制结果与细胞活性高度相关，即存在细胞毒性的化学品可产生潜在的“假阳性”活性。因此，每个体外试验相关的细胞毒性实验也被选取作为具有“假阳性”活性化合物的阈值过滤器。因此，TOX21_ARE_BLA_agonist_viability作为细胞毒性实验也被选取以防产生假阳性的现象。
蛋白转录活性异常(Transcription Factor Activity)作为KE3存在两个实验(ATG_AR_TRANS_dn,ATG_AR_TRANS_up)。两个体外实验分别测定蛋白转录活性上升趋势(ATG_AR_TRANS_up)或下调趋势(ATG_AR_TRANS_dn)，因此后续模型构建中分别利用进行不同转录趋势的两个模型构建。
转录异常(Gene Expression)作为第KE4存在包括7个试验的激活转录与抑制转录两种测试。化合物激活雄激素受体介导有害结局路径转录包括三个体外试验(OT_AR_ARELUC_AG_1440,TOX21_AR_BLA_Agonist_ratio,
TOX21_AR_LUC_MDAKB2_Agonist)以确定化合物是否存在拟性活性；化合物抑制雄激素受体介导有害结局路径转录包括两个体外试验(TOX21_AR_BLA_Antagonist_ratio,TOX21_AR_LUC_MDAKB2_Antagonist)和相关的两个细胞毒性试验(TOX21_AR_BLA_Antagonist_viability,TOX21_AR_LUC_MDAKB2_Antagonist_viability)以确定化合物是否存在非细胞毒性条件下的抗性活性。
同理，细胞增殖(Cell Proliferation)作为KE5存在包括2个试验(ACEA_AR_agonist_80hr,ACEA_AR_antagonist_80hr)和2个细胞毒性实验(ACEA_AR_agonist_AUC_viability,ACEA_AR_antagonist_AUC_viability)的激活细胞增殖与抑制细胞增殖两种测试。
除此之外，器官发育异常作为KE6，本实施例中选择大鼠Hershberger实验测定器官水平上拟雄活性和抗雄活性对男性相关器官的发育毒性。啮齿动物Hershberger实验被同时确定为U.S.EPA的实验指南(EPA 890.1400)和OECD的实验指南(OECD 441)。在U.S.EPA/OECD指南中，Hershberger实验利用去势雄鼠模型。大鼠在出生后42天左右被阉割，并允许至少7天的术后恢复期，以降低内源性雄激素(睾酮)水平。测定结果是基于5个雄激素依赖的附属性器官(androgen-dependent accessory sex tissues,ASTs)的发育重量变化的Hershberger测定法，包括腹侧前列腺(VP)、精囊(SV)(加上液体和凝固腺)、提肛球海绵体(LABC)肌肉、成对的尿道球腺(COW)和龟头(GP)。当两个以上AST存在明显的器官增重，则化合物存在拟雄效应；反之，当两个以上AST存在明显的器官减重，则化合物存在抗雄效应。拟雄效应和抗雄效应都是化学品产生的不同类型的干扰活性，都会导致男性器官的异常发育，进而产生男性发育毒性(MDT)。最终21个in vitro和in vivo试验被选取用来筛查通过雄激素受体介导有害结局路径作用而产生男性发育毒性的化合物数据集。
将收集的活性数据进行整理，以优化化合物结构与活性数据信息，具体包括：
首先，去除无结构信息(显示“NA”或“FAIL”)、高分子型、离子型、混合型化合物；
然后，利用公式(1)将化合物的活性(activity value,AV)标准化；
公式(1)中，Activity value代表活性强度值，Ki代表抑制常数，Kd代表离解常数，AC50代表半数活性浓度，IC50代表半数抑制浓度，EC50代表半数效应浓度，uM表示微摩尔量。在公式(1)下，针对每个实验，定义活性强度大于等于3的为有活性(AV≥3)，活性强度小于3(AV＜3)的为无活性。
数据整理后，针对每个研究事件，对收集的化合物进行活性分类，包括：
活性(Active)：针对某一研究事件，至少有一个实验存在活性，同时该活性若为抗性体外实验数据，其活性强度需大于同实验下的细胞毒性活性强度；
非活性(Inactive)：针对某一研究事件，所有实验测定结果都为无活性，或仅有抗性体外实验数据，且其抗性活性强度小于等于同实验下的细胞毒性活性强度；
拟性(Agonist)：针对某一研究事件，化合物存在活性前提下，存在拟性活性数据；
抗性(Antagonist)：针对某一研究事件，化合物存在活性前提下，存在抗性活性数据，且抗性活性强度大于细胞毒性活性强度。
最终，收集到的化合物信息如表2所示。
表2基于雄激素受体介导有害结局路径的男性发育毒性预测模型数据收集结果
表2中，“1”表“有活性”，“0”表“无活性”。
从表2中，可以看出，7个研究事件可以描述11个毒性终点，分别为配体-受体结合(receptor binding)、共因子招募(COA recruitment)、DNA结合(chromatin binding)、蛋白转录活性上调(ranscription factor activity-agonist)、蛋白转录活性下降(transcription factor activity-antagonist)、转录激活(gene expression-agonist)、转录抑制(gene expression-antagonist)、细胞增殖-激活(cell proliferation-agonist)、细胞增殖-抑制(cell proliferation-antagonist)、器官发育异常-增重(Hershberger test-agonist)和器官发育异常-减重(Hershberger test-antagonist)，对于每一毒性终点的详细信息见表2。
在本实施例中，所采用的数据可靠有效，从而提高了后续所建模型对人类发育毒性预测结果的可靠性和精确性。
步骤S200：根据化合物活性数据集，结合第三数量的分子描述符库，针对每一毒性终点，通过第四数量的不同种机器学习算法，分别训练出第四数量的对应的QSAR模型，并筛选出每一毒性终点对应的预测效果最佳的QSAR模型，获得的第二数量的QSAR模型的集合形成第一预测模型；
其中，第三数量是指分子描述符库的数量，第四数量是指选用的机器学习算法的种类的数量。
因此针对每个毒性终点，利用表2中的数据，对每个毒性终点进行建模。结合图4所示，在本实施例中，选用了三种分子描述符库和五种机器学习算法。更具体的，这三种分析描述符库包括OEState、Mold2和Dragon v.7，选用的五种机器学习算法包括K近邻算法(KNearest Neighbor，KNN)、朴素贝叶斯算法(Bayes,NB)、随机森林(Random Forest,RF)、支持向量机(Support Vector Machine,SVM)和决策树(Decision Tree,DT)。因此，针对七个研究事件的11个毒性终点，每一个毒性终点训练出5个QSAR模型，总共训练出55个QSAR模型。针对每个毒性终点训练出的5个QSAR模型进行筛选，选出预测效果最佳的QSAR模型，那么对于11个毒性终点，总共筛选出11个QSAR模型，利用这11个QSAR模型的集合构成第一预测模型。
具体的，结合图4，可以按照如下步骤进行。
步骤S201：利用KNIME platform软件中的“Partitioning Mode”模块，将数据集以4:1的比例分为训练集(training set)(80％)和测试集(test set)(20％)，其中，训练集用作模型构建和内部验证，测试集用作外部验证。
步骤S202：利用ChemBioDraw Ultra 14.0软件，检查化合物的结构信息(以SMILES表示)是否正确。如果能够确认化合物的结构信息是正确的情况下，也可以省略该步骤。
步骤S203：选择OEState、Mold2、Dragon v.7三个分子描述符库进行化合物的结构信息数据计算。这三个分子描述符库总共包含6049个1D、2D、3D的分子描述符，利用化合物的SMILES信息在Online Chemical Modeling Environment(OCHEM)在线平台进行计算。主要步骤包括四步结构优化，分别为Standardization、Neutralize、Remove salts、Cleanstructure。对所获分子描述符进行相关描述符的选取，主要步骤包括去除低变化描述符(low variance filter)、去除高度相关性的描述符(high correlation filter)和关键特征描述符的选取(feature importance selection)。
步骤S204：对于每个毒性终点，分别选择五种机器学习算法，训练得到五个的QSAR模型。由于存在11个毒性终点，因此，总共训练得到55个QSAR模型。
在本发明中五种机器学习算法分别为K近邻(K Nearest Neighbor,KNN)、朴素贝叶斯(Bayes,NB)、随机森林(Random Forest,RF)、支持向量机(Support VectorMachine,SVM)、决策树(Decision Tree,SVM)。因此，针对每个毒性终点，训练得到的5个QSAR模型分别为K近邻模型、朴素贝叶斯模型、随机森林模型、支撑向量机模型以及决策树模型。
步骤S205：采用第五数量的不同指标，对QSAR模型的预测效果进行验证，对每个毒性终点，选择具有最佳预测能力的QSAR模型。
其中，第五数量是指指标的数量，在本实施例中，可以采用8种不同指标，对所获得的QSAR模型进行验证，这8种指标分别为：真阳性(True Positive,TP)、假阳性(FalsePositive,FP)、真阴性(True Negative,TN)、假阴性(False Negative,FN)、敏感度(Sensitivity)、特异度(Specificity)、精确度(Accuracy)和曲线下面积(Area UnderCurve,AUC)，在具体的实施过程中，也可以仅采用这8种指标中的一种或多种对最佳预测模型进行筛选。表3给出利用8个指标对所训练的55个QSAR模型的进行验证的结果。
表3基于雄激素受体介导有害结局路径的男性发育毒性预测模型第一预测模型评估结果
并且，在本实施例中，以精确度(Accuracy)为最终准侧，筛选出针对每个毒性终点具有最佳预测效果的QSAR模型，如图5所示。
从图5中可以看出，对于配体-受体结合(receptor binding,Model 1)、共因子招募(COA recruitment，Model 2)、DNA结合(chromatin binding，Model 3)、蛋白转录活性上调(transcription factor activity-agonist，Model 4)、转录激活(gene expression-agonist，Model 6)、转录抑制(gene expression-antagonist，Model 7)、细胞增殖-激活(cell proliferation-agonist，Model 8)和细胞增殖-抑制(cell proliferation-antagonist，Model 9)利用随机森林(RF)算法训练出的QSAR模型的预测效果最好。
而对于蛋白转录活性下降(transcription factor activity-antagonist，Model5)，利用决策树(DT)算法的训练出的QSAR模型预测效果最好。对于器官发育异常-增重(Hershberger test-agonist，Model 10)和器官发育异常-减重(Hershberger test-antagonist，Model 11)，利用支持向量机算法(SVM)训练出的QSAR模型的预测效果最好。
在进行QSAR模型训练的过程中，本实施例中，利用5折交叉验证(five-fold crossvalidation)对构建QSAR模型的训练集进行内部验证，测试数据的稳定性。
额外的，从具有大鼠Hershberger test体内实验结果的48个化合物中，发现其中9个化合物存在基于雄激素受体介导有害结局路径的七个研究事件的所有实验结果。因此，进一步用该9个化合物对第一预测模型中的11个QSAR模型的预测能力进行验证。如表4和图6(A)所示，第一预测模型中的11个QSAR模型能精确预测9个参考化合物的基于雄激素受体介导有害结局通路的七个研究事件的实验结果，其精度高达92％。证明了第一预测模型中的11个QSAR模型能对化合物的基于雄激素受体介导有害结局路径上七个研究事件进行精确而快速的定性预测。
表4基于雄激素受体介导有害结局路径的男性发育毒性预测模型第一预测模型中9个代表化合物的预测验证
表4中，无色表示实测-预测结果一致；灰色表示实测-预测结果不一致；括号外为实验值，括号内为预测值；“1”表征为有活性，“0”表征为没活性。
对于所建立的第一预测模型针对每个研究事件的预测结果还能给出化学品的人类发育毒性机制，为绿色化学品的开发提供有效机制信息。
步骤S300：利用若干种具有体内实验数据的化合物，以及利用所述第一预测模型对具有体内实验数据的化合物进行预测的预测结果，利用朴素贝叶斯算法训练得到第二预测模型；
在本步骤中，利用48个具有体内实验数据的化合物及利用第一预测模型对具有体内实验数据的化合物进行预测的预测结果，预测结果见表5所示。通过朴素贝叶斯算法，将预测结果进行复合叠加，训练得到第二预测模型，在本实施例中，得到的是基于权重的男性发育毒性综合预测模型。利用朴素贝叶斯算法训练得到的第二预测模型实质上是一个权重模型。
表5基于雄激素受体介导有害结局路径的男性发育毒性预测模型利用的48个化合物及其第一预测模型的预测结果信息
需要注意的是，在这一步，由于所运用的体内动物实验数据为大鼠Hershberger试验，该实验不仅能检测化合物的拟雄效应，同时也能检测化合物的抗雄效应。因此，最终在第二预测模型中的存在独立的两个模型：(i)拟雄效应预测模型；(ii)抗雄效应预测模型。当化合物分别经过两个模型预测后，若存在拟雄和/或抗雄活性，则化合物存在男性发育毒性。同时，由于在第二预测模型所能运用的实验数据为48个具有体内实验数据的化合物，为了保证模型的应用域，在该阶段，将利用所有的数据作为训练集进行QSAR模型构建，并进行内部验证。同样的，利用步骤S200中八个指标中的一个或多个对模型的预测效果进行验证。表6图6(B-C)显示了第二预测模型中的拟雄活性和抗雄活性预测模型结果。
表6基于雄激素受体介导有害结局路径的男性发育毒性预测模型第二预测模型评估结果
结果表明，对于抗雄活性预测模型来说，可以完全预测出48个化合物是否具有抗雄活性，其预测精度可达100％；而对于拟雄活性预测模型，尽管五个拟雄化合物能被完全预测出，即不存在假阴性结果，但是36个无拟雄活性的化合物被预测为假阳性，即被错误预测成具有拟雄活性的化合物，其精度只有25％。尽管在监管角度上，0％的假阴性率保证了拟雄模型的有效性，但是过度的错误预测也达不到预测目的。
因此，本发明的实施例中，详细研究了48个化合物的化学结构及其特征，发现5个拟雄化合物都为甾体类(steroid)化合物(testosterone propionate,17-methyltestosterone,trenbolone,methyl-1-testosterone,testosterone)，而其他43个无拟雄活性化合物都不是甾体类化合物。因此，新增一条筛查条件“被预测为有拟雄活性的化合物是否为甾体化合物？”，过程如图7。
建模过程中，该筛查条件利用化学结构相似性(chemical similarity)进行预测和筛查。详细地，将五个拟雄化合物作为模板化合物(positive controls),预测化合物的结构与模板化合物的结构进行一一匹配。化学结构相似性利用Tanimoto SimilarityScore进行表征。而相似性打分是利用PaDEL-descriptor software包含的12个分子结构指纹库进行打分。12个分子结构指纹库分别为Fingerprinter,Extended Fingerprinter,Estate Fingerprinter,GraphOnly Fingerprinter,MACCS Fingerprinter,PubchemFingerprinter,Substructure Fingerprinter,Substructure FingerprintCount,KlekotaRoth Fingerprinter,KlekotaRoth FingerprintCount,AtomPairs2DFingerprinter,和AtomPairs2DFingerprintCount。Tanimoto Score的输出值均在0-1之间，score越大，化学相似性越高。因此，在本实施例中，设定Tanimoto Score的cutoff值为0.8，当被测化合物与5个拟雄化合物中的至少一个化合物相似度≥0.8，则证明，被测化合物为满足拟雄活性的甾体类化合物。新增筛查条件后的拟雄预测模型预测能力得到极大提升，如表6和图6(D)，可以完全预测出48个化合物是否具有拟雄活性，其预测精度可达100％。
步骤S400：将待测化合物输入到所述第一预测模型中进行定性预测，并将定性预测结果输入所述第二预测模型，完成对化学品的人类发育毒性进行预测。可以看出，待预测化合物实质上是要通过两层预测，第一层是通过第一预测模型进行定性预测，第二层是通过第二预测模型进行综合预测。
在本实施例中，按照图2和6所示流程，采用化合物flutamide对本实施例中所提出的方法进行测试验证，化合物首先输入到第一预测模型中进行定性预测，第一预测模型中的11个QSAR模型预测待测化合物在基于雄激素受体介导有害结局路径上关于七个研究事件的活性/非活性；然后，将第一预测模型的预测结果输入到第二预测模型中，在本实施例中，因为第二预测模型中包括拟雄活性和抗雄活性两个预测模型，因此，通过第二预测模型分别预测化合物的拟雄活性和抗雄活性。最终，当化合物存在拟雄活性和/或抗雄活性时，该化合物存在男性发育毒性。图8给出了本实施例的综合显示结果。分别显示了(i)待测化合物信息，包括化合物名称、CASN、结构信息、男性发育毒性预测结果；(ii)详细的男性发育毒性预测结果(表格格式)；(iii)详细的男性发育毒性预测结果(镭射图格式)。
进一步地，实施例中上述过程，仅能预测小分子有机物，无法预测重金属类、混合物、离子态结构化合物的男性发育毒性，因为存在致毒机制差异性问题。因此在本实施例中，利用训练集中化合物的26种物理化学性质来描述相关QSAR预测模型的应用域(Application Domain,A D)。26种物理化学性质包括1D、2D、3D，分别为nAcid，ALogP，AMR，apol，naAromAtom，nAromBond，nAtom，nHeavyAtom，nBonds，nBondsD，nBondsT，nBondsQ，bpol，ETA_Alpha，FMF，nHBAcc，nHBDon，TopoPSA，VABC，MW，AMW，XLogP，TPSA，nRing，nRotB和nRotBt。每个分子描述符所代表的详细物理化学性质列于表7。化合物的物理化学性质用PaDEL-Descriptor软件的1D&2D&3D分子描述符库进行计算。
表7与雄激素受体介导有害结局路径的男性发育毒性预测模型应用域选择的26种物理化学性质信息
应用域(AD)的构建过程包括四个步骤。首先，将用于建模的训练集所有化合物计算出各自的26种物理化学性质，并利用KNIME platform软件的normalizer模块对每种性质数值进行Min-Max Normalization标准化，其数值在0-1之间。其次，利用RComplexheatmappackage中的欧式距离(Euclidean distance)计算出该模型中训练集化合物的距离矩阵(Distance Matrix,DM)，DM即为QSAR模型的AD。再次，计算出被测化合物的26种物理化学性质数值，并基于相应QSAR模型的AD(即训练集化合物的26种物理化学参数)进行数值标准化，将标准化的数值与模型的DM进行相似性打分，并找出与被测化合物最相近的化合物和相似性距离(similarity distance)。最后，当被测化合物与最相近化合物的相似性距离≥0.5，则判定被测化合物在AD内，反之，不在AD内。
对于第一预测模型中的每个QSAR模型都存在一个AD。因此，整体上，一种基于雄激素受体介导有害结局路径的男性发育毒性预测方法的应用域是11个QSAR预测模型的并集。
从图8可以看出，被测化合物flutamide首先进行tier 1的11个QSAR模型预测，其预测结果与是否在AD内(括号内注明，“in”表在AD内，“out”表不在AD内)都已标注出。可发现对于化合物flutamide而言，它都在11个QSAR模型的AD中，证明了其预测的有效性。
目前基于传统动物实验(大鼠Hershberger实验)的有效实验数据极少。2018年，OECD的研究者进行了极其系统的文献搜索进行大鼠Hershberger实验数据的整理，以获得有效的动物实验数据库(https://www.regulations.gov/document/EPA-HQ-OPPT-2009-0576-0008)。结果发现，在符合OECD/US EPA 规范的前提下，接近3200条研究数据，只有134个化合物初步满足OECD/US EPA的规范。进一步发现，只有48个化合物至少存在两个及其以上的研究结果，且研究结论一致；而有24个化合物至少存在两个及其以上的研究结果，单结论不统一，并无法确定其在个体器官水平上的发育毒性。基于48个具有有效动物实验结果的化合物而构建起的QSAR模型，虽然在AD内能有效对化合物进行毒性预测(Accuracy＝1,图4)，但是AD范围的过于狭小问题极大局限了其预测模型的实际应用。本实施例基于雄激素受体介导有害结局路径的男性发育毒性预测方法是基于雄激素受体介导有害结局路径，利用多维QSAR模型构建起的二层预测模型。该模型不仅仅可以解析化合物的干扰机制(干扰哪个关键事件)，同时多维模型结合的方式也极大的扩展了仅基于动物实验数据的狭小应用域，大大增强了模型在实际应用中的预测范围与预测能力。
上述实施例仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和等同替换，这些对本发明权利要求进行改进和等同替换后的技术方案，均落入本发明的保护范围。
