摘要:本文主要讨论了如何评价虚拟筛选方法的性能,如何计算命中率、富集因子(EF)、绘制ROC曲线以及计算ROC曲线下面积(ROC AUC)。还讨论了为什么要选择合理的数据集进行方法学验证,并介绍了DUD与EDUD这两个社区数据集的利用。

作者:肖高铿

一. 背景

虚拟筛选的终极目的是将活性化合物从海量的数据库中尽可能地富集出来,在大规模虚拟筛选前需要考察一个虚拟筛选方法的性能:只有显著有效的方法我们才敢将之进一步用于大规模虚拟筛选。评估虚拟筛选方法性能常用回溯性验证:收集一个靶点的活性化合物与decoy化合物,进行虚拟筛选,用命中率(hit rate)、富集因子(Enrichment Factor,EF)与ROC曲线(Receiver Operating Characteristic Curve, ROC Curve)来评估虚拟筛选的性能。有的虚拟筛选软件比如Ligandscout与ROCS自带统计学工具评估虚拟筛选性能,大部分软件需要自己去评估统计学性能,因此有必要掌握如何计算EF与绘制ROC曲线。

1. 命中率(Hit rate)

Hit rate = 100% × (Hitssampled/N sampled)

其中N sampled是打分靠前的N个化合物;Hitssampled是打分靠前的N个化合物中真正为活性化合物的数量。

2. 富集因子(Enrichment Factor, EF)

EF = (Hitssampled/Nsampled) ÷ (Hitstotal/Ntotal)

其中Hitstotal为数据库中真正为活性化合物的总数;Ntotal为整个数据库中化合物的总数。

一般认为:当EF值大于1时,这个方法具有显著地活性化合物富集能力。

EF数学形式可以变换为:

EF = (Hitssampled/Nsampled) × (Ntotal/Hitstotal)
= Hit rate × EFMax

EFMax=Ntotal/Hitstotal

对于给定的decoy与ligand数据集,Ntotal/Hitstotal是常数,称为EFMax。不同的测试集其EFMax值不同,其EF值也不同,计算的EF值因此是不具可比性的。

EFX%:因为我们只关心早期的富集能力,尤其是打分最高前0.5%, 1%, 2%的时候EF特别感兴趣,此时EF通常写为EF0.5%,EF1%,EF2%

3. ROC曲线与ROC AUC的绘制原理

就像高效液相分离性质相近的两个组份--目标化合物与杂质--一样,虚拟筛选要分离的两个组份是:活性(ligand)与非活性(inactive,decoy)化合物(如图1所示)。分子对接、药效团筛、分子形状与静电技术以及2D指纹图谱等等之类的虚拟筛选技术就像是柱子、虚拟筛选的参数就像是液相分离的条件,我们的目标是要使用合适的条件将活性化合物与非活性化合物这两个组份充分分离。无论是那种虚拟筛选技术,都是对化合物进行打分,打分值就像保留时间,通过对打分值的统计,我们可以评估一个虚拟筛选方法是否可以有效的将活性化合物从海量的数据库中分离、富集起来。

如何进行虚拟筛选的方法学验证-墨灵格的博客

图1.虚拟筛选就象用高效液相分离性质相近两个组份--目标化合物与杂质--一样,只不过虚拟筛选要分离的两个组份是:活性(ligand,红色)与非活性(inactive/decoy,绿色)化合物。

ROC曲线(ROC Curve)是Operating Characteristic Curve的简称, 其绘制原理如图2所示。当我们根据虚拟筛选的打分值设定截断值,将优于截断值的化合物选出用于进一步活性测试,而将比截断值差的化合物抛弃掉;在这过程中有的阳性化合物被虚拟筛选判为阴性化合物,此时产生假阴性;阳性化合物被正确识别为阳性,称为真阳性;而阴性化合物而言,同理也有假阳性与真阴性。绘制不同截断值下的假阳性率与真阳性率图,即得到ROC曲线(绿色的那根)。

好的虚拟筛选方法,其ROC曲线应该尽可能靠近左上角,同时其曲线下面积(ROC AUC)也越大,理想的方法其ROC AUC=1。靠运气、随机命中方法的ROC曲线应该靠近对角线(虚线),此时ROC AUC=0.5。如果一个虚拟筛选方法接近对角线,说明该方法与靠运气的随机筛选无显著区别,就不值得用于这个项目中。如果一个虚拟筛选方法的ROC曲线位于对角线的下方,说明这个方法只能命中非活性化合物。

如何进行虚拟筛选的方法学验证-墨灵格的博客

图2. ROC曲线绘制的原理

4. 绘制ROC曲线与计算ROC AUC的软件

少部分软件比如Ligandscout与OpenEye的ROCS自带ROC曲线绘制功能,大部分情况下需要自行绘制。免费的统计学软件R,KNIME等具有ROC曲线绘制的专用软件包,还有SPSS等通用软件一般也提供有ROC曲线绘制功能。

二. 实例介绍

2.1. Ligandscout的ROC曲线绘制

具体操作流程见《Ligandscout教程–基于药效团的虚拟筛选》:设立decoy数据集与ligandscout训练集导入虚拟筛选模块的数据库,在虚拟筛选模块点击绘制ROC曲线,虚拟筛选完毕,即可获得ROC曲线与EF值,见图3。
如何进行虚拟筛选的方法学验证-墨灵格的博客

图3. Ligandscout的ROC曲线与EF值计算

2.2. ROCS的ROCS曲线绘制

具体见ROCS教程:在vROCS点击perform a ROCS validation(见图4 A),导入活性与decoy化合物,完成虚拟筛选即可分析活性与decoy的分离性能(图4B)。
如何进行虚拟筛选的方法学验证-墨灵格的博客

图4 A. vROCS的的perform a ROCS validation

如何进行虚拟筛选的方法学验证-墨灵格的博客

图4 B. 验证完毕vROCS报告EF值与ROC曲线

2.3. 用R软件包绘制FRED分子对接的ROC曲线

R软件包提供了多种ROC曲线绘制工具,比如enrichvs,ROCR,pROC等。以enrichvs为例,说明如何绘制ROC曲线。

  1. 数据整理
  2. DUD[3]下载PDE5数据集,用FRED将ligand与decoy的分别进行对接计算(具体方法见分子对接教程–用FRED进行虚拟筛选[7]),将打分结果合并生成一个文本文件,至少包含两列: 一列为打分值(越低越好),一列为标签(1表示为活性化合物ligand的打分,0表示对应这decoy的打分):

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    
    Title score label
    ZINC04199924 -16.976294 1
    ZINC04199924 -16.976294 1
    ZINC04199929 -16.28743 1
    ZINC04199923 -16.262144 1
    ZINC04199923 -16.262144 1
    ZINC01487393 -15.785425 1
    ZINC04199929 -15.64194 1
    ZINC04199920 -15.478016 1
    ZINC00579735 -15.425815 0
    ZINC03808646 -15.39209 0
    ZINC00579735 -15.37167 0
    ZINC00579735 -15.37167 0
    ZINC03811774 -15.346581 0
    ZINC03811774 -15.291958 0
    ZINC04199919 -15.268067 1
    ZINC04199930 -15.24563 1
    ZINC02109268 -15.172167 0
    ZINC04199929 -15.042315 1
    ZINC03811776 -14.932562 0
    ZINC03792330 -14.876609 1
    ZINC03786304 -14.802482 0
    ...
    #此处省略了很多2000多行

    假设上面的内容保存在pde5_score.txt的文件里。

  3. 将数据载入
  4. 在R里键入如下命令:

    1
    
    pde5=read.table("d:/fred/dud/pde5_score.txt",header=T)
  5. 绘制ROC曲线
  6. 在R里键入如下命令:

    1
    
    plot_enrichment_curve(pde5$score, pde5$label,decreasing=FALSE, col="blue")

    结果见图5:

    如何进行虚拟筛选的方法学验证-墨灵格的博客

    图5. FRED对接计算的打分值绘制ROC曲线图(DUD PDE5数据集)

  7. 计算pROC AUC
  8. 计算top 1%的ROC AUC:

    1
    2
    
    auac(score$score,score$label,decreasing=FALSE,top=0.01)</br>
    [1] 0.08435961
  9. 计算ROC AUC
  10. 1
    2
    
    auc(score$score,score$label,decreasing=FALSE)</br>
    [1] 0.7996552

2.4 用python来绘制FRED分子对接的ROC曲线

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
#!/usr/bin/env python
# coding: utf-8
 
import pandas as pd
import numpy as np
from sklearn import metrics
import os
 
# 读入docking后的打分文件
 
df = pd.read_csv('pde5_score.csv')
scores = df.score
y = df.label
 
# 计算假阳性,假阴性与阈值
fpr, tpr, thresholds = metrics.roc_curve(y, -scores, pos_label=1)
 
 
# 计算ROC AUC
auc = metrics.auc(fpr, tpr)
print('AUC ROC:',auc)
 
#绘图
import matplotlib.pyplot as plt
plt.plot(fpr,tpr)
plt.plot([0, 1], [0, 1], color="darkorange", linestyle="--",label="ROC curve (area = %0.2f)" % roc_auc)
plt.ylabel('True Positive Rate')
plt.xlabel('False Positive Rate')
plt.title("Receiver operating characteristic")
plt.legend(loc="lower right")
plt.show()
如何进行虚拟筛选的方法学验证-墨灵格的博客

图6 用Python绘制的ROC曲线

三. 讨论:为什么要设计decoy数据集来进行方法学验证

以分子对接为例,分子对接常用于:虚拟筛选发现苗头化合物、预测结合模式、根据打分值进行先导化合物物优化。为了考察分子对接打分能否可靠的用于先导物优化,Enyedy等人(2008)考察了打分值能否区分活性化合物与非活性化合物。结果发现,就ROC AUC而言,分子量与ClogP跟GLIDE的GLIDE_XP与GLIDE_SP一样具有同样的活性与非活性化合物区分能力,见图7A,B。
如何进行虚拟筛选的方法学验证-墨灵格的博客

图7 A. KDR(1YWN)对接计算的ROC曲线

如何进行虚拟筛选的方法学验证-墨灵格的博客

图7 B. KDR(1YWN)对接计算的ROC曲线

由此可见,我们不能仅根据分子对接与打分值来指导先导的结构优化,专家的解释还是非常重要的;对接软件的打分值与可能化合物的物理性质诸如分子量、极性表面积、油水分配系数、柔性键数量、氢键受体数、氢键供体数有很强的相关性。如果活性化合物与数据库化合物在这些物理性质上很大的差异,那么虚拟筛选性能的评估结果就会很受质疑:很可能是物理性质自身差异的原因导致化合物被富集出来,而不是分子对接或打分函数本身。

同时,也说明合理的验证用数据集很重要:首先这个数据集应该消除化合物物理性质带来的偏差。Huang等人[3]开发的DUD基准数据集正是专为虚拟筛选方法验证而设计。DUD数据集在使用过程中也发现一些缺陷,并被不断改进。DUDE[4]不仅是DUD的曾强版,还可以帮助我们生成活性化合物的decoy数据集。还有其它的数据集可供使用,比如Bauer等人[5]开发的DEKOIS 2.0数据集,以及最近Cleves等人[6]的DUDE+数据集等。

四. 附件

PDE5数据集FRED打分的结果:pde5_score.csv

五. 文献

  1. Yabuuchi, H., et al. (2011). "Analysis of multiple compound-protein interactions reveals novel bioactive molecules." Mol Syst Biol 7: 472.
  2. Enyedy, I. J.; Egan, W. J. Can We Use Docking and Scoring for Hit-to-Lead Optimization? Journal of computer-aided molecular design 2008, 22 (3-4):161–168.
  3. Huang, N.; Shoichet, B. K.; Irwin, J. J. Benchmarking Sets for Molecular Docking. J. Med. Chem. 2006, 49 (23), 6789–6801. https://doi.org/10.1021/jm0608356.
  4. Mysinger, M. M.; Carchia, M.; Irwin, J. J.; Shoichet, B. K. Directory of Useful Decoys, Enhanced (DUD-E): Better Ligands and Decoys for Better Benchmarking. J. Med. Chem. 2012, 55 (14), 6582–6594. https://doi.org/10.1021/jm300687e.
  5. Bauer, M. R.; Ibrahim, T. M.; Vogel, S. M.; Boeckler, F. M. Evaluation and Optimization of Virtual Screening Workflows with DEKOIS 2.0 – A Public Library of Challenging Docking Benchmark Sets. J. Chem. Inf. Model. 2013, 53 (6), 1447–1462. https://doi.org/10.1021/ci400115b.
  6. Cleves, A. E.; Jain, A. N. Structure- and Ligand-Based Virtual Screening on DUD-E + : Performance Dependence on Approximations to the Binding Pocket. J. Chem. Inf. Model. 2020, 0 (0), acs.jcim.0c00115. https://doi.org/10.1021/acs.jcim.0c00115.
  7. 肖高铿.分子对接教程–用FRED进行虚拟筛选. 墨灵格的博客. 2016.-09-15. http://blog.molcalx.com.cn/2016/09/15/openeye-tutorial-fred.html