これまで以下の手順で化合物の絞り込みを行ってきました。
記述子といった指標等を用いた絞り込み
活性化合物に頻出の部分構造を用いた絞り込み
共結晶構造に基づくファーマコフォアによるスクリーニング
フィンガープリントを用いた活性化合物群への類似性に基づくスコア化
創薬レイドバトル2018 の提出リストは
候補化合物 ID リスト (1): 500 化合物
候補化合物 ID リスト (2): (1) から更に 10 化合物に絞った優先順位付きリスト
です。すでに上記絞り込みにおいて、手順3.で残り約500個、手順4.でスコアによる順位づけも完了したため、リストはできた！！！・・・と言いたいところですが、スコア上位10化合物にいまいち自信が持てません。
残る課題としては、タンパク質側の情報を明示的には取り扱っていない、ということでしょうか？ 手順3、においてファーマコフォア作成時に共結晶構造を参考にはしたものの、情報の一部を抽出したもので間接的な取り扱いです。やはり最後は直接対決！！
1. AutoDock 早速RDKitをimportして、、、、といきたいところですが、RDKitにはドッキングの機能はないようなので別の方法を探さないといけません。今回は日本語で解説してくださっている記事、資料が検索したら出てきたという安直な理由でAutoDock vinaを使いたいと思います。*1
AutoDock The Olson Laboratory により開発されているLigand-Protein Dcokingソフトウェアだそうです。説明によると、結合自由エネルギーの予測と遺伝的アルゴリズム に基づく探索により、タンパク質とフレキシブルなリガンド間でドッキングを実施します。今回使用しようとしてるAutoDock vina *2 。また、ドッキングの実施や結果解析を手助けするGUI ツールとしてAutoDockTools というソフトウェア開発も行われています（こちらは他のツールとまとめてMGL Tools
これらを使えばいい感じにタンパク質にリガンドが結合するか否かわかるはず！
・・・でも、そもそもドッキングって具体的に何を計算しているのか？？
2. ドッキング？ 2-1. ドッキングについて ドッキングに関して、日本化学会情報化学部会誌（CICSJ、現在はケモインフォマティクス 部会？）にちょうど良い文献（ タンパク質-リガンドドッキングの現状と課題 *3 ）があったので参照します。
同文献によると、ドッキング計算（protein-ligand docking）は
タンパク質の立体構造を顕に 扱うSBDD(structure-based drug design)で特に重要な役割を果たす
計算機により標的タンパク質とリガンドの立体的な結合構造 と結合親和性 を予測する手法
だそうです。これを利用した仮想スクリーニング（virtual screening, VS）がSBVS（の一つ？）ということになります。
原子レベルで相互作用機序を理解できる ため、合理的な医薬品設計の手助けとなる一方、タンパク質のような非常に多数の原子を含む構造を扱うのは計算コストが高すぎるという問題があるため、簡易化したモデルが使われます。
Rigid Docking Flexible Docking Flexible Side-Chains Docking アミノ酸 側鎖の自由度も考慮（Paritally Flexible Protein Model）Flexible Protein Docking
このうち、標準的なものは「3. フレキシブルドッキング」で、AutoDockもここに含まれています。
ドッキングの目的である「結合構造、結合親和性の予測」にあたって、重要となるのが効率的な探査アルゴリズム と高精度なスコア関数 の２つです。前者は安定な結合構造を予想、探索するもので、これにより得られた結合構造がスコア関数により順位づけされます。文献中ではどちらも様々な手法の例が挙げられていましたが、AutoDockの場合、探査アルゴリズム として遺伝的アルゴリズム 力場に基づいたスコア関数 が使われています。
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 に分類）
力場に基づいたスコア関数 経験的なスコア関数 統計に基づいたスコア関数
1. はファンデルワールス (vdW)相互作用や静電相互作用といった物理化学的な相互作用 の計算を積み上げて構成されており、比較的構成が簡単 になる一方で、あまりにもモデル化されすぎているため精度があまり高くない という問題があるそうです。そこで、1. の関数に類似したものをもちいつつ、現実の実験データ（既知のタンパク質-リガンド結合自由エネルギー）を取り込んでパラメータをフィッティングさせた経験的なスコア関数 2. があり、こちらの方がΔGの再現性 が上がる傾向にあるそうです。1.、2. に対して、3. はさらに実験データに重みを置いており、既知の結晶構造を統計的に処理 することで得た確率分布からスコア関数が作成されているそうです。
雰囲気的には「1. → 2. → 3.」の順に「理論 → 実験」の比重が変化している感じです。*5
では、AutoDockのスコア関数ではどのような物理化学的相互作用が考慮されているのか？
ΔG = ΔEvdw + ΔEelectrostatic + ΔEH-bond + ΔEdesolvation + ΔSconf
各項は以下に対応しているとのことです。
ΔEvdw ファンデルワールス 相互作用 （Lennard-Jones (12,6) ポテンシャル）ΔEelectrostatic 誘電率 *6 ）ΔEH-bond ΔEdesolvation ΔSconf エントロピー 変化の項
上記のうち、最初の３つの項は物理化学的項目としてよく見かけるものですが、最後の２つドッキングならではという気もします。ΔEdesolvation エントロピー 変化を考慮するものがΔSconf
このあたり字面を追ってもよくわからなかったのですが、ワシントン大学 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. step 2. step 3.
では、具体的にやっていきましょう.
3-2. 使用するタンパク質の選択 実際にどのデータを使ってドッキングを行うか？記事 で比較した際に、この複合体はリガンドが大きくため、タンパク質側の結合ポケットが広がっている可能性が高い、と考えられるからです。
なぜポケットが大きい方が良いか？「タンパク質計算科学 基礎と創薬 への応用」の記述を引用します。(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".
以上の操作でだいぶスッキリした見た目になりました。
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
安易に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 ファイルを処理します。PDB ファイルを読み込んだのち、Tools の Surface /Binding Analysis から DockPrep を選択すると新しい画面が開きます。
この画面の指示に従って、
溶媒の除去 (Delete Solvent)
水素原子の付加 (Add hydrogens)
電荷 の付加 (Add charges)
といった必要な処理を選択し、Mol2ファイルで書き出します。
リガンドのファイルも同様に水素原子と電荷 の付加の処理をしました。カルボキシル基には水素原子が付加されておらず、アミノ基には付加されているので、それぞれ負電荷 、正電荷 を帯びるように処理されているようです。（下図、黄色丸部位）
3-5-2. ドッキングのジョブを投げる準備 一度「close session」した後、処理したタンパク質とリガンドのmol2ファイルを読み込みます。Surface /Binding Analysis から AutoDock Vinaを選択し、開くサブ画面に設定を書き込んでいきます。
設定の画面例
上の図で行なっている設定は、
output fileの名前（今回は5nix_redock とした）をつける
タンパク質（receptor）とリガンドを指定
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
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」とすることで一つずつに分けます。
共結晶構造中のリガンド（赤色）とともにそれぞれ描画してみます。
エネルギー的に近い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. step 2. step 3.
という流れです。
まずは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)
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)
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
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)
コンフォマ ーが４つ生成されました。py3Dmolを使った３Dの描画結果を以下に貼っておきます。
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):
o3d = rdMolAlign.GetO3A(mol,tmp)
o3d.Align()
fmParams = {}
for k in Ffactory.GetFeatureFamilies():
fparams = FeatMaps.FeatMapParams()
fmParams[k] = fparams
keep = ('Donor' ,'Acceptor' ,'NegIonizable' ,'PosIonizable' ,'Aromatic' )
featLists = []
for m in [mol, tmp]:
rawFeats = fdef.GetFeaturesForMol(m)
featLists.append([f for f in rawFeats if f.GetFamily() in keep])
fms = [FeatMaps.FeatMap(feats=x, weights=[1 ]*len (x), params=fmParams) for x in featLists]
score = fms[0 ].ScoreFeats(featLists[1 ]) / min (fms[0 ].GetNumFeatures(), len (featLists[1 ]))
return score
先のテスト化合物で得た４つのコンフォマ ーについて、Feat Mapのスコアを求めてみます。
test_scores = []
for conf in test_confs:
score = FeatMapScorer(conf, temp_8HW)
test_scores.append(score)
print (test_scores)
スコアが最大となるコンフォマ ーのみを取り出したいので、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 を表示するムービとして眺めることができました。
いくつかははみ出しているように見えますが、全体的に良い重なりとなっているように見えます。
試しにタンパク質も同時に表示させて、違う角度から眺めてみます。
アミノ酸 残基の側鎖を表示させていない状態でも、化合物がタンパク質にぶつかっている（重なっている）様子が見て取れます。やはりリガンドのコンフォマ ーに似せるだけでは、最良なポーズ、とは行かないようです。AutoDock Vinaによるドッキングは、この化合物たちの中から上手く結合ポケットに収まるようなポーズを見つけることができるでしょうか？
5. AutoDock vinaとODDTでSBVS スクリーニング対象化合物群の３次元構造の生成が完了したので、これを用いてPD-L1に対するSBVSを実施したいと思います。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)
生成された結合ポーズは４つ。リガンドの構造のリストとなっており、それぞれの構造に結合エネルギーなどの情報が付与されているようです。ODDTのドキュメントに従って情報を抜き出してみます。
type (dockin_reslut)
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' ]
Chimeraで計算した時の最小値（最安定）と同じ値となっていそうです。
残りの３つのエネルギーはChimeraの時と異なるようですが、おそらく初期構造をランダムで発生させているため、初期構造に依存して異なる結果が出ているのだと思います。
ODDTは使用方法がpybel(OpenBabel?)に近い書き方となっているそうです（使ったことがないのでわかりません・・・）。ドッキングの結果で得られるリガンドのタイプを確認してみます。
type (dockin_reslut[0 ])
この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個全て終わってました。
5-3. 全ドッキングポーズを出力・確認 全50個の化合物について、各化合物で生成されたドッキングポーズを全て格納したリストを作成します。
All_Result = []
for poses in Dock_resluts:
for pose in poses:
All_Result.append(pose)
len (All_Result)
449個のポーズが得られました。プロパティが保持されているかどうか確認してみます。
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']
"""
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で読み込めるかどうか確認してみます。
DockAllResultOutput = oddt.toolkit.Outputfile('sdf' , 'ODDT_Dock_All_Result.sdf' )
for mol in All_Result:
DockAllResultOutput.write(mol)
DockAllResultOutput.close()
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つ結合が伸びている、といったような不具合が生じているようです。
テキストエディタ でぽちぽちsdfの結合部分を修正し、再度読み込みます。
All_Dock_df = PandasTools.LoadSDF('./ODDT_Dock_All_Result_mod.sdf' )
All_Dock_df.head()
今度はエラーを出さずに読み込めました。
SDFの書き出し、構造の確認が完了したので結果のドッキングポーズをpymolで眺めてみます。
先と同様にタンパク質ともともとの共結晶構造中のリガンド、そしてドッキング結果のSDFを読み込み、ムービーを再生します。
おお！！バーチャルスクリーニング感あふれてる！！
全結合ポーズの中では、結合ポケットの端（タンパク質と溶媒の境界近く）に位置するものも多く見られ、結合ポケット奥深く（共結晶構造リガンドと重なるような位置）まで入り込んでいるものは少ないように見えます。先にみた「タンパク質計算科学」のコメント通り、結合ポケット奥深くに入るポーズは、アミノ酸 残基側鎖とのぶつかり合い(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で読み込みます。
DockBestModeResultOutput = oddt.toolkit.Outputfile('sdf' , 'ODDT_Dock_BestMode_Result.sdf' )
for mol in Best_Mode_Result:
DockBestModeResultOutput.write(mol)
DockBestModeResultOutput.close()
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 )
Dock_Top_10 = BM_df_st.iloc[:10 ]
これらドッキングスコア ベスト10の化合物群についてもPyMolのムービーで眺めて見ましょう。
おお！なんだか思っていたよりも良い感じ！！
スコアの良い結合ポーズでは、結合ポケットの奥深くまで化合物が入り込んでおり、反発を抑えつつ上手くタンパク質と相互作用ができていそうな雰囲気がただよっています。
意外な点としては、複数の化合物でビフェニル部分構造が、元の共結晶構造の該当の部分構造とは重ならない形のポーズが、最もスコアが良い結合構造として選ばれていることです。これらでは疎水性のビフェニル が親水性の溶媒接触 面側 に位置してしまっているようにも見え、ビフェニル部分構造を基準に選んだことが結合親和性に対しては逆効果に働いてしまっているようにも思えます。部分構造の結合位置を指定してドッキングする、といった細かい設定を行っていけば、より元のリガンドに近いような結合ポーズを出せるのかもしれません。
では、これらの化合物がどのような構造式のものなのか？ 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)
リガンドベースで類似性のスコアのみで評価した時とは異なる分子が取り出されているようです。両者を並べてみます。
丸で囲った部分が今回のドッキング（SBVS）と前回のLBVSで共通して選ばれている構造でTop 10のうち、３つとなりました。類似性のみのLBVSとはまた異なった特徴を持つ化合物を選んでくることができているようです。
6. まとめ 今回は、最後の難関ドッキングによるSBVSを行って見ました。色々とソフトウェアを行ったり来たりしてかなり見通しの悪い感じの記事になってしまいましたが、なんとか最終的にそれっぽい結果を出すことができました。
LBVSで絞り込んだ50個から最終的に10個を選び、LBVSのスコアとは異なるTop 10化合物となりました。SBVSの特徴として、既知の活性化合物とは構造的に異なるもの を見つけ出す可能性がある、との記載がありましたので、その意味では異なる結果が出たことでドッキングをした甲斐があった！・・・と思いたい。
（異なる構造が選びたければ先に部分構造や類似性で絞り込むのは順番的に良くないのでは？といった、つらいツッコミはご容赦いただければ幸いです。）
得られた結合ポーズは、当初想定していたのとは異なる結合様式もあった、など、まだまだ調整・工夫の余地がたくさんありそうです。ほとんどデフォルトのコピー&ペーストですし。。。とはいえ、あまりに”ヒトの目”でみて出したいポーズを狙いすぎても折角ドッキングソフト を使っている意味が薄まってしまいますので、今回はこの結果で良し！！
ということで一通りそれっぽいことができたので満足。
付け焼き刃で調べながら行なったので、間違い等多数あると思います。ご容赦、ご指摘いただければ幸いです。
P.S. 進捗悪すぎてすでに提出期限オーバー感があるのですが・・・今からの提出でも許してもらえるでしょうか？？？
