スコア関数を変えてドッキング結果を再解析する話

さて、前回の記事でドッキングによるSBVSもなんとか終わったので後は化合物IDのリストを提出するだけ・・・、というところですがとても気になる文献を知ってしまいました。

Performance of machine-learning scoring functions in structure-based virtual screening. Wójcikowski, M., Ballester, P.J., and Siedlecki, P., Sci.Rep. 2017(7)46710 pubmed

こちらの文献は前回使用したODDT(Open Drug Dicovery Toolkit)と同一の筆頭著者によるもので、どうやらドッキングのスコアを機械学習を用いて改善したという内容のようです。しかもODDT同様に、このスコア関数（RF-ScoreVS）を利用できるとのこと！！

前回の記事で参照した「タンパク質計算科学基礎と創薬への応用」*1 において、ドッキングのスコア関数は依然として精度に課題がある（実験値との乖離がある）といった趣旨の記述がありました。AutoDock VinaはAutoDockよりも精度が改善している、とはいえ10年ほど前のソフトウェアです。より改善したものがあるのであれば使ってみたい！

（・・・というのは建前でMachine Learningって言って格好つけたかった）

1. RF-Score-VS （Sci.Rep. 2017(7)46710）文献紹介
2. RF-Score-VS を使ってみる
3. アッセイ候補提出ファイルの作成
4. まとめ
5. 総括

1. RF-Score-VS （Sci.Rep. 2017(7)46710）文献紹介

まずは文献を読んでみようと思います。Scientific ReportsはOpen Accessなのが良いですね。(CC BY 4.0)
すでに記憶があやふやなのでスコア関数について少し復習します。

1-1. スコア関数の種類

スコア関数にはどのようなものがあるか？
先の記事で参照した日本化学会情報化学部会誌の文献（タンパク質-リガンドドッキングの現状と課題 *2 によると、大きく３つ

①力場に基づいたスコア関数（force-field-based scoring function）
②経験的なスコア関数 (empirical scoring function)
③ 統計に基づいたスコア関数(knowledge-based scoring function)

があり、さらに近年の動きとして、

④複数のスコア関数を利用するコンセンサススコア(concensus scoring)
⑤機械学習に基づいたスコア関数(machine-learning based scoring function)

が開発されている、とのことでした。

「④ コンセンサススコア」について「タンパク質計算科学」の説明を参照します。
「コンセンサス（consensus、合意）」をとる目的は、単一の手法の結果のみに頼らず、複数のスクリーニング結果を統合することでヒット率の向上、当たり外れの低減を目指すことです。結果を統合する方法として、各化合物に対して複数の手法でスコアを得たのち

積集合: 複数の手法で共通に選ばれた化合物をヒットとする
和集合: 複数の手法で少なくとも一つに選ばれた化合物をヒットとする
rank by scoring: 各手法のスコアを「0~1」の範囲にスケールした上で和をとり、新しいスコアとする
rank by rank: 各手法における化合物の順位（ランク）の平均によってランキングする

といった手法が例に挙げられていました。

では、「⑤ 機械学習に基づいたスコア関数」はどのようなものか？　論文の中身に進みます。

1-2. 既存のスコア関数の問題点

関数の構造上の問題

そもそもなぜ機械学習の導入にいたったか？論文中では従来から用いられているスコア関数を"classical scoring functions (SFs)" （おそらく「①力場」や「②経験的」、「③統計」によるスコア関数のこと？）としており、結合親和性の予測精度の向上に限界が見られていることの理由として、以下の問題点を指摘しています。

コンフォメーションのエントロピーをうまく説明しない
溶媒和によるエネルギーの寄与をうまく説明しない

また、従来のスコア関数は線形回帰モデル(linear regression model)を下敷きにしているため、

大量のタンパク質構造や結合に関する情報を吸収してモデルに取り込むことに限界がある(universal SFの限界)
標的のタンパク質（protein family）に焦点を合わせてスコア関数のパラメーターをチューニングしたくても、回帰モデルが変更できる形では提供されていない

といった問題を指摘しています。・・・難しい。

よくわからないので別の文献のFigureを引用します。（WIREs Comput. Mol. Sci.2015(5)405）*3

f:id:magattaca:20190407180233p:plain

Fig 中に示されているようにClassocal SFでは線型結合の式ですが、Machine-learnig SFでは関数となっています。そもそもモデル化の構造が異なっているということで、異なるパフォーマンスが期待できるよね！ということのようです。

モデル作成上の問題

さらに、classical SFのモデル作成時の問題点として、トレーニングセットとテストセットに重複(overlap)があるということが挙げられています。機械学習であればcross validationにより、より厳密なバリデーションの取り扱いが可能になります。

また、classical SFは結晶構造の情報を元に作成されていますが、これは結合するというプラスのベクトルの事例のみの情報を用いて構築していることになります。一方でスコア関数を適用するVirtual Screeningの目的は、スクリーニング結果として得られた「実際に結合する正例(プラスのベクトル/ +ve）」と「結合しない負例（マイナスのベクトル/ -ve）」が入り混じった状態である一連の結合ポーズの中から、正しいtrue binderを見つけ出すということです。つまり、Virtual Screeningでは負例（inactiveのリガンド / -ve）の比率の方が高いのに、その情報が関数作成段階に考慮されていないと指摘されていました。

（なるほど確かに言われてみればそう・・・なのでしょうか？勉強不足のため十分に理解できていません。）

ということで、機械学習を導入することで、これらの問題点を解決して行きましょう！　

1-3. 著者らによる機械学習スコア関数（RF-Score シリーズ）

文献中では機械学習スコア関数の例としてRF-Score、NNScore、SFCscoreといった名前が挙げられています。このうち、RF-Scoreは今回の文献の著者の一人、De. Pedro Ballester（CRCM）が筆頭著者として2010年に報告しているものです*4。RFはRandom Forestの略のようですが、こちらをversion 1 （RF-Score v1）として、RF-Score v2*5が報告されています。特にRF-Score v3の文献はタイトルが「Improving AutoDock Vina Using Random Forest: The Growing Accuracy of Binding Affinity Prediction by the Effective Exploitation of Larger Data Sets.」、とAutoDock Vinaを使った以上とても気になる文献です（アクセスできなかった・・・）。

これらのスコア関数の提案に対する議論、「overfittingしているのでは？」「新しいターゲットにも適用できるの？(applicablity)」に対して、非常に大きなデータセットを使って学習を実施し、virtual screeningのパフォーマンスが上げられることを示した、というのが2017年のSci. Rep.の概要となります。

1-4. RF-Scoreの仕組み

では、ベースとなったRF-Scoreはどのようなものなのか？ちょうど良いスライド資料があったので参照していきたいと思います*6。（元文献を読む能力がなくてすみません）

RF-Scoreではランダムフォレストを使ってスコア関数をモデル化していますが、元々の発想としては"classical SF"のうちの「3. 統計に基づいたスコア関数(knowledge-based scoring function)」を改善する、ということから始まっているようです。

統計的スコア関数の流れは以下のようになっています。

①多数の結晶構造（PDBのタンパク質-リガンド複合体構造(P-L complex)）を比較
②原子（あるいは官能基/フラグメント）にタイプを割り当てる
③中心の原子と周囲の原子のtype-type間の距離のヒストグラムを作成する（存在確率の空間分布）
④存在確率（P）を結合自由エネルギー（ΔF）に変換する（ΔF = -k_BT lnP）
⑤全てのP-Lアトムペアのエネルギーを合計して全体のエネルギーとする。

上記の「④ヒストグラムのエネルギー関数への変換」では"reverse Boltzmann"の手法が使われています。ここで、「温度Tの平衡状態においてタンパク質とリガンドの各原子は独立な粒子である」という仮定が入りますが、これは仮定としては粗すぎる、というのが問題となるようです。
理由として

分子中の原子と原子は結合している（Molecular connectivity）から、原子と原子の距離は独立とはいえない
排除体積の影響（Excluded volume）
平衡状態の仮定に物理的根拠がない
温度(T)による構造変化は小さく、ボルツマン分布により暗に考慮されているほどではない

といったことが挙げられています。以上をふまえて、著者らは「ヒストグラム（距離の分布 / distance distribution）と結合自由エネルギーとを関係付ける」ために、Boltzmannの公式のような公式を仮定するのではなく、「構造情報と結合親和性の関係性を機械学習を使って学習すれば良い」と考えたのがRF-Scoreの出発点となっているそうです。

上記の考えに基づき、以下の操作を行って作成したモデルがRF-Score (v1) となります。

①PDBの構造情報(PDBbind Benchmark)と質の高い結合親和性データ（K_d, K_i）
②タンパク質とリガンドの、閾値以下(d_cutoff)の距離にある原子をみつける
③アトムタイプペアをフィーチャーとして出現率を求める
④各タンパク質-リガンド複合体をフィーチャーのベクトルとして表す
⑤ベクトルを入力（説明変数）、結合親和性を出力（目的変数）としてランダムフォレストによる機械学習を実施

加えてさらに、スコア関数の精度向上にむけて公共データと様々な製薬企業の企業内データを用いてスコア関数の精度向上を目指したコンソーシアム（Scoring Function Consrtium）のデータに基づきフィーチャーの選択を改良した（RF-Score v2 / Data-driven feature selection）や、AutoDock Vinaのフィーチャーを取り入れる(RF-Score v3)、といった改良を加えてモデルの精度を高めてきた、というのが流れとなります。

1-5. RF-Score-VSでの取り組み

先に述べたようにRF-Score-VSでは、学習のデータセットに工夫を加えることでVirtual Screeningにおけるスコアの改善が目指されています。

データセットの工夫

まず、データセットとしてDUD-E(A Database of Useful Decoys: Enganced)データセットが使われています。DUD-E*7はスクリーニングソフトの性能テスト（ベンチマーク）として利用されているDUDの改良版であり、DUDは市販化合物の３次元データベースであるZINCのデータを既知の標的タンパク質ごとに、その既知の活性化合物（active）と、類似した化合物（inacitive）との組み合わせとして整理したデータセットであるとのことです*8。

先述のように、RF-Scoreの学習セットとしてはPDBの情報が用いられていました。これに対してRF-Score-VSの利用している学習セットには

より大きなデータセット
結晶構造ではなく、active リガンド(+ve)とinactive リガンド(-ve)のドッキングポーズ

を利用するという違いがあります。結晶構造はactiveリガンドのポーズしか含まないので、こちらを学習セットとして用いるよりも、activeとinactive両方の結合ポーズを学習データとして用いる方がVirtual Screeningに用いるスコア関数としてはより実際の使用状況に類似しているだろう、というのが背景にあるようです。

学習データ分割方法の工夫

また、興味深いことに学習データセットのcross-validationを行うにあたって、

①標的タンパク質ごと（per-terget）
②訓練データ、テストデータともに全ての種類のタンパクを含む（horizontal split）
③訓練データとテストデータで、標的タンパク質が重複しないようにした上で分割する（vertival split）

という３種類の分割方法が用いられています。

論文中の図を引用します。

f:id:magattaca:20190407181127p:plain

②horizontal splitは「すでに活性化合物（active ligand）が知られている標的についてvirtual screeningを行う」という状況を、③vertical splitは「既知の活性化合物がない状態でvirutal screeningを行う」という状況を模倣する、というのを目的としているそうです。

Virtual Screeningが行われる状況の場合分けを想定して、訓練データとテストデータの分割方法を分けるというのは、なるほど〜という感じでとても面白いと思いました（・・・よく使われる方法なのでしょうか？）。

1-6. RF-Score-VSの性能

評価指標

評価指標（metrics）としてはEnrichment factor(EF)というものが用いられています。

$$ EF_{x\%} = \frac{\frac{ActiveMolecules(x\%)}{Molecules(x\%)}}{\frac{ActiveMolecules(100\%)}{Molecules(100\%)}} $$

ランク付けした化合物の上位x%を取り出して、その中に本当の活性化合物が含まれている割合を求めます。この割合と、全化合物（x = 100）中に含まれる活性化合物の割合の比をとったものがEFとなります。Virtual Screeningの性能指標としてよく用いられるものには、EF以外にROC（Receiver Operating Characteristc）曲線のAUC（Area Under the Curve）があるそうです。
EFを用いた理由として、Virtual Screeningを行って選んだ化合物を実際にアッセイするとしたら、上位数パーセントを選んで行うことになる。その数パーセントに"当たり"が含まれていることが大事なので、その状況をより反映しているEFの方が適している、といった説明がなされていました。

確かに創薬レイドバトルも「スコアでランキングをつけて上位10個をアッセイ候補として提出」ということになっていました。この論文とても勉強になる。

性能

ではAutoDock Vinaをドッキングに用い、得られた結合ポーズを各種スコア関数でスコアリングして比較した結果を引用します。横軸が Enrichment Factor 1% （EF_1%）で、縦軸にスコア関数が並んでいます。

f:id:magattaca:20190407181224p:plain

分割方法として② horizontal splitを用いた、RF-Score-VSに置いて非常に良い結果が得られており、従来のclassical SFと比較して高いEF1%となっています。① per-targetは平均値としては、性能が向上しているものの、標的によって性能のばらつきが大きいという結果になっています。また、③ vertical splitは ② horizontalよりも性能は落ちていますが、classical SFよりも性能が向上していると記載されていました。

AutoDock Vinaのスコア関数を改善したという「RF-Score v3」が上記の結果では「AutoDock Vina Score」と比べて良くなっているように見えないのが気になる点ですが、評価対象のデータセットの質が変わった（結晶構造　→　ドッキングポーズ）ことが理由なのかもしれません。

他のドッキングソフト(Docking engine)を使った比較結果、DUD-Eとは異なるデータセット DEKOIS 2.0 に対して適用した際の結果についても記載されていましたので、ご興味のある方は論文を参照していただければと思います。

要するに、とにかく良い結果が出るっぽい！！！　使いましょう！

2. RF-Score-VS を使ってみる

2-1. ドッキングスコアを再計算

とても有り難いことにRF-Score-VSに関して、データセット（ http://github.com/oddt/rfscorevs ）、学習済みのモデル（http://github.com/oddt/rfscorevs_binary ）ともにGitHubで公開されています。

各OS向けの学習済みのbinaryが用意されているので、ready-to-use RF-Score-VSをローカルに落として来ました。
こちらのbinaryにはRD-SCORE-VS v2 が入っている、とのことです。

使い方は簡単で、ターミナルでダウンロードして来たフォルダ（rd-score-vs_v1）に移動します。

以下でヘルプを参照できます。

$ ./rf-score-vs --help

ドッキングに使ったタンパク質のpdbファイル（今回は5nix_chainAB.pdbという名前）と、AutoDock Vinaのドッキングで得たリガンド構造を含むsdf（ODDT_Dock_All_Result_mod.sdfという名前）を同じフォルダに入れてコマンドを実行します。

./rf-score-vs_v1 --receptor 5nix_chainAB.pdb ODDT_Dock_All_Result_mod.sdf -O ODDT_All_ligands_rescored.sdf

引数 --receptor に続いてpdb、sdfを入力し、 -O　で指定したファイル名で、再計算されたスコア（フィールド名: RFScoreVS_v2）が書き込まれたsdfが出力されました。

今回は全てのドッキングポーズ（50化合物ポーズ合計449個）についてスコアを再計算しています。では結果を確認していきましょう。

2-2. スコアを確認

RDKitでSDFを読み込んでおきます。

from rdkit import Chem
from rdkit.Chem import AllChem, PandasTools 
import pandas as pd

#PandasToolsで読み込み、スコアのTypeをfloatに変換する
Rescored_df = PandasTools.LoadSDF('./ODDT_All_ligands_rescored.sdf')
Rescored_df['RFScoreVS_v2'] = Rescored_df['RFScoreVS_v2'].astype(float)
Rescored_df['vina_affinity'] = Rescored_df['vina_affinity'].astype(float)

f:id:magattaca:20190407181604p:plain

AutoDock Vinaのスコアは結合自由エネルギ-の予測値（kcal/mol）で与えられるのに対し、RF-Score-vsのスコア再計算結果は pK_d/i の単位で与えられます。
学習の段階でDUD-Eの化合物の活性を pK_d/i = 6 よりも弱い（小さい）場合にinactiveとしています。また、トレーニングセットのDecoy（構造類似の不活性化合物?）には一律に pK_d/i = 5.95 を設定しています。従って予測値が6よりも大きければ活性あり、となりそうです。

スコアの分布をみてみましょう。ヒストグラムを描きます。

import matplotlib.pyplot as plt
%matplotlib inline
plt.hist(Rescored_df['RFScoreVS_v2'], bins=20, range=(5.8, 6.8))

f:id:magattaca:20190407181652p:plain

要約統計量を確認します。

Rescored_df['RFScoreVS_v2'].describe()
"""
count    449.000000
mean       6.005943
std        0.075275
min        5.950000
25%        5.960443
50%        5.977587
75%        6.023838
max        6.438919
Name: RFScoreVS_v2, dtype: float64
"""

最小値は 5.95で、inactiveのスコア設定にあった値が出ています。第1四分位数 5.96、第3四分位数 6.02と、3/4近くの結合ポーズが pK_d/i=6以下となっています。活性値の判断基準（6以上）からすると、活性なしと判別されたポーズが 3/4 近いということになりそうです。（・・・この解釈で良いのだろうか？？）

AutoDock Vinaのスコアと一緒にプロットしてみます。

x = Rescored_df['RFScoreVS_v2']
y = Rescored_df['vina_affinity']

plt.scatter(x,y)
plt.title("RF-Score-VS vs AutoDock Vina")
plt.xlabel("RF-Score-VS")
plt.ylabel("AutoDock Vina")
plt.grid(True)

f:id:magattaca:20190407181810p:plain

AutoDock Vinaのスコアは小さいほど良く、RF-Score-VSは大きいほど良い結果、となりますので左上から右下に向かう対角線方向に動いて入れば、予測結果に相関があることになりそうです。

図を見るとRF-Score-VS （x = 5.5, inactive）の線上に様々なy値の点が分布ししており、今回の再スコアリングではかなり異なるスコアリング結果が得られていそうです。特に興味深い結果なのが、y軸で最小値に近い点（前回記事でスコアが良いとして選択していた点）が、再スコアではかなりスコアが悪いという結果になることです。こんなに違う結果になるとは。。。

2-3. 化合物ごとのベストポーズを取得

上のプロットは全499個の結合ポーズ全てなので、ここから各化合物ごとに最も良いスコアとなっているポーズのみを取り出したいと思います。元々の化合物のIDを「original_id」というフィールドに格納しているので、これを使ってgroupbyすれば良さそうです。

# 各化合物ごとのポーズ数を確認してみる
Rescored_df.groupby('original_id').size()
"""
original_id
PB27242793     9
PB300184784    9
PB57131569     9
PB90021090     9
・・・以下略
"""

len(Rescored_df.groupby('original_id').size())
#50

original_idでグループ化したものから、RFScoreVS_v2の値が最大となっているものを取り出して新しいデータフレームとします。*9

id_group = Rescored_df.groupby('original_id')
best_scores = Rescored_df.loc[id_group['RFScoreVS_v2'].idxmax(), :]

# 化合物数を行数で確認
best_scores.shape
# (50, 18)

再度、二つのスコアでプロットしてみます。（コード省略）

f:id:magattaca:20190407182119p:plain

だいぶんスッキリした見た目になりました。要約統計量を確認します。

best_scores['RFScoreVS_v2'].describe()
"""
count    50.000000
mean      6.079934
std       0.111637
min       5.967962
25%       6.004503
50%       6.054512
75%       6.096116
max       6.438919
Name: RFScoreVS_v2, dtype: float64
"""

今回は第１四分位数で6.00、第２四分位数で6.05と、各化合物のベストスコアを抜き出したことで安心できる分布になりました。これなら希望が持てるかもしれません・・・

プロットに合わせて構造を確認してみたいと思います。@iwatobipen先生の記事Plot Chemical space with d3js based libraryを参考にさせていただき、d3jsを使います。

import mpld3  
from mpld3 import plugins
mpld3.enable_notebook()

#描画を作成するためのモジュール
from rdkit.Chem.Draw import rdMolDraw2D
from IPython.display import SVG

#svgを格納するための空のリスト
svgs_list = []

for i in range(len(best_scores)):
    view = rdMolDraw2D.MolDraw2DSVG(300,300)
    
    ROMol_col = best_scores.columns.get_loc('ROMol')
    mol = best_scores.iloc[i, ROMol_col]
    
    prepro_mol = rdMolDraw2D.PrepareMolForDrawing(mol)
    
    view.DrawMolecule(prepro_mol)
    view.FinishDrawing()
    
    svg = view.GetDrawingText()
    svgs_list.append(svg.replace('svg:,',''))

x = best_scores['RFScoreVS_v2']
y = best_scores['vina_affinity']

fig, ax =plt.subplots()
ax.set_xlabel("RF-Score-VS")
ax.set_ylabel("AutoDock Vina")
points = ax.scatter(x, y)

tooltip = plugins.PointHTMLTooltip(points, svgs_list)
plugins.connect(fig, tooltip)

f:id:magattaca:20190407182752g:plain

gifではわかりにくいですが、スポットにカーソルを合わせると構造式が表示されます。

2-4. RF-Score VSにおけるスコアTop 10の取得

Top 10を取り出してSDFとして書き出しておきます。

#スコアでソートし、上から10個を取り出す
best_scores_10 = best_scores.sort_values('RFScoreVS_v2', ascending=False).iloc[:10]
PandasTools.WriteSDF(best_scores_10, 'RF-ScoreVS_Top10.sdf', properties=best_scores.columns)

描画します。前回と同じなのでコードは省略しますが、レジェンドのRFがRF-Score VSのスコアを、ADがAutoDock Vinaのスコアを表します。

f:id:magattaca:20190407183106p:plain

AutoDock Vinaのスコアで選択したTop 10の結果と並べてみます。

f:id:magattaca:20190407183147p:plain

両スコア関数で選択される化合物の違いを見てみてると、ぼんやりとですがRF-Score-VSではリガンドの中心付近にヘテロ芳香族環が位置していたり、AutoDock Vinaによる選択と比較して、分子全体が丸まっているように見えるものなどもあります。

3. アッセイ候補提出ファイルの作成

新しいスコア関数 RD-Score VSによる再スコアリングと結果の取り出しが完了しました。AutoDock Vinaのスコアリングとそれぞれ違いがあって面白いので、各々Top 5を取り出して合わせて10化合物とし、アッセイ候補としたいと思います。（ダイバーシティー！）

提出フォーマットはタブ区切り、ということなのでcsvモジュールを使ってファイルを作ります。まずはリストを作成。

top_10_list = []

# RFScoreVSを基準に5個取り出し
RF_sup = Chem.SDMolSupplier('./RF-ScoreVS_Top10.sdf')

for i in range(5):
    temp_list =[]
    mol = RF_sup[i]

    rank_num = i + 1
    compound_id = mol.GetProp('original_id')
    score = mol.GetProp('RFScoreVS_v2')
    
    temp_list = [rank_num, compound_id,score]
    top_10_list.append(temp_list)

# AutoDock Vinaを基準に5個取り出し
ADVina_sup = Chem.SDMolSupplier('./SBVS_Top10.sdf')

for i in range(5):
    temp_list =[]
    mol = ADVina_sup[i]

    rank_num = i + 6
    compound_id = mol.GetProp('original_id')
    score = mol.GetProp('vina_affinity')
    
    temp_list = [rank_num, compound_id,score]
    top_10_list.append(temp_list)

10化合物についてIDとスコアの組み合わせのリストのリストが作成できたので、タブ区切りのファイルを書き出します。

import csv
with open('./top10_magattaca.txt', 'w') as f:
    writer = csv.writer(f, delimiter='\t')
    writer.writerows(top_10_list)

念のため読み込んで確認

with open('./top10_magattaca.txt') as f:
    print(f.read())
"""
1 Z393761442  6.438919
2 Z243276918  6.404471
3 Z200290200  6.386205
4 Z368242106  6.385228
5 Z278994950  6.207442
6 Z237879862  -10.9
7 Z19455181   -10.9
8 Z872606870  -10.8
9 Z2788907736 -10.4
10    PB90021090  -10.4
"""

無事目的のファイルが作成できていそうです。果たしてこの中に活性のあるものはあるでしょうか？？？

同様に、500化合物のリストも同じファイル形式としておきます。スコアとしてはLBVSを行なったときの類似性評価のTanimoto係数によるスコアを記載しておきます。　（色々こねくり回して494個と500個に少し足りませんがご容赦いただけると信じて）

Top500_sup = Chem.SDMolSupplier('./SCR_compounds_SimScore.sdf')
top_500_list = []
rank_num = 0

for mol in Top500_sup:
    rank_num +=1
    compound_id = mol.GetProp('original_id')
    score = mol.GetProp('Similarity_Score')
    temp_list = [rank_num, compound_id,score]
    top_500_list.append(temp_list)

完了！　これを提出フォームに登録すればOK！

4. まとめ

今回はODDTに関して調べているうちに行き当たった論文、Sci.Rep. 2017(7)46710 に魅かれて機械学習を用いたドッキングスコア関数について調べてみました。実際に適用して見るとAutoDock Vinaによるスコアリングと異なる傾向が見え、とても面白い結果だと感じました。文献中ではまだまだ修正・改善が必要だ、と書かれていましたが、既存のスコア関数のモデルの問題点など、非常に学ぶことが多い論文だと思いました。

前回の記事でドッキングを行なった際には、とりあえず動いてデータが出れば良し！ということで進めていきましたが、AutoDock Vinaが少し古いソフトウェアなのかな？これだけで良いのかな？という懸念がありました。ドッキングそのものは時間がかかるので、何度も条件を変えて試すのは大変ですが、結合ポーズが出ていればスコア関数を修正して評価精度をあげることができる、というのを学べて良かったです。（新しいっぽいことをやって格好つけたかった）

とりあえずデフォルト+α のことが少しできたので満足！

5. 総括

提出も完了したので、最後にこれまで行なってきた処理手順をまとめます。

1. 前処理（脱塩、構造標準化）
1. 指標による絞り込み（主にRule of 5の適用）
1. 部分構造による絞り込み（オルト置換ビフェニル構造）
1. ファーマコフォアスクリーニング（ファーマコフォアを変えて２回実施） ➡︎ 化合物 Top 500 (厳密には494個)
1. LBVS（既知活性化合物に対するフィンガープリントベースの類似性評価による順位づけ） ➡︎ 化合物 Top 50を次工程へ
1. SBVS (AutoDock Vina + RF-Score VS) 　　　　 ➡︎ 化合物 Top 10をアッセイ希望に選択

以上の６ステップです。書いてしまうとこれだけ、という感じもありますが、各工程でびっくりするくらいに行き詰まったので予想以上に時間がかかってしまいました。（当初、３月上旬提出だったのに、ずるずる遅れてしまいました。遅れている上に余計な寄り道ばっかりしてしまいましたし・・・。締め切りを４月中に延ばしてくださった創薬ちゃんに感謝です。）

展望も戦略もないままはじめて右往左往していましたが、行き詰まるたびにTwitterで解決方法や、方針、文献諸々を一からご教示いただき、なんともありがたい気持ちでいっぱいです。この場を借りて皆様にお礼申し上げます。今後ともよろしくお願い申し上げます。（図々しくてすみません）

あとは、選んだ化合物の中に活性を持つものがあることを祈って！果報は寝て待ちましょう。
今回の記事も間違いがたくさんあると思います。ご指摘いただければ幸いです。

*1:タンパク質計算科学基礎と創薬への応用　神谷成敏・肥後順一・福西快文・中村春樹　著　共立出版 2009年初版１刷

*2:Uehara, S., Tanaka, S., CICSJ Bulletin 2016(34)10

*3:Machine-learning scoring functions to improve structure-based binding affinity prediction and virtual screening Ain, Q.U., Aleksandrova, A., Roessler, F.D., and Ballester P.J. WIREs Comput. Mol. Sci.2015(5)405

*4:A machine learning approach to predicting protein-ligand binding affinity with applications to molecular docking. Ballester, P.J. and Mitchell., J.B.O., Bioinformatics. 2010(26)1169

*5:Does a more precise chemical description of protein-ligand complexes lead to more accurate prediction of binding affinity? Ballester, P.J., Schreyer, A., Blundell. T.L., J. Chem. Inf. Model. 2014(54)944)、RF-Score v3 ((Li, H., Leung, KS., Wong, M.H., Ballester P.J., Mol. Inform. 2015(34)115

*6:ML to improve scoring functions binding and virtual screening
RF-Score: a Machine Learning Scoring Function for Protein-Ligand Binding Affinities

*7:Mysinger,M.M.,Carchia,M.,Irwin,J.J.,Shoichet,B.K.,J.Med.Chem.2012(55)6582）

*8:「タンパク質計算科学」の記述参考

*9:pandasのgroupbyを使ってグループ内で最大値（最小値）を持つ行を取得する