Modeller による Cαモデルを鋳型とした全原子モデルの構築

# Cαモデル #
　まず，この記事では下の例に示すようなCαのみのモデルを「Cαモデル」と（勝手に）呼んでいます。
　
　一般的にタンパク質の構造データは，PDBに登録されているなどして，少なくとも主鎖は完全に見えており，ファイルにもしっかりと座標が記述されています。では，どういう時にCαモデルに遭遇するかというと，様々なソフトウェアを使った結果などです。例えば，[EOM](https://www.embl-hamburg.de/biosaxs/eom.html)は結果としてPDB形式の構造データを出力しますが，モンテカルロ・サンプリングを行なった領域については，Cα以外の主鎖は取り除かれています。

``` pdb
ATOM      1  CA  MET A   1      -8.138 -34.531  18.695  1.00  0.00           C
ATOM      2  CA  GLY A   2      -4.518 -35.506  18.075  1.00  0.00           C
ATOM      3  CA  LYS A   3      -2.397 -35.029  21.191  1.00  0.00           C
ATOM      4  CA  LYS A   4      -5.259 -34.104  23.514  1.00  0.00           C
ATOM      5  CA  THR A   5      -6.053 -37.791  23.972  1.00  0.00           C
ATOM      6  CA  LYS A   6      -5.056 -38.100  27.625  1.00  0.00           C
...
```

　この例以外にも，「Cαの位置を拘束してモデリングをしたい」などの需要はあるかと思います。その際に，[MODELLER](https://salilab.org/modeller/)を使って，主鎖を補完し，側鎖を発生させて全原子モデルを構築する方法について記します。



# 結論 #
　方法自体は非常に簡単です。簡単すぎて，MODELLER の[FAQ](https://salilab.org/modeller/9v6/FAQ.html#12)にも，
>Is it possible to use templates with the coordinates for  $ {C}_\alpha$ atoms only?
>Yes. You do not have to do anything special.

と説明すらされません。
　結論としては，MODELLERを使ったことがある方なら，「テンプレートにCαモデルを指定して，モデル構築すればそれで全原子になる」という説明で十分かと思います。



# 実践：1MBN #
定番の[1mbn](https://www.rcsb.org/structure/1mbn)でまずはやってみます。MODELLERを実行するために必要な，

* テンプレートのpdbファイル
* アライメントファイル
* MODELLERスクリプト

を揃えます。

### テンプレートpdb ###

元のpdbファイルはもちろん全原子なので，あえてCαモデルのpdbファイル作成します。方法はなんでも良いですが，一例として，

``` bash
$ grep '^ATOM.*CA.*' 1mbn.pdb > ca_1mbn.pdb 
```
とでもすればCαモデルのpdbが用意できます。

### アライメントファイル ###

MODELLERでは，pdbファイルの残基とここで指定するsequenceが完全に一致しないとエラーになってしまうので，先ほど作った`ca_1mbn.pdb`から残基配列を作るのが良いと思います。

``` pir
>P1; TARGET
sequence:TARGET:1:::::::
VLSEGEWQLVLHVWAKVEADVAGHGQDILIRLFKSHPETLEKFDRFKHLKTEAEMKASEDLKKHGVTVLTALGAILKKKGHHEAELKPLAQSHATKHKIPIKYLEFISEAIIHVLHSRHPGDFGADAQGAMNKALELFRKDIAAKYKELGYQG*

>P1;ca_1mbn
structure:ca_1mbn:1:A:153:A::::
VLSEGEWQLVLHVWAKVEADVAGHGQDILIRLFKSHPETLEKFDRFKHLKTEAEMKASEDLKKHGVTVLTALGAILKKKGHHEAELKPLAQSHATKHKIPIKYLEFISEAIIHVLHSRHPGDFGADAQGAMNKALELFRKDIAAKYKELGYQG*
```

### MODELLERスクリプト ###

最小限のスクリプトです。

``` python
from modeller import *
from modeller.automodel import *

log.verbose()
env = environ()

env.io.atom_files_directory = ['.', './pdb']
a = automodel(env,
       alnfile  = 'seq.pir',
       knowns   = ['ca_1mbn'],
       sequence = 'TARGET')
a.starting_model= 1
a.ending_model  = 1

a.make()
```

### 結果 ###
　ちゃんと全原子になっています。RMSD=0.454。特に，側鎖の衝突等もみられません。黄色がオリジナルで，マゼンタが最適化なし，ピンクが最適化あり（後述）になります。最適化によって側鎖の向きが修正されているのがわかると思います。

![qiita_fig1.png](https://qiita-image-store.s3.amazonaws.com/0/175412/95766a19-492a-cb37-460d-81f3885d2ea1.png)

 すでに十分な出来栄えですが，何か問題があったり，さらに精度向上を望む場合は，MODELLERの最適化処理を入れてみましょう。上記スクリプトに追記で，

``` python
a.ending_model  = 1

# 最急降下法の最適化を，500回繰り返す。下のと合わせて，30分くらいかかる。全部でRMSD＝0.350
a.library_schedule = autosched.slow
a.max_var_iterations = 500

# MDの最適化を，very_slowに設定して，100回繰り返す。10分くらいかかる。これだけでRMSD＝0.386
a.md_level = refine.very_slow
a.repeat_optimization = 100

a.make()
```


# 実践：1C53 #

　次に，[1C53](https://www.rcsb.org/structure/1c53)で試します。pdbファイルを見てみると分かりますが，CAしか座標が記述されていない構造データです。この主鎖・側鎖を発生させて全原子モデルを構築していきます。なお正解は分かりません。アライメント，スクリプトは，1mbnの例をそのまま1C53のデータに置き換えたものを使います（最適化はなし）。
　また，比較として[REMO](https://zhanglab.ccmb.med.umich.edu/REMO/)を使います。REMOの使い方は，`perl remo.pl ca_model.pdb seq.dat`とするだけなので，非常に簡単で高速です。

下図は，REMOの結果。なんとなく気持ち悪い箇所ができてしまう。この辺りは，REMO→側鎖除く→SCWRLなどとすればより精度よくできるかもしれません。

![qiitafig4.png](https://qiita-image-store.s3.amazonaws.com/0/175412/d6e8de3e-35d6-f2ac-1bcf-eaa553227a4f.png)


一方，MODELLERでは，至って平凡な全原子モデルを構築することができました。とはいえ，MODELLERは，スクリプトやアライメントファイルを用意しなければならないので，精度なりの手間がかかってくることは否めません。その辺りは，REMOとどちらが適切かケースバイケースかと思います。