チュートリアル：基質-酵素複合体のモデリング

:::note warn
本チュートリアルはLinux PCで行われることを想定しています。また、コマンドプロンプトやcondaなどの基本的な使い方については全く触れません。
:::
　本チュートリアルではPobA（*p*-hydroxybenzoate hydroxylase）のalpha-fold model（AF model）とFAD、*p*-hydroxybenzoate（PHA）との複合体モデルを作成します。PobAはPHAの3位炭素を水酸化するFAD依存性の酸化酵素です（図1）。PobAを対象として選定した理由は以下の2つです。
1. PobA-FAD-PHA複合体の[X線結晶構造](https://www.rcsb.org/structure/1IUW)が解かれている
2. FADとPHAとでligandの異なる扱い方を学ぶことができる

図1. PobAの反応機構（一部）
![fig1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3302733/e4749094-663b-3281-8a28-bdfd6d8b6bce.png 'fig1')

# 1. AF modelをダウンロード
　[AlphaFold Protein Structure Database](https://alphafold.ebi.ac.uk/entry/P20586)からPobAのAF modelをダウンロードします。ダウンロードしたPDBファイルのファイル名を`PobA.pdb`に変更しておきます。
# 2. FADについて
　基質-酵素複合体モデルを作成するにあたり、基質の結合部位（基質ポケット）を予想する必要があります。一般に基質ポケットの予想はホモログ酵素や変異実験の結果を考慮して行います。

　PobAの反応にはFADが関与するため、FADの近傍に基質PHAが結合すると予想できます。そこで、FAD → PHA の順でポケットを予想してみます。
 
## 2-1. FADポケットを予想
　Foldseekを用いて構造ベースでホモログ酵素を検索し、それらの情報からFADポケットを予想してみます。[Foldseek](https://search.foldseek.com/search)のページでUPLOAD PDBボタンから`PobA.pdb`をアップロードし、Databases & search settingsでPDB100以外の✔を外した後、Serchボタンをクリックします。結果を見ると、配列相同性が非常に高いホモログ酵素が上位にランクインしています。それらのFADポケットを参考にしても良いのですが、今回はあえて少し配列相同性の低いホモログ酵素を参考にしてみます（図2）。一つずつ構造を確認していくと、[5xgv](https://www.rcsb.org/structure/5xgv)にFADが含まれていることがわかります。

図2. Foldseekの結果
![fig2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3302733/7c5564c2-c168-b57a-a255-20ff4f11769d.png 'fig2')

## 2-2. FADを配置
　5xgvはpyrroindomycinsの生合成に関与するPyrE3（Diels-Alderase）です。PobAのFADポケットとPyrE3のそれとが同じであると仮定し、PobA-FAD複合体の作成を試みます。

　PyrE3に結合したFADをPobAへ移植するために、これら2つのタンパク質を重ね合わせる必要があります。重ね合わせにはPyMOLのalign機能が利用できますが、Foldseekのような構造ベースで検索されたホモログ酵素は、配列相同性が高くないために、綺麗に重ならないことが多いです（図3）。そこでTM-alignを利用して、構造ベースでタンパク質を重ね合わせます（図4）。
 
1. PDBから[5xgv](https://www.rcsb.org/structure/5xgv)をダウンロードする
2. PyMOLで`5xgv.pdb`を開く
3. `5xgv → A → remove waters`により水分子を削除する
4. Command Input Areaに`sele chain B`と入力しEnter
5. `(sele) → A → extract object`により5xgvのB鎖を抽出する
6. `5xgv → A → delete object`により5xgvのA鎖を削除する
7. `File → Export Molecule… → Save…`によりファイル名`PyrE3_chainB.pdb`として保存する
8. [TM-align](https://zhanggroup.org/TM-align/)のページで2つの参照ボタンから`PobA.pdb`→`PyrE3_chainB.pdb`の順でアップロードする
9. Run TM-align ボタンをクリックする
10. 頭文字がGのPDBファイルをダウンロードする
11. ファイル名を`PobA_vs_PyrE3.pdb`に変更する

図3. PyMOL alignでの重ね合わせ（灰色:5xgv、水色:PobA）
![fig3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3302733/616a41bd-e2de-860a-3da0-5032dd4351e6.png 'fig3')

図4. TM-alignでの重ね合わせ（灰色:5xgv、緑色:PobA）
![fig4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3302733/1627c0a7-95ab-378c-2061-979c5db55995.png 'fig4')


　続いてPyrE3に結合したFADをPobAに移植します。

1. PyMOLで`PobA_vs_PyrE3.pdb`を開く
2. Command Input Areaに`sele chain A`と入力しEnter
3. `(sele) → A → extract object`によりPobAを抽出する
4. FADをクリックする
5. `(sele) → A → extract object`によりFADを抽出する
6. `PobA_vs_PyrE3 → A → delete object`によりPyrE3を削除する
7. `File → Export Molecule… → Save…`によりファイル名`PobA_FAD.pdb`として保存する

## 2-3. FADを評価
　PyrE3に結合したFADをPobAに移植した構造が妥当であるかを評価します。それにあたり、PyrE3とFADとの相互作用に関与するアミノ酸残基がPobAにおいても保存されているのかを確認します。

1. PyMOLで`PobA_vs_PyrE3.pdb`を開く
2. Command Input Areaに`sele chain A`と入力しEnter
3. `(sele) → A → extract object`によりPobAを抽出する
4. FADをクリックする
5. `(sele) → A → find → polar contacts → to other atoms in object`によりFAD-PyrE3間の相互作用を表示する
6. `(sele) → A → modify → around → residues within 4 A`によりFADの4Å以内にあるアミノ酸残基を選択する
7. `(sele) → S → licorice`により表示する

　PyrE3においてFADと相互作用するGLU32、ARG33、ARG42、ASP275が、PobAにおいても保存されていることがわかります（図5）。このことから、PyrE3に結合したFADをPobAに移植した構造は妥当であると判断します。

図5. FAD周辺のアミノ酸残基（灰色:5xgv、緑色:PobA）
![fig5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3302733/da032ce1-c6c0-3a11-d6eb-e23d9d114658.png 'fig5')


# 3. PHAについて
　FADを配置できたので、次にPHAを配置します。Foldseekの結果（[6ju1](https://www.rcsb.org/structure/6JU1)、[1d7l](https://www.rcsb.org/structure/1D7L)など）から、PHAもFADと同様の手法で配置できそうですが、今回はあえてドッキングによりPHAを配置してみます。
## 3-1. PHAをドッキング
　初めにPHAの立体構造を取得します。PyMOL等で自作できますが、今回はPubChemからダウンロードします。

1．PHAの立体構造を[PubChem](https://pubchem.ncbi.nlm.nih.gov/compound/4-Hydroxybenzoate)からSDF形式でダウンロードする
2．PyMOLでダウンロードしたSDFファイルを開く
3．`all → A → hydrogens → remove`により水素原子を削除する
4．`File → Export Molecule… → Save…`によりファイル名`PHA.pdb`として保存する

　続いてAutodock Vinaを利用し、`PHA.pdb`をligand、`PobA_FAD.pdb`をreceptorとしてドッキングします。この作業はコマンドプロンプトで行います。なお仮想環境mgltoolsには[mgltools](https://anaconda.org/bioconda/mgltools)が、仮想環autodockには[Autodock vina](https://anaconda.org/conda-forge/vina)がインストールされています。
```bash:
# 1. 仮想環境mgltools入る
conda activate mgltools

# 2. PobA_FAD.pdbをpdbqt形式に変換（receptor preparation） 
pythonsh /path/to/envs/mgltools/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_receptor4.py \
    -r PobA_FAD.pdb \
    -o receptor.pdbqt \
    -U nphs_lps_waters \
    -A hydrogens \
    -v

# 3. PHA.pdbをpdbqt形式に変換（ligand preparation） 
pythonsh /path/to/envs/mgltools/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_ligand4.py \
    -l PHA.pdb \
    -o ligand.pdbqt \
    -A hydrogens \
    -v

# 4. auto dock toolsでgrid boxを設定する
adt
```
1. 画面右下の”Remind Me Later”が赤くなり次第、クリックする
2. `File → Read Molecule`で`receptor.pdbqt`を開く
3. `Grid → Grid Box...`でgrid boxを表示する
4. grid boxの位置をFADのC4位付近まで移動し、サイズを調整する
5. メモ帳でconfig_dockを開き、grid boxの位置とサイズとを書き込む（図6）

```:config_dock
#CONFIGURATION FILE (options not used are commented)

#INPUT OPTIONS
receptor = receptor.pdbqt
ligand = ligand.pdbqt
 
#CHANGE THE FOLLOWING DATA WITH YOUR BOX DIMENSIONS
#Size of the box (values bxf, byf and bzf)
size_x = 20
size_y = 20
size_z = 20

#SEARCH SPACE CONFIGURATIONS
#CHANGE THE FOLLOWING DATA WITH YOUR BOX CENTER COORDINATES  
#Center of the box (values bxi, byi and bzi)
center_x = -21.111
center_y = 28.750
center_z = -13.972

#OUTPUT OPTIONS
#out =
#log =

#OTHER OPTIONS
cpu = 4
exhaustiveness = 8
num_modes = 10
```

図6. Grid boxの設定
![fig6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3302733/0126774b-710e-228d-e0f1-3cfa14479eed.png 'fig6')

```bash:
# 5. 仮想環autodockに入る
conda activate autodock

# 6. config_dockの条件でドッキングを実行する
vina --config config_dock --out output.pdbqt
```
ドッキングが正常に終了すると`output.pdbqt`が出力されます。

## 3-2. PHAを評価
　PobA_FADにPHAをドッキングした構造が妥当であるかを評価します。PyMOLで`receptor.pdbqt`と`output.pdbqt`とを開き、ドッキング結果を確認します。PyMOLの画面右下ある`➤`ボタンをクリックして一つずつ構造を確認したところ、ドッキング結果を3パターンに大別できることがわかりました（図7）。パターン①、②ではPHAの反応点とFADのC4a位とが近いことがわかります。また①は、PHAのカルボン酸とPobAのARGが相互作用している点から、②に比べてより合理的だと考えます。これらのことから、①がPobA_FADにPHAをドッキングした構造として妥当であると判断します。この構造を`Export Molecule…`から`PobA_FAD_PHA.pdb`として保存します。

図7. PHAのドッキング結果
![fig7.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3302733/a3d07eef-f920-23a9-5050-f2390fe42bd2.png 'fig7')


# 4. 答え合わせ
　最後に、作成したモデル`PobA_FAD_PHA.pdb`と[X線結晶構造](https://www.rcsb.org/structure/1IUW)とを比較します。両者をPyMOLのalign機能で重ね合わせると、良く一致することがわかります（図8）。

 図8. PobA_FAD_PHA.pdbとX線結晶構造との重ね合わせ（灰色:5xgv、緑色:PobA）
 ![fig8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3302733/a16a1060-634d-f8dc-12e7-1306af07a1fd.png 'fig8')

# 参考
　これらのWebページが大変参考になります。併せてご参照ください。

https://qiita.com/Ag_smith/items/96f998c9d8703959744d

https://qiita.com/Ag_smith/items/63799b3beaa07c3990fe

https://computational-chemistry.com/top/blog/2017/04/26/autodock-vina/