SBVSしようとする話 - magattacaのブログ

これまで以下の手順で化合物の絞り込みを行ってきました。

記述子といった指標等を用いた絞り込み
活性化合物に頻出の部分構造を用いた絞り込み
共結晶構造に基づくファーマコフォアによるスクリーニング
フィンガープリントを用いた活性化合物群への類似性に基づくスコア化

創薬レイドバトル2018の提出リストは

候補化合物 ID リスト (1): 500 化合物
候補化合物 ID リスト (2): (1) から更に 10 化合物に絞った優先順位付きリスト

です。すでに上記絞り込みにおいて、手順3.で残り約500個、手順4.でスコアによる順位づけも完了したため、リストはできた！！！・・・と言いたいところですが、スコア上位10化合物にいまいち自信が持てません。

残る課題としては、タンパク質側の情報を明示的には取り扱っていない、ということでしょうか？　手順3、においてファーマコフォア作成時に共結晶構造を参考にはしたものの、情報の一部を抽出したもので間接的な取り扱いです。やはり最後は直接対決！！
・・・ということで、ドッキング（SBVS?）に挑戦しようと思います。
絞り込んだスクリーニング化合物群の中に、無事 PD-L1 のお眼鏡に叶う候補が含まれていることを祈ってシミュレーションしていきましょう。

1. AutoDock
2. ドッキング？
3. ドッキング実施
4. SBVSにむけたスクリーニング化合物の準備
5. AutoDock vinaとODDTでSBVS
6. まとめ

1. AutoDock

早速RDKitをimportして、、、、といきたいところですが、RDKitにはドッキングの機能はないようなので別の方法を探さないといけません。今回は日本語で解説してくださっている記事、資料が検索したら出てきたという安直な理由でAutoDock vinaを使いたいと思います。*1

AutoDockはThe Scripps Research InstituteのThe Olson Laboratoryにより開発されているLigand-Protein Dcokingソフトウェアだそうです。説明によると、結合自由エネルギーの予測と遺伝的アルゴリズムに基づく探索により、タンパク質とフレキシブルなリガンド間でドッキングを実施します。今回使用しようとしてるAutoDock vinaはAutoDockを高速化、精度を高めた改良版です*2。また、ドッキングの実施や結果解析を手助けするGUIツールとしてAutoDockToolsというソフトウェア開発も行われています（こちらは他のツールとまとめてMGL Toolsに含まれています）。

これらを使えばいい感じにタンパク質にリガンドが結合するか否かわかるはず！・・・でも、そもそもドッキングって具体的に何を計算しているのか？？

2. ドッキング？

2-1. ドッキングについて

ドッキングに関して、日本化学会情報化学部会誌（CICSJ、現在はケモインフォマティクス部会？）にちょうど良い文献（タンパク質-リガンドドッキングの現状と課題 *3 ）があったので参照します。

同文献によると、ドッキング計算（protein-ligand docking）は

タンパク質の立体構造を顕に扱うSBDD(structure-based drug design)で特に重要な役割を果たす
計算機により標的タンパク質とリガンドの立体的な結合構造と結合親和性を予測する手法

だそうです。これを利用した仮想スクリーニング（virtual screening, VS）がSBVS（の一つ？）ということになります。

原子レベルで相互作用機序を理解できるため、合理的な医薬品設計の手助けとなる一方、タンパク質のような非常に多数の原子を含む構造を扱うのは計算コストが高すぎるという問題があるため、簡易化したモデルが使われます。
「簡易化」 ≒ 「タンパク質、リガンドのもつ立体構造の柔軟性をどの程度まで考慮に入れるか？」で、柔軟さを許容すればするほど計算コストが高くなります。複雑さの段階に応じて以下の４つのモデルが紹介されていました。

Rigid Docking: 鍵と鍵穴モデル(Lock-and-Key Model)に従い、タンパク、リガンド共に剛体として扱う
Flexible Docking: リガンド側のみ柔軟性を考慮（タンパク質は剛体, Flexible Ligand Model）
Flexible Side-Chains Docking: リガンドに加えて結合ポケットのアミノ酸側鎖の自由度も考慮（Paritally Flexible Protein Model）
Flexible Protein Docking: タンパク、リガンド共に柔軟性を考慮（Induced-Fit Model）

このうち、標準的なものは「3. フレキシブルドッキング」で、AutoDockもここに含まれています。

ドッキングの目的である「結合構造、結合親和性の予測」にあたって、重要となるのが効率的な探査アルゴリズムと高精度なスコア関数の２つです。前者は安定な結合構造を予想、探索するもので、これにより得られた結合構造がスコア関数により順位づけされます。文献中ではどちらも様々な手法の例が挙げられていましたが、AutoDockの場合、探査アルゴリズムとして遺伝的アルゴリズム、スコア関数として力場に基づいたスコア関数が使われています。

2-2. 遺伝的アルゴリズム？

AutoDockの遺伝的アルゴリズムについて情報処理学会(IPSJ) 第75回全国大会論文集（2013年 p.901-902 電子図書館リンク）に記載があったので引用します。*4

AutoDockでは、最適解を探索するアルゴリズムとして遺伝的アルゴリズム (GA: Genetic Algorithm) を用いている。遺伝的アルゴリズムとはデータを遺伝子として表現し、複数の個体を適応度によって選択し、交叉、突然変異などの遺伝的操作を行いながら最適解を探索するアルゴリズムである。

また、愛媛大学　村上研究室のページ（遺伝的アルゴリズム）によると、特徴として

解空間構造が不明
決定的な優れた解法が発見されておらず、
全探索が不可能と考えられるほど広大な解空間を持つ問題に有効

という点があげられるそうです。確かにタンパク質と低分子の結合構造という、無数の組み合わせがありそうな問題を解く上で有効そうです。

先のIPSJ論文集に記載されているGAの流れを少し書き換えると以下のような雰囲気です。

ランダムで初期個体として複数の結合構造を発生（最初の世代）（原子座標：遺伝子、リガンド:個体、染色体 ）
スコア関数に基づき各結合構造のエネルギースコアを計算（適応度の評価）
スコアの良い構造を適応しているとして選択(selection)
選択した構造の間で部分構造を入れ替える（交差(crossover)）
さらに一定の確率で部分構造を変化させる（突然変異(mutation)）（局所最適(local minimum)に収束するのを防いで、大局的最適(global minimum)に近づけるため？）
5.までで残った構造を新しい世代として、再度2.→5.のステップを行う
最大世代数 or 最大評価回数に達したところで反復計算を終了
最終世代の結合構造の中で、最安定エネルギーのものを最適解として出力

生命科学の用語をアルゴリズムに用いるというのは面白いですね。これがミームというやつか・・・（ちがうか）

2-3. スコア関数

探索方法の流れが把握できたので、エネルギーを計算するスコア関数について眺めてみます。先ほど省略しまししたが、CICSJの文献ではスコア関数として３種類が挙げられていました。（AutoDockのスコア関数は1に分類）

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

1.はファンデルワールス(vdW)相互作用や静電相互作用といった物理化学的な相互作用の計算を積み上げて構成されており、比較的構成が簡単になる一方で、あまりにもモデル化されすぎているため精度があまり高くないという問題があるそうです。そこで、1.の関数に類似したものをもちいつつ、現実の実験データ（既知のタンパク質-リガンド結合自由エネルギー）を取り込んでパラメータをフィッティングさせた経験的なスコア関数 2.があり、こちらの方がΔGの再現性が上がる傾向にあるそうです。1.、2.に対して、3.はさらに実験データに重みを置いており、既知の結晶構造を統計的に処理することで得た確率分布からスコア関数が作成されているそうです。

雰囲気的には「1. → 2. → 3.」の順に「理論 → 実験」の比重が変化している感じです。*5

では、AutoDockのスコア関数ではどのような物理化学的相互作用が考慮されているのか？
以下のような式となっているそうです。

ΔG = ΔE_vdw + ΔE_{electrostatic} + ΔE_H-bond + ΔE_desolvation + ΔS_conf

各項は以下に対応しているとのことです。

ΔE_vdw : ファンデルワールス相互作用（Lennard-Jones (12,6) ポテンシャル）
ΔE_{electrostatic}: 静電相互作用（距離依存性誘電率 *6 ）
ΔE_H-bond : 水素結合による相互作用 (Lennard-Jones (12,10) ポテンシャル/ Goodford directionality)
ΔE_desolvation : 脱溶媒効果
ΔS_conf : リガンドのエントロピー変化の項

上記のうち、最初の３つの項は物理化学的項目としてよく見かけるものですが、最後の２つドッキングならではという気もします。
リガンドがタンパク質に結合するにあたって、溶媒（水）に囲まれた環境から、溶媒が外れ（脱溶媒）タンパク質の結合サイト（より疎水性の高い環境）に移行するという過程をたどると考えます。この時のエネルギー変化を考慮するためのものがΔE_desolvationです。また、溶液中では自由な立体構造（配座、コンフォメーション）をとっているリガンドですが、結合に際してがポケットに合うように固定化されこの自由度が失われます。この分のエントロピー変化を考慮するものがΔS_conf となっているようです。

このあたり字面を追ってもよくわからなかったのですが、ワシントン大学 Jay Ponder Labの講義資料Tutorial & Introduction to AutoDockでなんとなくイメージができました。*7

AutoDock vinaはAutoDockの改良版、とのことなのでスコア関数も変わっているかもしれませんが基本的には同じだと思います。（調べてなくてすみません・・・）

ぼんやりとドッキングの仕組みがわかってきたので実際の作業に移りましょう。

3. ドッキング実施

3-1. AutoDock vina の処理流れ

計算化学.com ドッキングシミュレーションのやり方(AutoDock vina）　にはAutoDock Vinaに必要なファイルとして「タンパクの可動部位」「タンパクの非可動部位」「リガンド」の３つの「PDBQTファイル」と、「inputファイル」の４つが必要と記載されていました。

PDBQTファイルはpdbファイルと異なり、原子の座標に加えて部分電荷(parital charge)とAutoDockのアトムタイプを含むファイル形式となっているそうです*8。より具体的にはPDBファイルに水素原子の負荷、部分電荷（ex. Gastiger部分電荷）の情報を付加してやれば良い、ということのようです。

上記に従ってタンパク質とリガンドのファイルを用意し、以下の流れでドッキングを行います。

step 1. Grid Box（ドッキングシミュレーションを行う探索範囲）の作成タンパク質の結合ポケットを範囲として指定します
（Boxの中心座標 + Boxのサイズ情報）
step 2.ドッキングを実施
step 3.結合親和性（エネルギー）、結合構造を得る

では、具体的にやっていきましょう.

3-2. 使用するタンパク質の選択

実際にどのデータを使ってドッキングを行うか？
今回は共結晶構造があるのでこちらのタンパク質情報を使います。複数報告されていますが、このうち「PDBid: 5NIX」を用いたいと思います。理由は以前の記事で比較した際に、この複合体はリガンドが大きくため、タンパク質側の結合ポケットが広がっている可能性が高い、と考えられるからです。

なぜポケットが大きい方が良いか？「タンパク質計算科学基礎と創薬への応用」の記述を引用します。(p. 217)*9

・・・注意すべきは、タンパク質-化合物ドッキングの精度は通常２Å以上の粗い精度であって、タンパク質と化合物の原子はたいていの場合、vdW衝突しているということである。・・・（中略）・・・したがって、vdWポテンシャルの衝突を減らすことが、ドッキングソフトの重要な技術となっている。

なるほど！それならば最初から空間が大きいもの選んどいた方が良いよね、、、っていう浅い考えです。（この本、積読してないで早く読めばよかった・・・）

3-3. Pymolでタンパク質を前処理

ドッキングに向けて「PDBid: 5NIX」の構造をPymolでちょっと処理しておきます。

まず、この構造は４量体（Chain A、B、C、D）となっているので、簡便のためChain C、Dを削除しておきます。

PyMOL> remove (chain C,D)
 Remove: eliminated 2093 atoms in model "5nix".

ついでに水も除きます。

PyMOL>remove resn hoh
 Remove: eliminated 133 atoms in model "5nix".

remove solventでも良いらしいです。参考

Sequenceを表示させて確認すると、EDO（エチレングリコール）というものもあったので、不要そうなので除去しておきます。

PyMOL>remove resn EDO
 Remove: eliminated 4 atoms in model "5nix".

以上の操作でだいぶスッキリした見た目になりました。

f:id:magattaca:20190403224018p:plain

3-4. Grid Box設定のための座標確認

ドッキングの中心（作成するGridの中心）を決めるために、座標を取得します。今回はファーマコフォアスクリーニングを実施済みの化合物群をドッキングに使用予定のため、ファーマコフォア作成時に利用したリガンド原子部位の中心をGrid中心とします。

PyMOL>pseudoatom center, sele
 Warning: 'center' is a reserved keyword, appending underscore
 ObjMol: created center_/PSDO/P/PSD`1 /PS1

f:id:magattaca:20190402232114p:plain

安易にcenterという名前をつけたら怒られました・・・。このpeudoatomの座標は以下です。

PyMOL>xyz = cmd.get_coords('center_', 1)
PyMOL>print xyz
[[ -5.486217  11.170054 -24.066675]]

全体を5nix_model.pdb、タンパク質のみ（Chain A,B）を5nix_chainAB.pdb、リガンドのみをmol2形式5nix_lig.mol2としてそれぞれ出力しておきました。

3-5. テストドッキング

PDBの共結晶構造から、余分なものを除いたタンパク質のみのPDBファイルと、リガンドのみのmol2ファイルが用意できたのでこれらを用いてドッキングの流れを確認してみます。

参考にしていたページでは準備にAutoDock Tools（MGL Tools）を使うとのことでしたが、私のMac上ではエラーが出てうまく動かなかったので、代わりにUCSF ChimeraのpluginとしてAutoDock vinaを使用します。（Chimeraはアカデミックフリーですが商用利用はライセンスが必要だそうです）*10

3-5-1. タンパク質/リガンドの前処理

AutoDock vinaに合わせて、先に用意したタンパク質PDBファイルを処理します。
ChimeraでPDBファイルを読み込んだのち、Tools の Surface/Binding Analysis から DockPrep　を選択すると新しい画面が開きます。

f:id:magattaca:20190402232521p:plain

この画面の指示に従って、

溶媒の除去 (Delete Solvent)
水素原子の付加 (Add hydrogens)
電荷の付加 (Add charges)

といった必要な処理を選択し、Mol2ファイルで書き出します。

f:id:magattaca:20190402232635p:plain

リガンドのファイルも同様に水素原子と電荷の付加の処理をしました。カルボキシル基には水素原子が付加されておらず、アミノ基には付加されているので、それぞれ負電荷、正電荷を帯びるように処理されているようです。（下図、黄色丸部位）

f:id:magattaca:20190402234105p:plain

3-5-2. ドッキングのジョブを投げる準備

一度「close session」した後、処理したタンパク質とリガンドのmol2ファイルを読み込みます。
Tools の Surface/Binding Analysis から AutoDock Vinaを選択し、開くサブ画面に設定を書き込んでいきます。

f:id:magattaca:20190402234444p:plain — 設定の画面例

上の図で行なっている設定は、

output fileの名前（今回は5nix_redockとした）をつける
タンパク質（receptor）とリガンドを指定
(mol2ファイルを読み込んでいるのに名前は処理前のもので表示された・・・？）
Gridの指定
AutoDock vinaを行う場所の設定（web？ or local？）

Gridの設定は、「Receptor search volume options」の部分で行なっており、Grid中心（Center）を先にpymolで確認した座標を参考に設定、Grid Boxの大きさを size x:20 y:20 z:20としてます。リガンドがBoxからはみ出していたのでGrid中心を少しX軸方向にcenterの座標の位置をずらし、微調整をしてあります（x:-7, y:11, z:-24）。

AutoDock vinaを実行する場所はpublic web serviceか、あるいはlocalを選択できます。今回はlocalで実施するためvinaを入れたpathを指定します。

以上で準備完了、あとは実行のみ

3-5-3. ドッキングの実行と結果

実行・・・Macが発熱します。

結果・・・４つの構造が残ったようです。

mode	affinity (kcal/mol)	RMSD l.b.	RMSD u.b.
1	-11.5	0.000	0.000
2	-11.5	1.155	1.585
3	-11.1	1.280	1.927
4	-10.8	3.173	6.337

pdbqtファイルがいくつか出力されているのでpymolで構造を眺めてみます。 pymolで出力5nix_redock.pdbqtを読み込むと4つの構造をstateとして含む状態で読み込まれます。少しみづらいので「Action → state → split」とすることで一つずつに分けます。

共結晶構造中のリガンド（赤色）とともにそれぞれ描画してみます。

f:id:magattaca:20190402235052p:plain

エネルギー的に近いmode 1、2、3に対してmode　4ではかなり結合ポーズが変化しています。ちょうどファーマコフォアには含めなかったベンジル側鎖がフリップしたような構造となっています。より安定なmode 1、2 はほぼ元々のリガンドに重なっており、今回の結果ではかなり再現性の良い結果となった、と言えそうです。（初期構造もグリッドも共結晶構造をもとに設定しているので当然？かもしれませんが。。。）

スコアに関してですが、香川大学のページの記載（Autodock4の場合）によると、

クラスターの平均エネルギーの差が約2.5 kcal/mol以下の場合，AutoDock force fieldの標準偏差以内であり，エネルギーからどちらが「正しい」かをいうことは出来ない

とのことだそうです。今回この標準偏差よりも安定化しているので、良い結合ポーズ・・・ということで良いのでしょうか？再度「タンパク質計算科学基礎と創薬への応用」の記述を引用します。(p. 225)

・・・通常、ΔGの再現は悪く、2.5~3kcal/mol程度の精度が出ればいいほうであるが、ヒット化合物のΔGは4~8kcal/mol程度、最適化された医薬品で10kcal/mol程度なので、ΔGの計算のみでヒット・リード化合物の活性地を評価するのは困難である。・・・

初版が10年ほど前の書籍なので、現在のソフトウェアの精度はこの時よりも上がっていると思われますが、AutoDock vinaの報告が J.Comp.Chem. 2010(31)455、なので今回の目安としては有効なのではないかと思います。

4. SBVSにむけたスクリーニング化合物の準備

テスト構造での検証、手順の流れがわかってきましたのでようやく本題のスクリーニングに移りたいと思います。スクリーニング対象化合物としては、とりあえずLBVSのスコアTop 50としたいと思います。（私のPCの性能的にこれくらいが限界、、、）

4-1. ファーマコフォアを利用したコンフォマーの発生

まずはリガンドの立体構造が必要なので初期構造を準備します。折角なのでRDKitを利用して、これまでのスクリーニングの結果を応用したいと思います。

手元にあるスクリーニング対象化合物は、活性化合物群に対する類似性をベースにスコアリングを実施したものです。したがって、共結晶構造中のリガンドとある程度の類似性があると思われます。そこで、安定かつ共結晶構造中のリガンドにできるだけ近い配座を用意したいと思います。

具体的には

step 1. ファーマコフォアをembed
step 2. ETKDGでembedした構造を最適化
step 3. 一つの化合物に複数embedできた場合、Feature-mapで最も共結晶のリガンドに近い構造を代表として選出する。

という流れです。

まずはLBVSのスコアを付与したSDFを読み込み、スコアでソートしてTop50化合物を取り出しました。そしてファーマコフォアの設定には、以前ファーマコフォアサーチで用いたものと同じ値を使用し、設定を行いました。この辺りの作業は繰り返しなので省略します。

前回との違いは、このファーマコフォアサーチの関数においてさらに立体構造の最適化のステップを加えることです。以前ファーマコフォアサーチに使った関数はembedできた数のみを返すようにしていましたが、これをconfomerのリストを返すように書き換えます。また、この関数では不斉の情報を取り除いていますが、幸い（？）今回の化合物リストに不斉の情報はなさそうなのでそのまま進めます。

def PConfs(mol):
    #立体の情報を除去
    smi = Chem.MolToSmiles(mol, isomericSmiles=False)
    mol_re = Chem.MolFromSmiles(smi)
    AllChem.Compute2DCoords(mol_re)
    
    #フィーチャーをそもそも持っているか？
    match, mList = EmbedLib.MatchPharmacophoreToMol(mol_re, Ffactory, pcophore)
    
    if match == True:
        #距離の検証
        #距離行列の取得
        bounds = rdDistGeom.GetMoleculeBoundsMatrix(mol_re)

        #ファーマコフォアのマッチング 
        pList = EmbedLib.GetAllPharmacophoreMatches(mList,bounds,pcophore) 
        
        #pListのそれぞれについてFeatureとマッチする原子のidを取得する
        num_match = len(pList)
        phMatches = []
        for i in range(num_match): 
            num_feature = len(pList[i])
            
            phMatch = []
            
            for j in range(num_feature):
                phMatch.append(pList[i][j].GetAtomIds())
                
            phMatches.append(phMatch)
        
        #原子のidとファーマコフォアをもとに3次元構造を埋め込む    
        confs = []

        for phMatch in phMatches:
            bm,embeds,nFail =EmbedLib.EmbedPharmacophore(mol_re, phMatch, pcophore,count=5, silent=1)
            for embed in embeds:
                #水素原子を付加
                embedH = AllChem.AddHs(embed)
                #構造最適化
                AllChem.EmbedMolecule(embedH, AllChem.ETKDG())
                
                #元々の分子のプロパティを追加する作業
                PropNames = list(mol.GetPropNames())
                prop_dict = {}
                #プロパティの値を取り出し
                for n in PropNames:
                    vl = mol.GetProp(n)
                    prop_dict[n] = vl
                #confomerに付与
                for k, v in prop_dict.items():
                    embedH.SetProp(k, v)
                
                #リストに追加
                confs.append(embedH)         
        return confs            
    else:
        return 0

スクリーニングに用いる50個の化合物の一番最初の分子を取り出し、関数の挙動を確認します。

test_m = top50_mols[0]
test_confs = PConfs(test_m)
len(test_confs)
# 4

コンフォマーが４つ生成されました。py3Dmolを使った３Dの描画結果を以下に貼っておきます。

f:id:magattaca:20190403000118p:plain

4-2. Feature Mapによるコンフォマーの選別

次に、得られたconfomerの中から共結晶中のリガンドコンフォマーに最も類似したものをFeature Mapを使って選択します。RDKitのブログ記事Using Feature Mapsを参考にしました。元の文献を読めていないので理解できていないのですが、複数のコンフォマーの中から良いものを選ぶことができそうな感じです。また、比較対象のリガンドとして、ファーマコフォアを満たすのに必要最小限の構造をもつように見える、「PDBid:5N2F」のリガンド（8HW）を選択しました。

RDKitのブログ記事をそのまま関数にしていきます。

from rdkit.Chem.FeatMaps import FeatMaps
from rdkit.Chem import rdMolAlign

def FeatMapScorer(mol, tmp):
    #比較対象(tmp)との間でアラインメントをとる
    o3d = rdMolAlign.GetO3A(mol,tmp)
    o3d.Align()
    
    #feature mapのパラメーターの辞書を作成（デフォルトの値を使用）
    fmParams = {}
    for k in Ffactory.GetFeatureFamilies():
        fparams = FeatMaps.FeatMapParams()
        fmParams[k] = fparams
    
    #featureのうちより一般的なものだけを残すための設定
    keep = ('Donor','Acceptor','NegIonizable','PosIonizable','Aromatic')
    featLists = []
    for m in [mol, tmp]:
        #featureをまず全部取得
        rawFeats = fdef.GetFeaturesForMol(m)
        #keepにある欲しいものだけリストに追加 
        featLists.append([f for f in rawFeats if f.GetFamily() in keep])
    
    #FeatMap オブジェクトを作成
    fms = [FeatMaps.FeatMap(feats=x, weights=[1]*len(x), params=fmParams) for x in featLists]
    
    #feature mapに対するスコアを取得
    score = fms[0].ScoreFeats(featLists[1]) / min(fms[0].GetNumFeatures(), len(featLists[1]))
    
    #scoreを返す
    return score

先のテスト化合物で得た４つのコンフォマーについて、Feat Mapのスコアを求めてみます。

test_scores = []

for conf in test_confs:
    score = FeatMapScorer(conf, temp_8HW)
    test_scores.append(score)

print(test_scores)
#[0.25512226843031627, 0.17970526398716216, 0.05047724449395227, 0.025276456816820638]

スコアが最大となるコンフォマーのみを取り出したいので、numpyのargmaxを利用します。

import numpy as np

test_score_array = np.array(test_scores)
test_best_conf_index = np.argmax(test_score_array)
test_best_conf = test_confs[test_best_conf_index]

上記test_best_confを描画し、無事立体構造が選別できていることが確認できたので、本番、50個の化合物群全てに適用します。

Top50_3D_list = []

for mol in top50_mols:
    confs = PConfs(mol)
    
    FMapScores = []
    
    for conf in confs:
        score = FeatMapScorer(conf, temp_8HW)
        FMapScores.append(score)
    
    ScoreArray = np.array(FMapScores)
    BestConfIndex = np.argmax(ScoreArray)
    BestConf = confs[BestConfIndex]
    
    Top50_3D_list.append(BestConf)

Top50_3D.sdfとして書き出しておきました。

4-3. PyMolでコンフォマーを眺める

本当にテンプレート構造に対して重なりの良いコンフォマーが選ばれているのか、PyMolで確認してみます。テンプレートにした元のリガンド（8HW）と、生成した50化合物のコンフォマーを含むSDFをPyMolで読み込みます。複数の化合物を含むSDFは、1つの化合物を1つのstateとするファイルとして読み込まれているようです。右下再生ボタンを押すと、順番に各stateを表示するムービとして眺めることができました。

f:id:magattaca:20190403230400g:plain

いくつかははみ出しているように見えますが、全体的に良い重なりとなっているように見えます。

試しにタンパク質も同時に表示させて、違う角度から眺めてみます。

f:id:magattaca:20190403230828g:plain

アミノ酸残基の側鎖を表示させていない状態でも、化合物がタンパク質にぶつかっている（重なっている）様子が見て取れます。やはりリガンドのコンフォマーに似せるだけでは、最良なポーズ、とは行かないようです。AutoDock Vinaによるドッキングは、この化合物たちの中から上手く結合ポケットに収まるようなポーズを見つけることができるでしょうか？

5. AutoDock vinaとODDTでSBVS

スクリーニング対象化合物群の３次元構造の生成が完了したので、これを用いてPD-L1に対するSBVSを実施したいと思います。
・・・しかし一つ一つChimeraで実行するのは大変。。。コマンドもイマイチ理解できないし・・・
と思っていたところ、またしても便利なソフトウェアがありました。その名もOpen Drug Dicovery Toolkit!!（略称 ODDT）。さまざまなオープンソースのケモインフォマティクスのツールがあって便利になったものの、それぞれのinput、outputが異なっていて面倒、、、という問題を解決するため、一つのパッケージに統合してもっと便利にしていきましょう！という趣旨のようです。

元の文献はこちら *11。ドキュメントも充実しています（リンク）。@iwatobipen先生のブログ記事でも使用方法が紹介されていました（Open drug discovery toolkit for python）。

早速使っていきましょう。

5-1. ODDTでchimeraでの操作を追試

いきなりスクリーニング本番を行なって失敗するとショックなので、まずは先にUCSF chimeraで行なったテストドッキングと同じことができるか、練習してみましょう。

まずはODDTをimport...

import oddt
from oddt import toolkit

ドッキングに使用する相手のタンパク質のpdbファイル（先にchimeraで読み込んだものと同じもの）を読み込みます。

protein = next(oddt.toolkit.readfile('pdb', './5nix_docking/5nix_chainAB.pdb'))
protein.protein = True

うまく読み込め、タンパク質であると認識してくれているようです。使い方はヘルプでも参照可能です。

help(oddt.toolkit.readfile)
"""
 Help on function readfile in module oddt.toolkits.ob:
readfile(format, filename, opt=None, lazy=False)
"""

次にリガンド（テストドッキングでは共結晶構造のリガンド mol2ファイル）を読み込みます。

ligand = next(oddt.toolkit.readfile('mol2', './5nix_docking/5nix_lig.mol2'))

Autodock vinaを使うためモジュールを読み込みます。

from oddt.docking import AutodockVina

AutodockVina.autodock_vina()で、ドッキングパラメータを指定します。引数はそれぞれ、

protein : ドッキングに使うタンパク質
size : Grid Boxのサイズ、前回同様に (x, y, z) = (20, 20, 20)
center : Grid Boxの中心、前回と同じ座標 (x, y, z) = (-7, 11, -24)
exhaustiveness、num_modes、energy_range: すべてchimeraのAutodockVina pluginと同様
executable: AutodockVinaのlocalのpathを指定　（相対pathの深さから、私のファイル管理がグチャグチャなことがバレますね・・・）

という感じです。

test_docking = AutodockVina.autodock_vina(protein=protein, size=(20, 20,20), center=(-7, 11, -24),
                                         exhaustiveness=8, num_modes=9, energy_range=3, n_cpu=3,
                                         executable='../../../../autodock_vina_1_1_2_mac/bin/vina')

パラメータを設定してインスタンス（？）ができたら、dock()でドッキングを実行です。引数にligandを指定してやります。

dockin_reslut = test_docking.dock(ligand)

4,5分で計算終了しました（chimeraの時と同じくらい）。水素の付加や部分電荷の付加といった作業を行なっていませんが、勝手にそのあたりの前処理もしてくれるみたいです。

len(dockin_reslut)
# 4

生成された結合ポーズは４つ。リガンドの構造のリストとなっており、それぞれの構造に結合エネルギーなどの情報が付与されているようです。ODDTのドキュメントに従って情報を抜き出してみます。

type(dockin_reslut)
# list
dockin_reslut[0].data.keys()
"""
['vina_affinity',
 'vina_rmsd_lb',
 'vina_rmsd_ub',
 'vina_rmsd_input',
 'vina_rmsd_input_min']
"""

親和性やrmsdといった情報があるようです。一つエネルギーを確認してみます。

dockin_reslut[0].data['vina_affinity']
# '-11.4'

Chimeraで計算した時の最小値（最安定）と同じ値となっていそうです。
4つの結果をPandasToolsのデータフレームで読み込むと以下のようになります。

f:id:magattaca:20190403001907p:plain

残りの３つのエネルギーはChimeraの時と異なるようですが、おそらく初期構造をランダムで発生させているため、初期構造に依存して異なる結果が出ているのだと思います。

ODDTは使用方法がpybel(OpenBabel?)に近い書き方となっているそうです（使ったことがないのでわかりません・・・）。ドッキングの結果で得られるリガンドのタイプを確認してみます。

type(dockin_reslut[0])
#oddt.toolkits.ob.Molecule

このobはOpenBanelの略のようです。

SDFに書き出してみます。RDKitと少し異なり、Outputfile()を用意し、構造を一つずつforループで回してwrite()していきます。

TestSDFile = oddt.toolkit.Outputfile('sdf', 'ODDT_docktest.sdf')
for mol in dockin_reslut:
    TestSDFile.write(mol)
TestSDFile.close()

一連の流れが確認でき、妥当な結果（近い結果）が得られていそうなので、実際のスクリーニングに進みます。

5-2. ODDTでスクリーニング本番

以下を実施、、、

protein_5nix = next(oddt.toolkit.readfile('pdb', './5nix_docking/5nix_chainAB.pdb'))

DockingParams = AutodockVina.autodock_vina(protein=protein_5nix, 
                                           size=(20, 20,20), center=(-7, 11, -24),
                                           exhaustiveness=8, num_modes=9, energy_range=3,
                                           n_cpu=4,
                                           executable='../../../../autodock_vina_1_1_2_mac/bin/vina')

Dock_resluts = []

for mol in oddt.toolkit.readfile('sdf', './Top50_3D.sdf'):
    dock_res = DockingParams.dock(mol)
    Dock_resluts.append(dock_res)

3時間ほどMacが騒音を出した後、停止しました・・・

len(Dock_resluts)
# 50

50個全て終わってました。

5-3. 全ドッキングポーズを出力・確認

全50個の化合物について、各化合物で生成されたドッキングポーズを全て格納したリストを作成します。

All_Result = []

for poses in Dock_resluts:
    for pose in poses:
        All_Result.append(pose)

len(All_Result)
# 449

449個のポーズが得られました。プロパティが保持されているかどうか確認してみます。

# keyを確認
Dock_resluts[0][0].data.keys()
"""
['ID',  'MW',  'MolLogP',  'NumHAcceptors',  'NumHDonors',  'NumRotatableBonds',  'PScore',  'PScore2', 'Similarity_Score', 'TPSA', 'original_id',  'vina_affinity', 'vina_rmsd_lb', 'vina_rmsd_ub', 'vina_rmsd_input', 'vina_rmsd_input_min']
"""
#valueを確認
Dock_resluts[0][0].data.values()
"""
['', '457.48600000000027', '3.5305000000000017', '5', '2', '5', '10', '4', '20.832658', '96.970000000000013', 'Z18052460', '-8.0', '0.000', '0.000', '9.41251', '9.408608']
"""

データはキチンと取れていそうなので、間違えて消してしまわないうちに書き出しておきます。一度sdfで書き出した後、RDKitで読み込めるかどうか確認してみます。

#ODDTでSDF書き出し
DockAllResultOutput = oddt.toolkit.Outputfile('sdf', 'ODDT_Dock_All_Result.sdf')
for mol in All_Result:
    DockAllResultOutput.write(mol)
DockAllResultOutput.close()

#RDKitのPandasToolsで読み込み
All_Dock_df = PandasTools.LoadSDF('./ODDT_Dock_All_Result.sdf')
"""
RDKit ERROR: [22:07:38] Explicit valence for atom # 2 C, 5, is greater than permitted
RDKit ERROR: [22:07:38] ERROR: Could not sanitize molecule ending on line 42625
RDKit ERROR: [22:07:38] ERROR: Explicit valence for atom # 2 C, 5, is greater than permitted
・・・以下エラーが続く
"""

エラーがたくさんでました。試しにMarvin viewでひらいてみると、開くことができました。エラーが出た構造を確認してみるとヘテロ環（芳香族）の扱いにおいて結合次数の割り当てに失敗し、炭素に5つ結合が伸びている、といったような不具合が生じているようです。

f:id:magattaca:20190403003131p:plain

テキストエディタでぽちぽちsdfの結合部分を修正し、再度読み込みます。

All_Dock_df = PandasTools.LoadSDF('./ODDT_Dock_All_Result_mod.sdf')
All_Dock_df.head()

f:id:magattaca:20190403003257p:plain

今度はエラーを出さずに読み込めました。

SDFの書き出し、構造の確認が完了したので結果のドッキングポーズをpymolで眺めてみます。先と同様にタンパク質ともともとの共結晶構造中のリガンド、そしてドッキング結果のSDFを読み込み、ムービーを再生します。

f:id:magattaca:20190403233016g:plain

おお！！バーチャルスクリーニング感あふれてる！！

全結合ポーズの中では、結合ポケットの端（タンパク質と溶媒の境界近く）に位置するものも多く見られ、結合ポケット奥深く（共結晶構造リガンドと重なるような位置）まで入り込んでいるものは少ないように見えます。先にみた「タンパク質計算科学」のコメント通り、結合ポケット奥深くに入るポーズは、アミノ酸残基側鎖とのぶつかり合い(vdW相互作用)で弾かれてしまい、なかなか安定なポーズを見つけづらいのでしょうか？　あるいはただ単に探索空間であるGrid Boxの、座標位置、サイズの設定が良くなかったのかもしれません。

5-4. 各化合物のベスト構造のみを出力・確認

つぎに、各化合物について、複数の結合ポーズの中で最安定のもののみを一つずつ取り出してみます。結合親和性(affinity)の値で取り出すと、値が重複した場合に取り出す構造にブレが生じるので、rmsd値が最小のものを取り出していきます。（AutoDockVinaの結果では最安定の結合ポーズと比較したrmsd値が計算されて出力されています。）

Best_Mode_Result = []

for poses in Dock_resluts:
    rmsd_tmp = []
    for pose in poses:
        rmsd = pose.data['vina_rmsd_ub']
        rmsd_tmp.append(rmsd)
    rmsd_arr = np.array(rmsd_tmp)
    best_idx = np.argmin(rmsd_arr)
    best_pose = poses[best_idx]
    Best_Mode_Result.append(best_pose)

こちらについても一度sdfで書き出した後、エラーの構造を修正して、RDKitで読み込みます。

#ODDTでSDF出力
DockBestModeResultOutput = oddt.toolkit.Outputfile('sdf', 'ODDT_Dock_BestMode_Result.sdf')
for mol in Best_Mode_Result:
    DockBestModeResultOutput.write(mol)
DockBestModeResultOutput.close()

# 構造を修正ののち、名前を変えてSDFを保存。RDKitで再度読み込み
BM_df = PandasTools.LoadSDF('./ODDT_Dock_BestMode_Result_mod.sdf')
BM_df['vina_affinity'] = BM_df['vina_affinity'].astype(float)

スコア順にソートして、ドッキングスコアの良い10個を取り出します。

#スコア順にソート
BM_df_st = BM_df.sort_values('vina_affinity', ascending=True)
#Top10を抜き出す
Dock_Top_10 = BM_df_st.iloc[:10]

これらドッキングスコア　ベスト10の化合物群についてもPyMolのムービーで眺めて見ましょう。

f:id:magattaca:20190403234232g:plain

おお！なんだか思っていたよりも良い感じ！！

スコアの良い結合ポーズでは、結合ポケットの奥深くまで化合物が入り込んでおり、反発を抑えつつ上手くタンパク質と相互作用ができていそうな雰囲気がただよっています。

意外な点としては、複数の化合物でビフェニル部分構造が、元の共結晶構造の該当の部分構造とは重ならない形のポーズが、最もスコアが良い結合構造として選ばれていることです。これらでは疎水性のビフェニルが親水性の溶媒接触面側に位置してしまっているようにも見え、ビフェニル部分構造を基準に選んだことが結合親和性に対しては逆効果に働いてしまっているようにも思えます。部分構造の結合位置を指定してドッキングする、といった細かい設定を行っていけば、より元のリガンドに近いような結合ポーズを出せるのかもしれません。

では、これらの化合物がどのような構造式のものなのか？　LBVSの時の類似性スコアとともに10個の化合物を描画してみます。

SBVS_scores = list(BM_df_st['vina_affinity'][:10])
LBVS_scores = list(BM_df_st['Similarity_Score'][:10])

score_legends = []
for s, l in zip(SBVS_scores, LBVS_scores):
    score = 'SB: ' + str(s) + '/ LB: ' + str(round(float(l), 1))
    score_legends.append(score)

Draw.MolsToGridImage(BM_df_st['ROMol'][:10], legends=score_legends)

f:id:magattaca:20190403004218p:plain

リガンドベースで類似性のスコアのみで評価した時とは異なる分子が取り出されているようです。両者を並べてみます。

f:id:magattaca:20190403005401p:plain

丸で囲った部分が今回のドッキング（SBVS）と前回のLBVSで共通して選ばれている構造でTop 10のうち、３つとなりました。類似性のみのLBVSとはまた異なった特徴を持つ化合物を選んでくることができているようです。

6. まとめ

今回は、最後の難関ドッキングによるSBVSを行って見ました。色々とソフトウェアを行ったり来たりしてかなり見通しの悪い感じの記事になってしまいましたが、なんとか最終的にそれっぽい結果を出すことができました。

LBVSで絞り込んだ50個から最終的に10個を選び、LBVSのスコアとは異なるTop 10化合物となりました。SBVSの特徴として、既知の活性化合物とは構造的に異なるものを見つけ出す可能性がある、との記載がありましたので、その意味では異なる結果が出たことでドッキングをした甲斐があった！・・・と思いたい。

（異なる構造が選びたければ先に部分構造や類似性で絞り込むのは順番的に良くないのでは？といった、つらいツッコミはご容赦いただければ幸いです。）

得られた結合ポーズは、当初想定していたのとは異なる結合様式もあった、など、まだまだ調整・工夫の余地がたくさんありそうです。ほとんどデフォルトのコピー&ペーストですし。。。とはいえ、あまりに”ヒトの目”でみて出したいポーズを狙いすぎても折角ドッキングソフトを使っている意味が薄まってしまいますので、今回はこの結果で良し！！　

ということで一通りそれっぽいことができたので満足。

付け焼き刃で調べながら行なったので、間違い等多数あると思います。ご容赦、ご指摘いただければ幸いです。

P.S. 進捗悪すぎてすでに提出期限オーバー感があるのですが・・・今からの提出でも許してもらえるでしょうか？？？

*1:

参考記事一覧:
1. 計算化学.com ドッキングシミュレーションのやり方(AutoDock vina）
2. 計算化学.com AutoDock のスコアファンクションはどれくらい正確か？CASF-2013 ベンチマーク
3. 香川大学　ケミカルバイオロジー研究室 Docking simulation
4. Slide Share AutoDock_vina_japanese_ver.3.0

*2:AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization and multithreading Trott, O., and Olson, A.J., J. Comput. Chem. 2010(31)455

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

*4:「GPU を用いたエネルギー計算の並列化によるドッキングシミュレーションの高速化」東京工業大学　佐々木崇浩宇田川拓郎関嶋政和

*5:タンパク質計算科学基礎と創薬への応用　神谷成敏・肥後順一・福西快文・中村春樹　著　共立出版 2009年初版１刷の記述も参考にしています。

*6:Mehler, E.L. and Solmajer, T. "Electrostatic effects in proteins: comparison of dielectric and charge models" Protein Engineering 1991(4)903

*7:Jay Ponder MOLECULAR MODELING (Chemistry 478) Spring Term 2014
Lecture 20 Small Molecular & Protein Docking
Lecture 21 Tutorial & Introduction to AutoDock
あるいは
Marc A. Marti-Renom / Structural Genomics Unit / Bioinformatics Department / Prince Felipe Resarch Center (CIPF)
Docking of small molecules. AutoDock.
などを参考にしました。

*8:What is the format of a PDBQT file?

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

*10:Chimeraの使用にあたって以下のファイルを参考にしました。
* 複合体構造モデリング
 東京大学大学院農学生命科学研究科アグリバイオインフォマティクス教育研究プログラム
寺田透
Autodock tutorial VINA with UCSF Chimera
José R. Valverde
CNB/CSIC

*11:Wójcikowski, M., Zielenkiewicz, P., and Siedlecki, P. Open Drug Discovery Toolkit (ODDT): a new open-source player in the drug discovery field Journal of Cheminformatics 2015(7)26