チュートリアル：AutoDock4でCovalent Docking

　[1pwc](https://www.rcsb.org/3d-view/1PWC)はβ-lactamaseのSer側鎖とpenicillin Gとが共有結合した複合体です。本記事ではAutoDock4 (ad) のcovalent dockingで、β-lactamaseの*apo*体から1pwcを再現します。

　adでcovalent dockingを行うにはMGLToolsに含まれていないpythonスクリプトを使用する必要があります。そこで[Covalent Docking - AutoDock](https://autodock.scripps.edu/resources/covalent-docking/)から`adCovalentDockResidue_v1.2.tar.gz`をダウンロードし、コマンドプロンプトで解凍します。本記事は、解凍後のフォルダの`README.txt`に沿ってcovalent dockingします。
```bash
tar -xvf adCovalentDockResidue_v1.2.tar.gz
```

# 1. ファイルの準備
　ドッキングに必要なligandとreceptorとを準備します。ligandのpenicillin GはPubChemから3次元構造をダウンロードしてきます。またreceptorのβ-lactamaseはPDBからダウンロードしてきます。
## 1-1. receptorの準備
　PDBからβ-lactamaseの*apo*体である[3PTE](https://www.rcsb.org/3d-view/3PTE)をダウンロードします。PyMOL等で`3pte.pdb`から結晶水を消去し、`receptor.pdb`として保存します。
## 1-2. ligandの準備
　結合する残基がSerである場合、構造にMeOHを付加したligandを作成します。初めに[PubChem](https://pubchem.ncbi.nlm.nih.gov/compound/Penicillin-G)から3次元構造（sdf形式）をダウンロードし、ファイル名を`PNM.sdf`に変更します。続いてPyMOL等でMeOHが付加した開環体をモデリングし、`PNM_MeOH.mol2`（mol2形式）として保存します。ダウンロードしてきた3次元構造ではプロトン化状態が中性付近のものでないため、その修正もここでしておきます。最後に`PNM_MeOH.mol2`をChem3D等でエネルギー最小化し、`PNM_MeOH_opt.mol2`として保存します（図1）。

 図1. PNM.sdfとPNM_MeOH_opt.mol2
 ![fig1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3302733/94bc008c-e04c-f356-f62f-b1778faefc51.png)


# 2. prepareCovalent.py
　先ほどダウンロードしたpythonスクリプト（`prepareCovalent.py`）で、receptorのSer側鎖とligandに付加したMeOHと重ね合わせます。

　まず、`PNM_MeOH_opt.mol2`においてMeOHのC, Oが何番目の原子であるのかを確認します（図2）。Cが6番目、Oが9番目であることから、`--ligindices`を`6,9`とします。

図2. PNM_MeOH_opt.mol2のMeOH部分
![fig2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3302733/1e441e8d-e97f-f5b7-713e-230d13466fc7.png)


　続いて、`receptor.pdb`においてどのSer残基がligandと共有結合するかを確認します。先行研究から`receptor.pdb`はChainA 62番目のSer側鎖でpenicillin Gと結合していることから、`--residue`を`A:SER62`とします。

　[mgltools](https://anaconda.org/bioconda/mgltools)をインストールした仮想環境で以下を実行すると`ligcovalent.pdb`が出力されます。
 
```bash
pythonsh /path/to/adCovalentDockResidue/adcovalent/prepareCovalent.py \
                  --ligand PMN_MeOH_opt.mol2 \
                  --ligindices 6,9 \
                  --receptor receptor.pdb \
                  --residue A:SER62 \
                  --outputfile ligcovalent.pdb
```

　`receptor.pdb`と`ligcovalent.pdb`とを描画すると、receptorとligandとが衝突していることがわかりますが、ここでは問題ありません（図3）。

図3. receptor.pdbとligcovalent.pdbとの衝突
![fig3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3302733/ecd31bce-f2b3-6661-7ae0-c9587bdfe14d.png)




# 3. PDBQTへの変換
　`receptor.pdb`と`ligcovalent.pdb`とをadで取り扱うためPDBQT形式に変換します。mgltoolsをインストールした仮想環境で以下を実行すると、それぞれをpdbqt形式に変換したファイルが出力されます。
```bash
# receptor.pdbの変換 (receptor.pdbqtが出力される)
pythonsh /path/to/envs/mgltools/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_receptor4.py \
                  -r receptor.pdb \
                  -o receptor.pdbqt \
                  -U nphs_lps_waters \
                  -v

# ligcovalent.pdbの変換 (ligcovalent.pdbqtが出力される)
pythonsh /path/to/envs/mgltools/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_ligand4.py \
                  -l ligcovalent.pdb \
                  -o ligcovalent.pdbqt \
                  -v 
```
# 4. flexible/rigid PDBQTの生成
　adのCovalent Dockingではligandをflexible residueとして扱います。そのためにreceptor.pdbqtとligcovalent.pdbqtとを、それぞれflexible partとrigid partとに分けます。
```bash
# receptor.pdbqtの変換 (receptor_rigid.pdbqt、receptor_flex.pdbqtが出力される)
pythonsh /path/to/envs/mgltools/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_flexreceptor4.py \
                  -r receptor.pdbqt \
                  -s receptor:A:SER62

# ligcovalent.pdbqtの変換 (ligcovalent_rigid.pdbqt、ligcovalent_flex.pdbqtが出力される)
pythonsh /path/to/envs/mgltools/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_flexreceptor4.py \
                  -r ligcovalent.pdbqt \
                  -s ligcovalent:A:SER62
```
　Covalent dockingにはreceptor.pdbqtのrigid partとligcovalent.pdbqtのflexible partとを使用します（図4）。
 
図4. AutoDock4のCovalent dockingで使用するファイル
![fig4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3302733/950c98e8-be7d-ca05-26a7-88b7fb82ea0c.png)


# 5. パラメータファイルの生成
　ドッキング条件を指定するパラメータファイルを出力します。以下を実行すると、
 `receptor.gpf`と`ligcovalent_receptor.dpf`が出力されます。
```bash
# .gpf (AutoGrid用のパラメーターファイル) の生成
pythonsh /path/to/envs/mgltools/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_gpf4.py \
                -r receptor_rigid.pdbqt\
                -x ligcovalent_flex.pdbqt\
                -l ligcovalent_flex.pdbqt\
                -y -I 20\
                -o receptor.gpf

# .dpf (AutoDock用のパラメーターファイル) の生成
pythonsh /path/to/envs/mgltools/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_dpf4.py \
                -r receptor_rigid.pdbqt\
                -x ligcovalent_flex.pdbqt\
                -l ligcovalent_flex.pdbqt\
                -o ligcovalent_receptor.dpf\
                -p move='empty'

# 空のファイルを作成しておく。
touch empty
```
　生成された`ligcovalent_receptor.dpf`をメモ帳等で開き、`unbound_model bound`を`unbound_energy 0.0`に書き換えます。さらに`ga_run`の値を10から100に書き換え、出力されるドッキングポーズを100個にします。

# 6. ドッキング
　最後にドッキングします。[autogrid4](https://anaconda.org/bioconda/autogrid)と[autodock4](https://anaconda.org/bioconda/autodock)とをインストールした仮想環境で以下を実行し正常に終了すると、ドッキングモデルが記載された`ligcovalent_receptor.dlg`が出力されます。
```bash
autogrid4 -p receptor.gpf -l receptor.glg
autodock4 -p ligcovalent_receptor.dpf -l ligcovalent_receptor.dlg
```
　結果は`adt`から確認します。mgltoolsをインストールした仮想環境でコマンドプロンプトに`adt`と入力し実行すると、AutoDockToolsの画面が現れます。

1. 画面右下の”Remind Me Later”が赤くなり次第、クリックする
2. `Analyze → Docking → Open...`により`ligcovalent_receptor.dlg`を開き、ドッキングポーズを表示する
3. `Analyze → Macromolecule → Open...`により`receptor_rigid.pdbqt`を開き、receptorを表示する
4. `Analyze → Conformations → Play, ranked by energy...`によりドッキングポーズをエネルギーが低い順にソートする
5. **一つずつドッキングポーズを確認し、確からしいポーズを選択する**
6. `右から2番目のボタン → Write Complex`から表示しているポーズを保存する（図5）。

図5. ドッキングポーズを保存するときの操作
![fig5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3302733/f54d36a8-b036-6bc1-ebe1-973e85e8109f.png)


# 7. 答え合わせ
　保存したドッキングポーズと[X線結晶構造](https://www.rcsb.org/3d-view/1PWC)とを比較します。それらをPyMOLのalign機能で重ね合わせて答え合わせします。今回のドッキングシミュレーションでは出力した100ポーズの内、rank91とrank100とがX線結晶構造に近いということがわかりました。しかしrank91ではエステル結合の向きが、rank100ではカルボン酸の位置が1pwcと異なっています。また今回のように、設定条件にも依りますが、正しいポーズのランクが必ずしも良いとは限りません。

図6. rank91、rank100、X線結晶構造の重ね合わせ（灰色:1pwc、緑色:rank91、黄色:rank100）
![fig6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3302733/4279e499-dfa4-1be7-d373-53c3aaa4519a.png)