Docking simulation

Updated on October 19, 2016.

AutoDock4により低分子化合物とタンパク質とのドッキングシミュレーションを行う.AutoDockは標的タンパク質のリガンド結合部位をグリッドとして表現し,このグリッドに対するリガンドの結合様式を遺伝的アルゴリズムやMonte Carlo Simulated Annealing法などにより最適化する.
まずタンパク質側鎖を動かさないrigid dockingの方法を解説する.
AutoDock4およびMGL Toolsをインストールして,PATHを通しておく.

Step 1: Building a small molecule

低分子の立体構造の作成にはAvogadroが使用できる.  ドッキングでは環構造は変化しないため,フレキシブルな環構造を持つ化合物についてはconformation searchを行う必要がある.
今回は,1PTR.pdbファイルから抽出したphorbol-13-acetate (PRBH.pdb,水素付加済み) を元のタンパク質結晶構造 (1PTR_apo.pdb) に対してドッキングし,シミュレーションの妥当性を検証する.1PTRはアミノ酸側鎖を3残基欠いているが,化合物の結合部位からは離れているため,修正せずにドッキングに使用する.

Step 2: Preparing Protein

作業ディレクトリ (日本語がPATHに含まれていると後でエラーが出るので,上位フォルダにも日本語を含まないこと) でPMVを起動.

$ pmv &
[Menu] File > Read Molecule > Open...  1PTR_apo.pdb
[Menu] Edit > Hydrogens > Add > OK (All Hydrogens, noBondOrder, Yes)
[Menu] File > Save as > Write PDB > Filenameを1PTR_H.pdbとして > OK

Step 3: Preparing Ligand

上部アイコン右から3番目のAutoDock ToolsボタンをクリックするとADT4.2メニューがあらわれる.

[ADT Menu] Ligand > Input > Open...  > PRBH.pdb > OK

non-polarな水素原子は削除され,PDBファイルを開いた場合はgasteiger chargeが付加される.部分電荷情報を持つmol2ファイルを使用した場合は,その電荷情報が使用されるため,予めMOPACやGAMESS等で構造最適化等を行い,結果ファイルをopenbabelで電荷情報を含むmol2形式に変換してもよい.
MOPACを利用した例(構造最適化):

$ babel PRBH.pdb PRBH.mop

PRBH.mopの1行目を “PM7 BONDS CHARGE=0 SINGLET XYZ PRNT=2 PRECISE GNORM=0.0”に 書き換える.

$ mopac PRBH.mop

30秒程度で終了.出力ファイルをmol2に変換する.

$ babel -imopout PRBH.out PRBH.mol2

AutoDock Toolsでリガンドを読み込んだ際に,シクロプロパン環の炭素原子が芳香族炭素であると誤判定されるので,保存した後に手動で修正する.

[ADT Menu] Ligand > Torsion Tree > Detect root...
[ADT Menu] Ligand > Output > Save as PDBQT...

テキストエディタでPRBH.pdbqtを開く.

ATOM     14  C   LIG d   1      11.019  29.040  23.863  0.00  0.00     0.072 A 
ATOM     15  C   LIG d   1      11.408  28.203  25.066  0.00  0.00    -0.050 A 
ATOM     16  C   LIG d   1      10.439  29.377  25.237  0.00  0.00     0.038 A

A (芳香族)と判定されている炭素原子をCに書き換えて保存.

ATOM     14  C   LIG d   1      11.019  29.040  23.863  0.00  0.00     0.072 C 
ATOM     15  C   LIG d   1      11.408  28.203  25.066  0.00  0.00    -0.050 C 
ATOM     16  C   LIG d   1      10.439  29.377  25.237  0.00  0.00     0.038 C

PMV上のリガンドは一旦削除して,修正したPDBQTファイルを再度読み込む.

[Menu] File > Read Molecule > Open... > PRBH.pdbqt

Step 4: Creating grid parameter file (gpf)

[ADT Menu] Grid > Macromolecule > Choose... > 1PTR_apo > Select Molecule > OK > Save > OK

Non-polarな水素原子が削除され,PDBQTファイルが保存される.Znイオンの電荷が0と警告が出る.PDBQTファイルをテキストエディタで開くと,Znイオンの電荷は0で,システインの硫黄原子の電荷もおかしいが, 化合物の結合部位からは離れているため修正しない.

[ADT Menu] Grid > Set Map Types > Choose Ligand... > PDBH > Select Ligand
[ADT Menu] Grid > Grid Box...

リガンド結合cleftが全て含まれるようにBoxのサイズと位置を調整する.

[Grid Options Menu] File > Close saving current
[ADT Menu] Grid > Output > Save GPF.. > File name “R.gpf” > Save

Step 5: Creating docking parameter file (dpf)

[ADT Menu] Docking > Macromolecule > Set Rigid Filename... > 1PTR_apo.pdbqt > Open
[ADT Menu] Docking > Ligand > Choose... > PRBH > Select Ligand > Accept
[ADT Menu] Docking > Search Parameters > Genetic Algorithm... > Number of GA RunsやMaximum Number of generationsなどの値を設定(今回は100, Small) > Accept
[ADT Menu] Docking > Output > Lamarckian GA(4.2)... > File name “R.dpf” > Save

Step 6: Running AutoGrid

[ADT Menu] Run > Run AutoGrid... >

ログのファイル名が設定されていないので,

autogrid4 -p ./R.gpf -l  &

autogrid4 -p ./R.gpf -l ./R.glg &

に書き換えてから,Launchをクリック.すぐ終了し,グリッドマップファイルが生成する.あるいは,上記のコマンドをTerminalから実行してもよい.

Step 7: Running AutoDock

[ADT Menu] Run > Run AutoDock... >

ログのファイル名が設定されていないので,

autodock4 -p R.dpf -l  &

autodock4 -p R.dpf -l R.dlg &

に書き換えてから,Launchをクリック.あるいは上記コマンドをTerminalから実行する.
もしAutoDock will say “autodock4: FATAL ERROR: Sorry, I can’t find or open AD4.1_bound.dat”.
というエラーが出た場合は,AD4.1_bound.datを探してworking directoryにコピーする.

$ cp /Library/MGLTools/1.5.7rc1/MGLToolsPckgs/AutoDockTools/AD4.1_bound.dat .

再度Autodock4を実行する.

$ autodock4 -p R.dpf -l R.dlg &

今回はMaximum Number of evaluationsがshort,Number of GA Runsが100なので3分程度で終了.Autodock vinaと異なりAutodockは並列で実行できない.

Step 8: 結果の解析

R.dlgをテキストエディタで開き,”HIST”と検索してファイルの最後部のヒトグラムの結果を見る.

________________________________________________________________________________
     |           |     |           |     |                                    
Clus | Lowest    | Run | Mean      | Num | Histogram                          
-ter | Binding   |     | Binding   | in  |                                    
Rank | Energy    |     | Energy    | Clus|    5    10   15   20   25   30   35
_____|___________|_____|___________|_____|____:____|____:____|____:____|____:___
   1 |     -6.66 |  99 |     -6.14 |  78 |##############################################################################
   2 |     -5.69 |  46 |     -5.38 |  15 |###############
   3 |     -5.50 |  87 |     -5.30 |   7 |#######
_____|___________|_____|___________|_____|______________________________________

Binding Energyと分布が確認できる.Energyよりも分布の方を重視する人もいるとのこと (リンク先PDF p.11).
Autodock FAQ: How do I evaluate AutoDock’s clustering results? のまとめ

  • 有意なサンプル数を得るため,ドッキングは少くとも50 runs,出来れば100 runsを行う.
  • リガンドのサイズに合ったRMSD tolerance (rmstol) を設定する (大きな分子の場合少くとも2オングストローム).
  • Docking search (ga_num_evals) が十分長いかどうか.これはリガンドやタンパク質側鎖の回転の自由度に依存する.理想的にはsearchが十分長ければ,最低エネルギーの1つのクラスターのみが得られる.
  • 判断が付かない場合は,結合様式を目で見てどちらが化学的に妥当かを考える.
  • クラスターの平均エネルギーの差が約2.5 kcal/mol以下の場合,AutoDock force fieldの標準偏差以内であり,エネルギーからどちらが「正しい」かをいうことは出来ない.

Step 5: 結果の可視化,複合体構造の書き出し

[ADT Menu] Analyze > Dockings > Open... > R.dlg > Open

後の複合体書き出しのためにMacromoleculeを選んでおく.

[ADT Menu] Analyze > Macromolecule > Choose... > 1PTR_apo > Select Macromolecule

クラスターのグラフ等を見ることもできる.

[ADT Menu] Analyze > Conformations > Play, ranked by energy...

新たに開かれたwindowの右から2番目の ”&” ボタンをクリック.

[Set Play Options  Window] Write Complex > 名前 > Save

保存したPDBQTファイルはpolarな水素原子しか含んでいないので,再度開き,水素原子を付加した後にPDB形式で保存する.