虚拟筛选及其方法学验证与ROC曲线python实现

前言

分子对接代表的基于结构的虚拟筛选已成为基于结构药物设计的标准方法。下面是有关虚拟筛选的笔记，另外一部分是机器学习二分类模型中常用的评估指标ROC曲线及AUC值，有错误请留言指出。
附注：主要从维基百科和百度百科的资料。

---

虚拟筛选

虚拟筛选(virtual screening，VS)也称计算机筛选，即在进行生物活性筛选之前，利用计算机上的分子对接软件模拟目标靶点与候选药物之间的相互作用，计算两者之间的亲和力大小，以降低实际筛选化合物数目，同时提高先导化合物发现效率。

从原理上来讲，虚拟筛选可以分为两类，即基于受体的虚拟筛选和基于配体的虚拟筛选。

基于受体的虚拟筛选从靶蛋白的三维结构出发，研究靶蛋白结合位点的特征性质以及它与小分子化合物之间的相互作用模式，根据与结合能相关的亲合性打分函数对蛋白和小分子化合物的结合能力进行评价，最终从大量的化合物分子中挑选出结合模式比较合理的、预测得分较高的化合物，用于后续的生物活性测试。相对基于配体的虚拟筛选，其优势可避免活性化合物结构微小改变导致的活性改变。不足之处：最大的问题是打分函数的准确性和适用性，一般考虑计算机计算速度，通常采用比较简单的打分函数，但简单的打分函数不能很好的考虑弱相互作用。其次是基于结构的虚拟筛选需要受体结构和指定结合位点，但许多重要的靶标都没有可用的受体结构。而且分子对接中仍存在需要解决的问题是：1.分子柔性，2.溶剂化效应，3.打分函数。

基于配体的虚拟筛选依据的是结构决定性质的原理，一般是利用已知活性的小分子化合物，根据化合物的形状相似性或药效团模型在化合物数据库中搜索能够与它匹配的化学分子结构。最后对这些挑选出来的化合物进行实验筛选研究。此类方法也有三种模型：1.药效团模型，2.定量构效关系，3.结构相似性方法。作为分子对接的重要补充手段,基于配体的虚拟筛选具有速度快,通用性好的优点。

分子对接是一种基于靶标蛋白结构的药物筛选方式。通过小分子化合物与靶标进行分子对接，综合分析得分及空间构象情况，包括静电作用、氢键作用、疏水作用、范德华作用等性质，可以探索配体小分子与受体生物大分子具体作用方式和结合构型，解释化合物产生活性的原因，为合理地优化化合物结构提供指导；也可筛选潜在活性化合物，为实验提供参考。

药效团筛选是一种基于小分子化合物的高效药物筛选手段。通过分析一个或多个活性小分子的药效特征，推纳概括出使得分子具有活性的重要药效基团特征。药效团筛选的计算量较小，可以在分子对接前进行，对几百万或几千万的小分子数据库进行药效团筛选只需要很短时间。客户只需要提供1个或多个活性分子，就可以构建公共药效团进行筛选，搜索含相同特征的小分子，并指导新活性分子的合成。

定量构效关系(quantitative structure activity relationship ,QSAR)，是一种借助分子的理化性质参数或结构参数，以数学和统计学手段定量研究有机小分子和生物大分子相互作用，有机小分子在生物体内吸收、分布、代谢、排泄等生理相关性质的方法。

结构相似性方法(structual similarity, SSIM）是通过各种描述符或指纹进行相似性匹配，从而判断化合物是否具有类似活性或治病机理。

下面内容是虚拟筛选中涉及到的相关方法，摘自《如何进行虚拟筛选的方法学验证》http://blog.molcalx.com.cn/2016/09/22/virtual-screening-methodology-validation.html

一、虚拟筛选方法学

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

1. 命中率(Hit rate)

Hit rate = 100% × (Hits^sampled/N ^sampled)

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

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

EF = (Hits^sampled/N ^sampled) ÷ (Hits^total/N^total)

其中Hits^total为数据库中真正为活性化合物的总数；N^total为整个数据库中化合物的总数。

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

EF数学形式可以变换为：

EF = (Hits^sampled/N ^sampled) × (N^total/Hits^total)
= Hit^rate × EF^Max

EF^Max=N^total/Hits^total

对于给定的decoy与ligand数据集，N^total/Hits^total是常数，称为EF^Max。不同的测试集其EF^Max值不同，其EF值也不同，计算的EF值因此是不具可比性的。

EF^X%：因为我们只关心早期的富集能力，尤其是打分最高前0.5%, 1%, 2%的时候EF特别感兴趣，此时EF通常写为EF^0.5%，EF^1%，EF^2%。

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

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

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

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

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

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

以分子对接为例，分子对接常用于：虚拟筛选发现苗头化合物、预测结合模式、根据打分值进行先导化合物物优化。为了考察分子对接打分能否可靠的用于先导物优化，Enyedy等人（2008）考察了打分值能否区分活性化合物与非活性化合物。结果发现，就ROC AUC而言，分子量与ClogP跟GLIDE的GLIDE_XP与GLIDE_SP一样具有同样的活性与非活性化合物区分能力。


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

同时，也说明合理的验证用数据集很重要：首先这个数据集应该消除化合物物理性质带来的偏差。DUD基准数据集(http://dud.docking.org)正是专为虚拟筛选方法验证而设计，DUDE(http://dude.docking.org)可以帮助我们生成活性化合物的decoy数据集：一组与活性化合物具有相似物理性质的随机化合物。
decoy在分子对接和虚拟筛选的性能评估有影响，需注意。

二、ROC曲线及AUC的python实现

因为ROC曲线和AUC是机器学习二分类模型中非常常用的评估指标，刚开始理解比较吃力，建议结合例子进行理解。下面面容参考《机器学习评价指标 ROC与AUC 的理解和python实现》https://blog.csdn.net/qq_20011607/article/details/81712811?utm_medium=distribute.pc_relevant_t0.none-task-blog-BlogCommendFromMachineLearnPai2-1.channel_param&depth_1-utm_source=distribute.pc_relevant_t0.none-task-blog-BlogCommendFromMachineLearnPai2-1.channel_param

首先一个用分类器预测完后，我们会得到一个二分类的混淆矩阵：

准确率
所有预测正确的样本占所有样本的比例，但这个指标在类别比例不平衡时不太恰当。

查准率(Precision)
直观理解为： 在所有预测为正例中，真正正例的比例。（有点绕，多读几遍就好）

查全率(Recall)
直观理解为： 在所有真实的正例中，预测为正例的比例。

F1-score
由于Precision和Recall是一对不可调和的矛盾，很难同时提高二者，也很难综合评价。故提出F1来试图综合二者，F1是P和R的调和平均（harmonic mean）

ROC曲线：很多学习器能输出一个实值或者概率预测，然后设定一个阈值，高于阈值为正类，反之负类。分类的过程就是设定阈值，并用阈值对预测值做截断的过程，当这个阈值发生变动时，预测结果和混淆矩阵就会发生变化，最终导致一些评价指标的值的变化。

在进一步了解ROC之前，首先定义两个指标：

真正例（True Positive Rate，TPR），表示所有正例中，预测为正例的比例：

假正例（False Positive Rate，FPR），表示所有负例中，预测为正例的比例：

以FPR为横坐标，TPR为纵坐标，那么ROC曲线就是改变各种阈值后得到的所有坐标点 (FPR,TPR) 的连线，画出来如下。红线是随机乱猜情况下的ROC，曲线越靠左上角，分类器越佳。

真实情况下，由于数据是一个一个的，阈值被离散化，呈现的曲线便是锯齿状的，当然数据越多，阈值分的越细，”曲线”越光滑。
AUC（Area Under Curve）就是ROC曲线下的面积。AUC值是一个概率值，当你随机挑选一个正样本以及一个负样本，当前的分类算法根据计算得到的Score值将这个正样本排在负样本前面的概率，就是AUC值。当然，AUC值越大，当前的分类算法越有可能将正样本排在负样本前面，即能够更好的分类。

用AUC判断分类器（预测模型）优劣的标准：

AUC = 1 是完美分类器，采用这个预测模型时，存在至少一个阈值能得出完美预测。绝大多数预测的场合，不存在完美分类器。
0.5 < AUC < 1，优于随机猜测。这个分类器（模型）妥善设定阈值的话，能有预测价值。
AUC < 0.5，比随机猜测还差；但只要总是反预测而行，就优于随机猜测。

简单说：AUC值越大的分类器，正确率越高。

---

用python的sklearn：绘制ROC曲线+求AUC

求 ROC：

import numpy as np
from sklearn.metrics import roc_curve, auc

y = np.array([1, 1, 1, 1, 1,
              2, 2, 2, 2, 2])
y_proba = np.array([0.1, 0.4, 0.4, 0.3, 0.5,
                    0.4, 0.6, 0.7, 0.8, 0.5])   # probability of prediction as positive

fpr, tpr, thresholds = roc_curve(y, y_proba, pos_label=2)
auc = auc(fpr, tpr)

运行后在窗口控制台分别键入fpr，tpr，thresholds，auc：

>>> fpr
array([0. , 0. , 0.2, 0.6, 1. ])
>>> tpr
array([0.2, 0.6, 0.8, 1. , 1. ])
>>> thresholds
array([0.8, 0.6, 0.5, 0.4, 0.1])
>>> auc
>>> 0.9

将 fpr 和 tpr 的值用pylab绘制即得ROC曲线。如下图：

整个ROC曲线及AUC值运行代码如下：

from sklearn import metrics
import numpy as np
from sklearn.metrics import roc_curve, auc
import pylab as plt

y = np.array([1, 1, 1, 1, 1,
              2, 2, 2, 2, 2])
y_proba = np.array([0.1, 0.4, 0.4, 0.3, 0.5,
                    0.4, 0.6, 0.7, 0.8, 0.5])   # probability of prediction as positive
fpr, tpr, thresholds = roc_curve(y, y_proba, pos_label=2)
auc = auc(fpr, tpr)

roc_auc = metrics.auc(fpr, tpr)
plt.figure(figsize=(6,6))
plt.title('Validation ROC')
plt.plot(fpr, tpr, 'b', label = 'Val AUC = %0.3f' % roc_auc)
plt.legend(loc = 'lower right')
plt.plot([0, 1], [0, 1],'r--')
x = [0, 1] #plt.xlim([0, 1])
y = [0, 1] #plt.ylim([0, 1])
plt.ylabel('True Positive Rate')
plt.xlabel('False Positive Rate')
plt.show()