2.3 危害特征描述原理及方法

危害特征描述是加工过程风险评估最关键步骤之一,评价方法与评价的危害对象一致,由于该步骤实际上主要是对基准-剂量关系的探讨,因此也叫做剂量-反应(效应)关系评估。同时很大程度上,研究结果对质效应(如酶活力抑制率等)较量效应(如活细胞减少数量等)而言更有意义,危害特征描述在某种程度上也简称为剂量-反应评估。

目前,除传统用于毒性评价的无明显损害作用水平(NOAEL)方法外,新开发出的方法包括基准剂量法(BMD)、化学物系数调整法、概率剂量反应评价法等,其中基准剂量法因其对数据量和使用对象专业数学水平要求相对较低,且大多情况下方法精准度较传统NOAEL而言极大提升(顾刘金等,2006;何贤松等,2013),因此该方法被接受认可,应用推广程度最高。

2.3.1 有阈值、无阈值与剂量-反应评估

2.3.1.1 有阈值与剂量反应评估

有阈值(threshold)(Schwartz et al,1995)指化学物质引起受试对象出现可指示最轻微的异常或改变时需要的最低剂量,也就是一种物质使机体(人或实验动物)刚开始发生反应的剂量或浓度。阈值又分为急性阈值及慢性阈值。急性阈值为一次暴露所产生不良反应的剂量点,而慢性阈值是长期不断反复暴露所产生不良反应对应的剂量点。

一种化学物对每种反应都可有一个阈值,因此一种化学物可有多个阈值。对某种反应,对不同的个体可有不同的阈值,同一个体对某种反应的阈值也可随时间而改变。就目前科学发展现状,对于某些化学物和某些毒效应还不能证实存在阈剂量(如遗传毒性致癌物和性细胞致突变物)。阈值并非实验中所能确定,在进行风险评估时通常用NOAEL或最低可见有害作用水平(LOAEL)作为阈值的近似值。在利用NOAEL或LOAEL时,应说明测定的是什么效应、何种群体和其染毒途径。当所关心的效应被认为是有害效应时,就确定为NOAEL或LOAEL(Calabrese and Baldwin,2003)。

在风险评估为风险管理提供科学依据时,最常用到的是制定最大残留限量(MRLs)标准,但在标准制定过程中,最关心的不是一条曲线,而是一个剂量点,而该点是判定有害与无害之间的一个临界点,即NOAEL或LOAEL。该点可通过剂量-反应模型获取或直接通过动物试验获取,然后通过该点来推导MRLs(Chou et al,1998)。

另外,阈值表征危害对个体产生的不良反应,但实际上该不良反应常常取决于观察指标、检测技术与检测方法、仪器灵敏度、试验设计等多重因素,于是很难严格规定该“安全”与“不安全”剂量点,故阈值实际上更具备理论研究价值和意义。阈值实际上具有两层不同的含义。一是基于科学含义,指在不良效应发生情况下的暴露水平,如有生理刺激,但未形成不良反应;二是阈值代表一个水平,在该水平上没有不良反应,但这仅仅指在该水平下的反应极为不明显,以至不能被观察和监测到,如NOAEL(Leisenring and Ryan,1992)。在类似情况下,相对于实际结论而言,这更基于分析者或采纳分析方法的原则所界定的阀值。在第一层含义中阈值也许需要考虑剂量-反应模型,阈值剂量的引入将剂量-反应模型分为不同情况:一方面在阈值以下,有效剂量为零;另一方面,在阈值以上,有效剂量是剂量-反应模型上剂量点减去阈值剂量。阈值通常不能提高模型的拟合效果,而阈值的置信度一般非常大。有阈值与无阈值是从毒理学角度划分风险评估范畴的一个重要分界点(Leisenring and Ryan,1992)。另外,对于有阈值化学物风险评估,不同机构赋予其不同的剂量点表征名,世界卫生组织(WHO)采用了每日允许摄入量(ADI),世界卫生组织/国际化学品安全机构(WHO/IPCS)采用了TI,美国环保署采用了参考摄入量(RfD)或参考浓度(RfC),美国毒物与疾病登记署(Agency for Toxic Substances and Disease Registry,ATSDR)采用了最小风险水平(minimum risk level,MRL),加拿大卫生部采用了日允许摄入量/浓度(tolerable daily intake/concentration,TDI/TDC)。

2.3.1.2 无阈值与剂量-反应评估

那些具有遗传毒性的致癌物,因其无典型的所谓引起受试对象出现可指示最轻微的异常或改变时需要的最低剂量,一旦产生暴露,便可指示目标不良效应,这类效应即无阈值(unthreshold),具有这类效应的物质即无阈值物(Christenson et al,2014)。一般该类物质安全水平通过数学模型来推导,虽然有几种数学模型,但通常采用线性多级模型(超线性多阶打击模型等)。这些模型在某个特定暴露水平具有可置信区间的上限值与下限值,在下限值中可包括零值,而所谓安全水平一般以保守方式表示为终生癌症超额危险的上限估计值10-6(即摄入含该无阈值物质的食品70年,每10万人中会产生一个癌症病例)。该数值不一定等于接触该浓度水平所引起的癌症病例,这是考虑大量不确定因素后的最大潜在危险水平,很可能实际风险水平比该水平还要低,也可能还要高,但低剂量暴露水平的风险程度往往无法由实验验证。

2.3.2 主流的剂量评估方法

2.3.2.1 基准剂量法

基准剂量(benchmark dose,BMD)方法是美国人Crump于1984年首先提出并得到美国环境保护署(U.S Environmental Protection Agency,EPA)和美国食品与药品管理局(U.S Food and Drug Administration,FDA)的极大重视,随后加以推广和应用,目前BMD法是当前国际组织和发达国家继Lehman和Fitzhugh于1954年提出商值法(即NOAEL法)后最积极推广并用于剂量-反应评估的重要技术程序之一。该方法有诸多优势,主要包括:一是BMD利用了所有的动物实验数据(二分类型数据),选用合适的剂量-反应模型,并通过统计方法(极大似然估计等)得到结果,因而对实验设计依赖性减小,消除了实验设计的随意性,提高结果可靠性、准确性;二是 BMD可根据基准剂量反应(BMR)水平的选择,了解BMD不同反应水平;三是需获取反应剂量95%可信区间下限值,必须纳入组数、每组受试对象数及终点指标观察值的离散度。如果资料质量水平不高,则可信区间会宽(即不确定性提高),BMD水平也相应降低,使得保护水平提高。

BMD是对于一个特定的观察终点,应用其所有剂量-反应数据来估计总体剂量-反应关系的形状。BMD是一个剂量水平,从估计的剂量-反应曲线上获得,与反应的特定改变有关,该反应被称为基准剂量响应(benchmark response,BMR)(Gordon and Malley,2014)。基准剂量下限(BMDL)是BMD的可信区间下限,该值通常被用作参考点。图2-1说明了BMD方法的基本概念。从图2-1可以看出,计算出来的BMDL(如BMR的5%)可解释为:基准剂量下限值5(BMDL5)=反应可能低于5%的剂量,此处的“可能”由统计学上的可信限来确定,通常是95%的可信限。以观察到的反应的平均值(三角形)及其可信区间进行作图。实体曲线是拟合的剂量-反应模型。根据该曲线可确定BMD,BMD通常是指与较低的且可测量的效应(称为BMR)相对应的剂量。虚线分别代表效应95%可信限的上限和下限,是剂量的函数。它们与水平线的交叉点所对应的横坐标分别是BMD的下限和上限,即基准剂量下限(BMDL)和基准剂量上限(BMDU)。

图2-1 基准剂量(BMD)方法的基本概念(Gordon and Malley,2014)

BMD指基准剂量;BMDL指基准剂量下限;BMDU指基准剂量上限;BMR-5%指基准剂量反应的5%

在利用基准剂量(BMD)方法推到每日允许摄入量(ADI)过程中会用到不确定系数,种间差异和人类变异性是能常用于导出每日允许摄入量/日允许摄入量(ADI/TDI)的100倍不确定性系数的依据,与种间差异和人类变异性有关的不确定性同样适用于基于动物数据的某项有效性量度(MOE),但是存在附加的不确定性,所述附加的不确定性与试验的/可观测的范围以下的剂量-反应关系的性质有关,与对于生成变异细胞以及后来克隆扩充并发展成为某种癌症具有决定性的过程中的遗传多态性的影响有关。因此,一项100倍的MOE不足以反映这样的事实:起始点基准剂量下限(BMDL)值不能认定为一项阈值或与作用方式有关的附加的不确定性。

ADI=BMDL/UF  (2-2)

式中,ADI为每日允许摄入量;BMDL为基准剂量下限值;UF为不确定系数。

线性低剂量外推的应用,使用5%反应的BMDL来评估一项百万分之一的一生风险,相当于一项50000的有效性量度(MOE)。该方法必须通过一定不确定系数(uncertainty factor,UF)加以修正后才能获取每日允许摄入量(ADI),有的研究不能直接获得NOAEL,所以在获得真正意义上的NOAEL时,必须考虑到UF。目前,UF有6个,这些UF合理被应用于公式(2-2)中,用于合理估计ADI。6个UF如下。

UF1——种内差异,主要为人类之间差异。这个时候数据直接来自人体试验,合理结果采取外推UF设定为10。

UF2——种间差异,主要是实验动物外推到人类之间差异。对于BMD而言,当人群暴露研究不能直接采纳时,从实验动物外推到人类时,采用UF为10;对呼吸暴露,人相对浓度NOAEL作为获取ADI方法时,该UF为3,因为人相对浓度时已考虑了毒代动力学的差异。

UF3——长期、中期及短期暴露差异。当从慢性(成本及可操作性等因素决定在动物实验中不会采纳)、亚慢性及急性动物实验推导时,采用UF为10。

UF4——LOAEL推导不可见作用最高剂量(NOAEL)带来的差异。采用UF为10。

UF5——当资料不完善或可信度不是相当高时,如数据是否来自GLP实验室,采用UF为10。

UF6——修饰系数(modifiying factor,MF),该系数在0~10之间,数据可信度差异(不同对比试验分析),数据是否经过不同实验室重复试验,数据是否全面等。

经过这些不确定系数一系列复杂计算修正后的NOAEL才能采用。UF一般基数为10。目前一般选定的UF为前两项,即UF1与UF2,分别为种间差异和种内差异,每项取10,所以不确定系数为100。

2.3.2.2 化学物系数调整法

1994年国际化学品安全机构(IPCS)召开会议后讨论了化学物质的调整系数概念(chemical-specific adjustment factors,CSAFs),该方法目的是为更真实反映人体对于不同物质的毒代动力学(toxicokinetics,主要表示化学物在人体内经过复杂代谢后的产物)及毒效动力学(toxicodynamics,主要表示化学物在人体内经过复杂代谢后对人体的作用)后对UF的贡献率[该方法被世界卫生组织/国际化学品安全机构(WHO/IPCS)所应用],取代原来非常传统的一律采用UF为10的做法,这样,减少了对经验数学模型的依赖(Organization,1994)。如图2-2所示。

图2-2 不确定系数(UF)引入毒代动力学与毒效动力学(Organization 1994)

H—不同人之间的差异;A—人类与实验动物之间的差异;D—毒代动力学;K—毒效动力学

2.3.2.3 结构活性方法

(1)二维定量结构-性质/活性关系(QSPR/QSAR)基础理论

化学物在生物体的吸收、分布、代谢和排泄(ADMEs)取决其理化性质。因此,只有理解了外源性化学物的理化性质,才可正确评估其生物活性、生物利用度、在生物体内转运过程及在生物系统不同隔室之间的分布情况。二维定量结构-活性关系(2D-QSAR)揭示化学物疏水性、电性和立体性等理化参数与生物活性之间的相关性。其模型框架如图2-3所示。由于所有参数与系统自由能相关,故为线性自由能相关模型(Carlsen,2009)。

图2-3 二维定量结构-活性关系(2D-QSAR)模型框架(Carlsen,2009)

LC50—半数致死量浓度;LD50—半数致死量;EC50—半数有效浓度;IC50—半数抑制率浓度

(2)数据定义

①化学物结构参数 化学物结构参数可称为自变量,该变量可通过实测也可通过计算获取,可以是物化参数,也可以是量化参数或图形参数。化学物结构参数根据模型框架,包括电性参数、输水性参数、立体参数等。

电性参数包括哈米特常数(Hammett)、原子电荷、最高占有分子轨道、最低空轨道、前线轨道电荷密度以及解离常数等。其中,Hammett用得最多,以σ表示,,其中kxkh均表示解离常数,但前者是取代基x存在时酸解离度(25℃,水),后者是苯甲酸解离度(25℃,水),该常数根据取代基对苯甲酸解离度影响提出,在芳香族取代基中应用最广泛。同时因为吸电子基使苯甲酸解离度增加,因此吸电子基为正值,反之,斥电子基为负值。常可利用间位或对位取代基常数,由于邻位取代基既有电性效应,又有熵效应,所以一般而言,不利用邻位取代基常数。

输水性参数包括分配系数和疏水性常数。许多方法可用于分配系数(lgKow,也称为lgP)的实验测量,lgP的预测可使用贡献法或定量结构-性质关系(QSPR)模型。针对QSAR关系建模和环境化学物规趋建模中所涉及的其他理化性质,也有许多实验测定和预测方法。

立体参数包括塔夫脱(Taft)立体参数、范德华立体参数以及最小立体差异参数等。其中Taft立体参数应用最多,由于水解反应速率受到取代基Y的诱导,以及共振及立体效应的影响,所以在酸性条件下,酯的水解几乎完全取决于立体因素。

②生物活性数据 生物活性数据可称为因变量,由试验测定,可以是连续型也可以是二分型数据。同样可以基于体内实验获取,也可通过体外实验测定获取。因变量通常以产生标准生物效应时化学物的剂量或浓度负对数(-lgC)表示,而标准生物效应则采用剂量效应曲线的敏感部位,如酶抑制剂抑制50%时的浓度水平(IC50)、半数致死剂量(LD50)、最小抑菌浓度(MIC)以及产生50%最大效应的化学物浓度水平(EC50)表示。另外,数据最好选取化学物作用靶器官的生物活性数据进行建模。半数致死剂量(LD50)最好来自啮齿类动物急性毒性实验数据,因为此类数据具有更好的统计学和毒理学双重意义。影响LD50的因素包括生物物种、品系、性别、生理状况及在实验室分析测试条件等。如采用Sherman品系成年大鼠,在严格且相同的实验条件下进行经口农药毒性试验,结果表明大多数农药对雌性大鼠的毒性强于对雄性大鼠的毒性。每种农药对雌、雄两种大鼠的LD50的置信限不发生重叠,这表明雄性和雌性大鼠之间存在着敏感性差异。

③置信限 由于二维定量结构-活性关系(QSAR)模型开发往往数据来源丰富,因此首先必须对数据有效性进行确认和验证。通常,生物数据需要经过转换才能用于QSAR建模。如化学物的LD50需转换为物质的量单位后才能用于结构活性建模。此外,使用经典数学统计方法(如回归分析)时,为避免统计误差,生物数据须转换为对数形式,然后才可用于建模。按照惯例,因为活性较高的化学物实测值较大,所以最好将这些数值转换为负对数。生物数据集经过对数转换后,最好能跨越几个数量级,这样利于QSAR建模。最后,如果LD50来自多个数据来源或来源不明,则需将数据分类。有些毒理学数据分类(如致癌性)基本遵从布尔表达方式(即致癌/非致癌)。

④指示符变量 指示符变量,亦称为一维描述符,可影响化学物生物活性的分子结构特征(如原子或官能团)。指示符变量是虚拟参数,当不能用结构参数或理化参数描述时,指示符变量往往可以解释该结构对呈现的活性的贡献率,对复杂化学物进行QSAR建模,指示符变量意义在于初筛结构-活性关系的相关性并在解释作用机理方面起到重要作用。

指示符变量只能用0或1来描述,表示编码化学物分子中一个结构元件。指示符变量也被称为虚拟变量或布尔描述符,这些描述符是用于描述分子最简单的方法。布尔描述符尤其适合编码化学物分子中取代基的位置或数目。Free-Wilson方法即是基于此类描述符的方法,该方法可定量分析取代基对化学物生物活性的贡献率。根据假定,分子母环上的取代基对于化学物生物活性的恒定贡献呈现简单加和性原则。方法基本步骤是生成一个由0和1所组成的数据矩阵,该矩阵中每一列对应于分子上特定位置的一个特定取代基,该取代基可视为一个独立变量(即自变量)。数据表中还有一列因变量数据(即生物数据)。因变量和自变量之间应用多元回归分析,使用经典的统计方法来实现拟合并检验拟合优度。模型的回归系数代表了取代基对化学物生物活性的贡献。目前为止,Free-Wilson方法已广泛应用于QSAR建模。

对于分子结构的描述,除简单的布尔描述外,分子上特定的原子或官能团出现的频率也可作为分子描述符。然而,出现频率必须足够高才能确保在回归分析中得到良好的回归系数。即使不使用回归分析或其他统计方法,此类描述符也可用于QSAR的预测。事实上,指示符变量可独立预测一些特定终点。在可预测特定终点的情况下,指示符变量被称为结构警报。

(3)分子描述符选择及确定

有机分子的整个结构还可被描绘为不带氢原子结构,由此推导出的数值描述符称为拓扑指数。目前有许多算法可计算拓扑指数,对于所有已有和未合成的化学物均可很容易使用这些描述符来计算。同时,将几种描述符有机结合起来可综合多变量描述化学物分子。

迄今为止,可计算出二维分子的描述符数以千计,但只有具备理化性质意义且彼此互不相关的描述符才可用于QSAR建模。在实际应用中,计算出的部分描述符,由于过度统计拟合,以至即使该描述符具有很好的预测能力,从机理角度而言,这些描述符依然毫无意义。一个集合里选取的分子描述符应尽可能彼此独立(即正交关系),如建模所用的描述符彼此相关性太强,则由于偶然相关性,很可能得到的是自圆其说,不能现实应用的模型。为了避免此问题,根据不同描述符性质,可选择使用主成分分析法或对应因子分析法。上述两种方法为线性多变量分析,即通过原始变量的线性组合而创建新的变量,创建的新变量分别称为主成分和因子,这些新变量呈现彼此正交关系,可用于描述符的数据矩阵,以显示和解释描述符、化学物和活性之间的关系。

(4)二维定量结构-性质关系(QSPR)建模方法

QSPR建模需使用不同数理统计方法,从现有数千个理化性质参数中筛选出“最好”的几个参数,然后基于这几个参数,生成QSPR/QSAR模型。典型的统计方法包括简单线性回归(simple linear regression,SLR)、逐步线性回归(multiple linear regression)、多级逐步线性回归(stepwise multiple linear regression)、遗传算法(genetic algorithms)、神经网络算法、聚类方法(cluster analysis,CA)、因子分析(FA)、主成分分析法(principal compinent analysis,PCA)以及偏最小二乘法(partial least squares,PLS)等。这些方法可根据结构、生物性质、取代基及理化性质对化学物进行分类,便于针对不同分类的化学物进行针对性研究。

①回归分析 由于变量之间关系复杂,无法得到精确数学表达式,或由于误差导致结构参数与生物活性参数之间关系存在某些不确定性,需要采用统计方法,在大量实验和观测结果中寻求隐藏在随机性后面的统计规律,这样的方法称为回归关系。对这类回归关系进行分析的方法称为回归分析。回归分析是对一组数据进行最小二乘拟合处理并建立参数之间的相关性过程,当有几种结构性质可能对生物活性有影响,产生贡献率时,可用多元回归来进行预测和分析,拟合函数的统计评价也是该分析方法的一部分,验证的方法包括Free-Wilson方法及Hansch方法等。

在该方法中,判定预测情况优劣的指标有样本数量n、复相关系数R、标准偏差S以及显著性检验F

样本数量n:在QSPR/QSAR建模过程中,要求样本数量不能过少,最好越多越好,且最好是参数数量的5倍以上,具体数目一般不能低于10个,10个以下的样本数量不具有太大意义。最后建立回归方程,在建立方程时可删除一些异常样本,但应遵循删除样本不能超过总样本数量的10%,同时对异常样本应一个一个删除,然后检验,慎重处理,如此反复,直到获取最好的回归模型。

复相关系数R:复相关系数R可描述方程中生物活性参数与结构参数之间线性关系的大小程度,复相关系数一般由公式(2-3)计算得出。

   (2-3)   

式中,n为样本数量;yexp,i为第i个样本的实验值;ycalc,i为第i个样本的计算值;为所有样本的均值。

复相关系数R的取值范围为0≤R≤1,且R越接近1,回归效果越好。

标准偏差S:标准偏差S是判断回归效果优劣的标准,反映数据离散程度,其计算方程见公式(2-4)。

   (2-4)   

公式(2-4)中与公式(2-3)中相同符号所表示的相关参数及其含义一致;m为回归方程中化学物结构参数的数量。

显著性检验F:显著性检验主要用于判定样本是否具有显著性差异,F越大,则差异性越大。在最后的方程中,化学物理化性质参数间呈现正交性,判断正交性的简单办法是将选定的化学物取代基进行两两对应的回归分析,相关系数应低于0.8。

回归包括一元线性回归和多元线性回归,一元线性回归[见公式(2-5)]表达生物活性参数(Y)和化学物结构参数(X)之间存在的关系。但实际上往往影响生物活性的参数是多个,即X=(X1X2X3,…),研究YX=(X1X2X3,…)之间定量关系的问题就是多元回归。多元线性回归可优化和预测同源先导化学物活性,分析化学物作用机理,推测受体模型结构,并获得因果模型且物理意义明确。

Y=aX+b  (2-5)

式中,ab为回归方程的回归系数。

判断生物活性与理化性质参数之间是否存在线性关系,可通过对回归系数进行t检验来判定,如果发现有的偏回归系数不显著,则要从回归方程中删除没有显著作用的理化性质参数,不可同时将几个不显著的参数同时删除,应当先删除t值最小的,重新计算回归方程后对新的方程再行检验,一一删除,直到回归方程中所有生物参数均显著为止。

均方根误差(RMSE)或均方误差(MSE):

   (2-6)   

式中,yk表示试验值;表示预测值;nk表示样本数量。

②聚类分析 聚类分析又称为类聚群分析、群分析等,是按照不同化学物、取代基团或不同结构信息参数之间的相似程度,用数学方法将这些信息进行多元统计的方法,该方法可对基于某化学性质参数为基准的不同化学物分类,使相似化学物、相似取代基等物质汇总一起,节省建模时间和成本。

③主成分分析法 主成分分析法在QSAR建模中,由于所选用的结构参数、理化参数之间往往存在不同程度相关性,使得提供的信息发生重叠,从而掩盖本质信息,试图透过重叠信息探明本质,可通过线性组的方法来实现。当x个化学物中用y个参数描述,而其中z个参数之间可能存在相关性,则可用主成分分析方法给出z个彼此无相关性的新综合参数,即原始参数的线性组——主成分。合理从z个主成分中挑选出能代表关键信息的主成分,即可获取由原始参数包含的绝大部分信息,根据主成分值绘制主成分图,以考察化学物更为合理的分类情况。主成分分析随着次序增加,其信息重要性依次降低。

④偏最小二乘法(PLS) 偏最小二乘法考虑化学物结构参数相关的同时,还考虑生物活性参数的情况,以此探明上百甚至上千个独立结构参数与生物活性参数之间的相关性。该方法首先用交叉验证法求取预测残差,以选出最大预测模型,然后用非交叉方法计算回归模型表达式,回归系数R2越大,标准差S越小,表示相关性越好,QSAR模型的预测能力越佳。目前,PLS是当前应用于QSAR建模最好的方法,其主要有如下优点:首先,多元线性回归方法要求提出相关的结构参数,过分强调回归方程的相关系数,而对于相关结构参数的信息提取有所忽视,而PLS方法对该参数之间非共线性要求并不苛刻,较合理地提取参数信息,使得PLS方程更大、更全面且具有更好的统计意义;其次,传统主成分分析方法中注册成分不包含生物活性参数的成分,而PLS 法主成分中包含该参数;另外,当结构参数超过样本总量时,PLS法能得出具有统计意义的方程,而传统线性回归方法则要求样本数量至少是结构参数的五倍以上才能获取;最后是PLS方法运用了交叉验证,以降低偶然相关性,比多元逐步回归偶然相关的可能性小,也使得模型过度拟合大为改进,因此预测能力及准确度加强。

基于一套训练数据集建立起来运行良好的QSPR模型未必能很好地预测训练集以外化学物的理化性质。因此,评估QSPR对未知化学物的预测能力非常重要。一般方法是运用建立该QSPR模型训练集以外的几种化合物的理化性质的数据集(该数据集应当非常精准)进行验证,纳入验证的数据集被称作测试集。测试集相关化学物理化性质类似于训练集相关化学物的理化性质。因此,通常做法是将化学物的所有数据集分成两组,较大的一组(50%~95%)作为训练集,较小的一组(5%~50%)作为测试集。如果测试集标准误差远大于训练集标准误差,则该QSPR模型预测能力不佳,反之则相反。如果化学物总数据集本身较小,不足以拆分成为训练集和测试集。在此情况下,可使用内部交叉验证,即从训练集中随机抽取一种化学物数据集进行验证,同样随机抽取另一种化学物的数据集用于开发QSPR模型,然后反过来,原来用于开发模型的数据集用于验证,而用于验证的数据集用于建模,直到所有化学物验证数据集和建模数据集组合全部通过上述方式尝试完毕。然后计算交叉验证的R2值,称为Q2,这是QSPR建模内部评价模型预测能力的一个指标,当Q2≤0.5时,结果可接受。

(5)模型筛选和优化

有许多方法可将生物数据与分子描述符之间建立相关性模型。如何选择最适合的统计方法,以及应当注意避免的主要误区至关重要。建立QSAR模型,首先最简单的方法是以图形表示建模过程中用到的生物数据及各种分子描述符。简单的散点图可帮助建模者辨识和舍弃没有意义的分子描述符,并有助于避免模型设计的偏置。

在结构-活性关系(SAR)或QSAR建模过程中,用于计算结构–活性模型的最佳统计方法选择非常重要,但往往被忽略。很多人在描述化学物分子活性与其分子描述符之间关系时常常假定为线性关系,而对于连续数据建模,通常使用线性回归分析和偏最小二乘分析,对于分类数据建模,通常使用线性判别分析。事实上,上述假定常常会出现错误,这是因为大量SAR模型通常呈现非线性关系,因此只有使用纯粹的非线性统计方法,才能正确地表征这种非线性结构-活性关系。然而,这些非线性统计工具的正确使用需要一定统计知识,只有将这些参数调整为正确值,才可计算出可接受的预测结果。因此,最好的策略是先用线性方法(回归分析或PLS分析),再用非线性方法(如三层感知器方法,见图2-4)检验预测结果质量是否得到提升(Katritzky and Gordeeva,1993)。

图2-4 三层感知器示意图(Katritzky and Gordeeva,1993)

(6)模型诊断

建模的重要步骤首先是验证模型的统计意义并检查模型是否能够正确地预测化学物结构-活性。如果SAR模型旨在预测二分类型数据结果(如阳性或阴性、有或无等),那么模型预测结果会有4种可能性,即真阳性(TP)、真阴性(TN)、假阳性(FP)和假阴性(FN)。假阳性结果是把真正的“无活性”(或称为“阴性”)错误地归类为“活性”(或称为“阳性”)。从这四种类型结果可计算出各种参数,最重要的有以下几种:

①灵敏度或真阳性率(TPR)=TP/(TP+FN)。

②假阳性率(FPR)=FP/(FP+TN)。

③特异性或真阴性率(TNR)=TN/(FP+TN)。

④准确性=TP+TN/(TP+FP+TN+FN)。

上述TPR与FPR坐标图被称为标准特征曲线(ROC)。对于不同分类模型的比较,ROC曲线非常重要。

如QSAR模型计算基于连续型数据[如lg(1/LD50)],则建模验证方法是用每一种化学物活性预测值减去该化学物相应活性实验值,得到残差,过大残差被称为离群值。离群值的存在表明建模者还需更多地了解QSAR模型(Escher and Hermens,2002)。此外,导致离群值的原因如下。

①用于QSAR建模计算的一组化学物选择不当。

②生物数据试验值不准确或错误。这可能由于建模终点和实验条件等方面的差异导致,还可能由于数据整理和优化过程中出现错误等。

③模型所选择的描述符不正确。

④化学物在其生化作用位点的分子机制不同于该数据及其他化学物。

⑤化学物生成了一个或多个代谢产物,这些产物在其生化作用位点的分子机制不同于该数据集其他化学物。

⑥选择的统计方法不恰当,不足以将生物活性和选定分子描述符之间确立函数关系。

只有找到正确的原因才能删除处于离群的化学物,方法是添加更多的化学物,添加其他分子描述符或选用另一种统计方法。改进QSAR模型非常耗时,尤其该模型来自非线性方法(如人工神经网络)。

(7)模型预测性能评估

模型建成后,都需要评估其预测性能。方法的第一步是内部验证,通常采用留一法交叉验证(LOO)或留n法(leave-n-out)交叉验证,后者通常优于前者。留一法顾名思义是从建模数据集删除一个化学物,其余化学物用于拟合模型,删除的化合物用于验证该模型的预测能力,然后周而复始,直到数据集的每个化学物都被轮番删除一次。使用该过程可生成各种统计数据,如预测残差的平方、交叉验证的R2(被称为Q2)、预测误差的标准偏差(SDEP)等。其他常用内部验证技术还包括不同类型的重采样方法,如自举法(bootstrapping)、随机测试法(randomization tests)以及置换测试(permutation tests)法等。对于利用线性回归分析或PLS分析针对同系物进行QSAR建模可使用这些内部验证方法对模型稳定性进行合理评估。对于线性判别分析建模的SAR模型,只要数据集类别不太失衡,就可使用这些方法对模型稳定性进行评估。相反,对于纯粹非线性方法建立非同系物SAR或QSAR模型,利用上述这些方法进行模型预测能力评估就不太准确。利用纯粹的非线性方法建模,实际上,即使针对相同数据集和参数进行建模,也可能得出不同模型,这些模型的预测略有差异。究竟何种模型最优,这需要折中选择,就是模型对于训练集化学物预测能力和对于外部测试集化学物预测能力之间权重平衡,以确保该模型总体预测能力最优。事实上,对于线性或非线性方法建立的SAR或QSAR模型,若要正确地评估模型预测性能,就必须使用外部测试集。但由于可用的生物数据测试数据有限,以至不具备设置外部测试集条件,因为少之又少的数据都不得已用来建模。在此情况下,只能退而求其次使用内部验证统计方法。

理想情况下,SAR或QSAR模型既可用来作为建模的训练集或学习集,又可作为性能评估的外部测试集。这些集合的选择,则根据具体的建模需求和目标而定。训练集和测试集的划定通常使用线性和非线性多变量方法,以确保这两个集合中的各化学物具备代表生物活性和化学结构的能力。为尽量准确评估结构–活性模型的预测能力,可尝试把测试集拆分为内在测试集(ISTS)和外在测试集(OSTS)。ISTS包含的化学物结构可代表训练集化学结构(如不同位置的同分异构体),用以评估QSAR模型的插值性能。OSTS包含的特殊化学物没必要具有代表性,用以评估模型边缘性或附带性能。很明显,当可用生物数据很多,不成为建模的限制因素,就可使用该策略进行模型性能评估(Lin et al,2003)。

(8)三维QSPR/QSAR基础理论及建模

在Hansch所建立的QSAR模型方法中,描述符的计算是根据化学物分子的“平面”表示图。因此没有必要了解受体-配体的相互作用。而目前QSAR建模已拓展到六维QSAR(6D-QSAR),甚至更高阶形式,但实质上而言,真正将线性平面建模拓展至立体建模且意义重大的是3D-QSAR。Cramer提出的比较分子场分析法(comparative molecular field analysis,CoMFA)认为,如果一组结构类似化学物以同样方式作用同一个受体化学物,它们与受体分子之间在各种作用场应具有一定相似性,而活性取决于每个化学物周围分子场的差异,而分子场在分子水平上影响生物活性的相互作用通常是非共价力,利用分子力场中的立体场和静电场可较好地解释分子的生物活性。目前3D-QSAR建模方法发展较为完善,根据分子场及其相关信息、对齐标准以及化学计量学基础上建立起来的相关方法也较多。

主要3D-QSAR/QSPR方法如下。

①分子场分析法(CoMFA) 分子场分析法基本假设是将一组分子叠合在一起,对其周围的范德华力场和库仑静电场恰当采样,所获取的信息足以解释这组分子的生物活性。采样是通过将分子叠合空间划分为若干固定空间间隔网格,然后用一个探针粒子在分子周围的空间中游走,计算探针粒子与分子之间的相互作用,并记录下空间不同坐标中相互作用的能量值,从而获得每个网格点上,每个分子与一个适当探针之间的相互作用能量。使用PLS分析这些通过计算得出的场参数属性(能量值)与建模研究的生物活性之间的关系。CoMFA关键步骤是初始的分子精准叠合,该步骤既耗费时间又需要有一定前期研究积累。进行分子叠合通常在建模数据集选取最活跃的分子作为模板。对于非刚性分子,活性构象的选择是进行分子叠合的主要障碍之一,因此首先需要进行系统的构象搜索,以定义将要使用的最低能量构象。CoMFA方法作为第二代分子场比较分析定量构效关系模型构建方法具有四大特点:模型生成速度快、全自动进行数据拟合、无需手动叠合分子结构、模型可用于虚拟筛选。

②相似性指数分析法(CoMSIA) 相似性指数分析法(comparative molecular similarity index analysis)是Klebe等于1994年提出,CoMSIA与CoMFA相对,不同之处在于分子场的能量函数采用与距离相关的高斯函数,而非库仑势能和Lemard-Jone势能函数形式。CoMSIA将一组分子叠合并嵌入在三维晶格中(此步骤与CoMFA相同),对分子中原子与一个探针之间的相似性评估使用多达五种场属性,包括立体场、静电场、疏水场、氢键受体场和氢键配体场,因此可更清晰地分析和显示分子叠合空间区域中不同场属性对于生物活性的贡献。

③共同反应性模式方法(COREPA) 为避免分子叠合构象对齐和探针选择(特定的药效团的原子/片段)之间存在误差的问题,出现共同反应性模式方法。该方法通过分析一组化学物分子构象分布的整体和局部反应参数(这些反应性参数与建模研究的终点相关)实现预测。

④距离几何学三维定量构效关系(DG3D-QSAR) 距离几何学三维定量构效关系严格来讲是一种介于二维和三维之间的QSAR方法。这种方法将药物分子划分为若干功能区块定义药物分子活性位点,计算低能构象时各个活性位点之间的距离,形成距离矩阵;同时定义受体分子的结合位点,获得结合位点的距离矩阵,通过活性位点和结合位点的匹配为每个分子生成结构参数,对生理活性数据进行统计分析。分子形状分析认为药物分子的药效构象是决定药物活性的关键,比较作用机理相同的药物分子的形状,以各分子间重叠体积等数据作为结构参数进行统计分析获得构效关系模型。