Amber Tutorial "1.4 Building Protein Systems in Explicit Solvent"の実施

#1. 参考情報
https://ambermd.org/tutorials/basic/tutorial7/index.php

#2. 環境
PC構成

* OS：CentOS 7.9.2009
* CPU：Intel Core i7-4770
* グラフィックボード：Nvidia GeForce GTX1650
* メモリ：16GB
* 外付けSSD:240GB

ソフト

* Ambertools20(Anacondaでインストール)

#3. PDBファイルの用意

https://www.rcsb.org/structure/2YX8

からPDBファイルをダウンロード

#4. PDBデータの修正

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/788299/a65f4573-0a1a-b9e0-ecb1-c192e06b032a.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/788299/eb9449a8-a24d-e13a-5720-64f923ec9dd6.png)


##4.1 MSE残基の変更

このPDBデータは結晶構造を取得する際にmulti-wavelength anomalous dispersion (MAD) という手法が使われていて、methionineがselenomethionineになっているそうなので、そこを修正する。

``` 
HETATM  168  N   MSE A  48       6.674   9.818  33.156  1.00 45.09           N  
HETATM  169  CA  MSE A  48       6.049  10.959  32.490  1.00 44.83           C  
HETATM  170  C   MSE A  48       4.687  10.621  31.903  1.00 44.85           C  
HETATM  171  O   MSE A  48       3.782  11.457  31.884  1.00 44.46           O  
HETATM  172  CB  MSE A  48       6.968  11.508  31.392  1.00 42.34           C  
HETATM  173  CG  MSE A  48       8.206  12.210  31.932  1.00 37.42           C  
HETATM  174 SE   MSE A  48       7.764  13.636  33.161  1.00 33.09          SE  
HETATM  175  CE  MSE A  48       7.589  15.051  31.873  1.00 30.88           C  
```
これを変更する

* HETATMをATOMに
* 残基名MSEをMETに
* 元素名SEをSに


```
ATOM    168  N   MET A  48       6.674   9.818  33.156  1.00 45.09           N  
ATOM    169  CA  MET A  48       6.049  10.959  32.490  1.00 44.83           C  
ATOM    170  C   MET A  48       4.687  10.621  31.903  1.00 44.85           C  
ATOM    171  O   MET A  48       3.782  11.457  31.884  1.00 44.46           O  
ATOM    172  CB  MET A  48       6.968  11.508  31.392  1.00 42.34           C  
ATOM    173  CG  MET A  48       8.206  12.210  31.932  1.00 37.42           C  
ATOM    174  SD  MET A  48       7.764  13.636  33.161  1.00 33.09           S  
ATOM    175  CE  MET A  48       7.589  15.051  31.873  1.00 30.88           C  
```

もう一箇所MSE残基があるので、そちらも変更する

##4.2 ジスルフィド結合の指定

システイン残基間がジスルフィド結合している。
PDBファイルの最初の方にその旨が書かれている

```
SSBOND   1 CYS A   27    CYS A   82                          1555   1555  2.03  
SSBOND   2 CYS A   40    CYS A   72                          1555   1555  2.03  
SSBOND   3 CYS A   57    CYS A  104                          1555   1555  2.03  
```

ジスルフィド結合は3つあり、それぞれ残基27と82、40と72、57と104が結合している

```
ATOM      1  N   CYS A  27      34.939   8.141  26.097  1.00 61.68           N  
ATOM      2  CA  CYS A  27      34.305   6.810  26.338  1.00 61.01           C  
ATOM      3  C   CYS A  27      33.741   6.231  25.060  1.00 63.98           C  
ATOM      4  O   CYS A  27      33.961   6.755  23.970  1.00 63.63           O  
ATOM      5  CB  CYS A  27      33.143   6.934  27.318  1.00 56.61           C  
ATOM      6  SG  CYS A  27      31.731   7.869  26.631  1.00 47.29           S  
```

CYSをCYXに変更する

```
ATOM      1  N   CYX A  27      34.939   8.141  26.097  1.00 61.68           N  
ATOM      2  CA  CYX A  27      34.305   6.810  26.338  1.00 61.01           C  
ATOM      3  C   CYX A  27      33.741   6.231  25.060  1.00 63.98           C  
ATOM      4  O   CYX A  27      33.961   6.755  23.970  1.00 63.63           O  
ATOM      5  CB  CYX A  27      33.143   6.934  27.318  1.00 56.61           C  
ATOM      6  SG  CYX A  27      31.731   7.869  26.631  1.00 47.29           S  
```



##4.3 protnation stateの指定

H++にアクセスしてプロトン化状態を推定する
http://biophysics.cs.vt.edu/index.php

process a structureをクリック
CYSを修正する前のPDBファイルを選択して「Process File」をクリック
pHに7.0を入力して「Process」をクリック
出力されたファイルをVMDなどで可視化すると、75番目のHISは両方プロトン化されており、（＝HIP）
97番目のHISはdelta位がプロトン化されている（＝HID）ことがわかる

元のPDBファイルを編集してHISの残基名を変更する

```
ATOM    394  N   HIP A  75      25.140   8.660  35.933  1.00 39.77           N  
ATOM    395  CA  HIP A  75      25.640   7.344  36.319  1.00 42.09           C  
ATOM    396  C   HIP A  75      26.176   6.605  35.097  1.00 40.49           C  
ATOM    397  O   HIP A  75      27.259   6.022  35.141  1.00 40.59           O  
ATOM    398  CB  HIP A  75      24.537   6.501  36.971  1.00 46.71           C  
ATOM    399  CG  HIP A  75      24.139   6.972  38.335  1.00 52.91           C  
ATOM    400  ND1 HIP A  75      25.003   7.650  39.171  1.00 56.38           N  
ATOM    401  CD2 HIP A  75      22.980   6.836  39.023  1.00 55.42           C  
ATOM    402  CE1 HIP A  75      24.392   7.912  40.312  1.00 58.17           C  
ATOM    403  NE2 HIP A  75      23.163   7.428  40.249  1.00 58.09           N  
```


##4.4 CONECTの削除
PDBファイル中の「CONECT」から始まる行を削除する

##ファイル名の変更
上記の変更を行ったファイルを「RAMP1.pdb」とします。

以上で準備は完了


#5. LEaPで処理

tleap.inを作成します

```
#tleap.in
#タンパク用の力場ff19SBの読み込み
source leaprc.protein.ff19SB
#水の力場OPCの読み込み
source leaprc.water.opc
#イオンの力場の読み込み
loadamberparams frcmod.ions1lm_126_hfe_opc
#PDBファイルの読み込み
ramp=loadpdb RAMP1.pdb
#ジスルフィド結合の指定
bond ramp.27.SG ramp.82.SG
bond ramp.40.SG ramp.72.SG
bond ramp.57.SG ramp.104.SG
#タンパク単体の状態で一旦パラメータファイルと座標ファイルを出力
SaveAmberParm ramp RAMP1_gas.prmtop RAMP1_gas.inpcrd
#中和するためにイオンの追加
addIons ramp Na+ 2
#溶媒を追加。形はOctahedron。Proteinから壁までの距離が10Å
solvateOct ramp OPCBOX 10.0
#bufferとなるイオンを追加
#150 mM NaCl salt buffer を再現するため、ボックスの体積からイオンの数を計算して加える
addIonsRand ramp Na+ 19 Cl- 19 
#溶媒とイオンを加えた状態で一旦パラメータファイルと座標ファイルを出力
SaveAmberParm ramp RAMP1_ion.prmtop RAMP1_ion.inpcrd
```
あとは、

```
tleap -f tleap.in
```
で実行。
もちろん

```
tleap
```
でインタラクティブに実行してもOK

ジスルフィド結合は、残基名をCYXにするだけでは、どことどこを結合するのか判断してくれない
（確かCYXの残基が2つだけだったら判断してくれた）

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/788299/fd53e537-411c-6be1-189c-2fa6a138f618.png)



#6. おまけ：GROMACS形式への変換
convert_to_gro.pyを作成

```python
import parmed as pmd

molname = "RAMP1_ion"
prmtop = molname + ".prmtop"
inpcrd = molname + ".inpcrd"
top = molname + ".top"
gro = molname + ".gro"
amber = pmd.load_file(prmtop, xyz=inpcrd)
amber.save(top)
amber.save(gro)
```

実行

```
python convert_to_gro.py
```

これでgroおよびtopファイルが出力されるので、GROMACS計算ができる