バーチャルスクリーニングの指標のテストその２

昨日の記事の続きです．なんの話かということは昨日の記事を参照のこと．
(バーチャルスクリーニングの指標のテスト)

今日はさらに指標を４つ紹介する．AU-ROC，AUAC，BEDROC，EFの４つ．

AU-ROC (Area Under the Receiver Operating Characteristic Curve)
みんな大好きROC曲線のAUC．
バーチャルスクリーニング（early recognition problem）でのROC曲線は，横軸にヒットじゃない化合物の数，縦軸にヒット化合物の数を引いてプロットする．
例えば「リガンドの数10，ヒットの数5，ヒットの順位は1,3,4,6,9だった」ときはこんな感じ．

で，この曲線の下の面積がAU-ROC (AUCと呼ばれることの方が多いが，) である．

で，このAU-ROCはまぁ一応積分計算なのであるが，離散問題なので簡単に計算できる．
${\mbox{AU-ROC}=1-\frac{\sum_{i=1}^n r_i}{n(N-n)}+\frac{n+1}{2(N-n)}}$
${n}$ はヒットの数， ${r_i}$ は( ${i}$ 個目の)ヒットの順位， ${N}$ はテストした化合物の数である．

有名でよく使われる指標だが，バーチャルスクリーニングの指標としてはあまり良くない，と言われている*1．

AUAC (Area Under the Accumulation Curve)
あんまり有名じゃないしあんまり使われてない指標．ROC曲線は横軸にヒットじゃない化合物の数をとっていたが，Accumulation曲線は横軸に（全ての）化合物の数をとる．
例えばさっきの例，「リガンドの数10，ヒットの数5，ヒットの順位は1,3,4,6,9だった」はこんな感じ．

Accumulation曲線のときは線は斜めに引く．この曲線下面積AUACも，もちろん簡単に計算できる．
${\mbox{AUAC}=1-\frac{1}{nN}\sum_{i=1}^n r_i}$
${n}$ はヒットの数， ${r_i}$ は( ${i}$ 個目の)ヒットの順位， ${N}$ はテストした化合物の数である．
やっぱりバーチャルスクリーニングの指標としてはあまり良くない，と言われている．
なお， ${n << N}$ だと AU-ROC≒AUAC となる．そういうわけでAU-ROCを使ってればOK，なのだ．

BEDROC (Bolzmann-Enhanced Discrimination ROC)
昨日紹介したRIEを0-1に正規化した感じの値．
${\mathrm{BEDROC}=\mathrm{RIE}\times\frac{\frac{n}{N}\sinh(\alpha/2)}{\cosh(\alpha/2)-\cosh(\alpha/2-\alpha n/N)}+\frac{1}{1-e^{\alpha(N-n)/N}}}$
$\alpha$ はRIEで使っているパラメータである．
BEDROCは結構良い指標らしい．最近は良く使われている．

EF (Enrichment Factor)
この界隈ではかなり有名．結構良い指標と言われている．
${\mathrm{EF}=\frac{\mathrm{count}(\mathrm{hits} | r_i \leq \chi N)}{\chi n}}$
${\chi}$ はEFのパラメータで，上位 ${\chi}$ %を見て計算する，という意味になる．
＃ちなみにmimeTeX記法でなぜか"#"が入力できなかった．"\#"でも駄目らしい．countは#と同じ意味．

例えば「200個のhitが10000個のligandの中にある．15%(30個)のhitが上位5%(500位以内)に存在した．」というときのEFは，
${\chi=0.05}$ ， ${\mathrm{count}(\mathrm{hits} | r_i \leq \chi N)=30}$ なので，
${\mathrm{EF}=30/(0.05\times 200)=3}$
となる．

EFの最大値は ${1/\chi}$ ，ランダムのときはほぼ1で正確には ${\lfloor \chi N\rfloor/\chi N}$ となる．

４つの指標をテスト
前回のプログラムをちょっと変更して今回紹介したパラメータを計算できるようにした．

#!/usr/bin/python

import sys
from math import *

if len(sys.argv) != 3:
    print 'Usage: %s [N] [list of hits (ex: 2,3,5)]' % sys.argv[0]
    quit()

N = float(sys.argv[1])
hits = sys.argv[2].split(",")
n = float(len(hits))

ec = 0 # EF count
arp = 0.0
rie = 0.0

alpha = 8.0 # RIE paramater
chi = 0.2 # EF parameter

for ri in range(1, int(N)+1):
    if str(ri) in hits:
        arp += ri
        rie += exp(-alpha*ri/N)
        if ri <= chi*N:
            ec += 1
    ri += 1

arp = arp / n
rie = rie / ((n/N)*((1-exp(-alpha))/(exp(alpha/N)-1)))
auroc = 1 - arp/(N-n) + (n+1)/(2*(N-n))
auac = 1 - arp/N
bedroc = rie * (n/N) * sinh(alpha/2)/(cosh(alpha/2)-cosh(alpha/2-alpha*n/N)) + 1/(1-exp(alpha*(N-n)/N))
ef = float(ec)/(chi*n)

print "ARP =", arp
print "RIE =", rie
print "AU-ROC =", auroc
print "AUAC =", auac
print "BEDROC =", bedroc
print "EF =", ef

今回はリガンドを10個，ヒットが4個で試してみる．
Aのヒット順位は(1, 2, 5, 9)，Bのヒット順位は(2, 3, 5, 6)だったとしよう． ${\chi=0.2, \alpha=8.0}$ で計算する．

$ python virscr2.py 10 1,2,5,9
ARP = 4.25
RIE = 2.05435170462
AU-ROC = 0.708333333333
AUAC = 0.575
BEDROC = 0.855180830287
EF = 2.5

$ python virscr2.py 10 2,3,5,6
ARP = 4.0
RIE = 0.978186814512
AU-ROC = 0.75
AUAC = 0.6
BEDROC = 0.402850472687
EF = 1.25

というわけで，

ARP：Bの勝ち
RIE：Aの勝ち
AU-ROC：Bの勝ち
AUAC：Bの勝ち
BEDROC：Aの勝ち
EF：Aの勝ち

である．指標によって変わるが，ARPやAU-ROC, AUACは少数でも後ろの方の順位のヒットがあると引きずられやすいため，Bの方が良くなっている，と言えるだろうか．逆にRIE, BEDROC, EFは，パラメータを使って上位にヒットがあることを考慮して値を出しているためAの方が良いという結果になると言えると思う．
おとなしくEF（かBEDROCあたり，もしくはもっと新しいもの）を使いましょうということだ．

*1:Truchon J-F, Bayly CI. Evaluating virtual screening methods: good and bad metrics for the "early recognition" problem. J Chem Inf Model, 47: 488-508, 2007.　とか