失敗しないGROMACSトラジェクトリの変換方法

2019年4月25日追記：

一部の記述を修正。

## 目的

  - GromacsのMD計算で得られたトラジェクトリファイルを適切に変換し、VMDやPyMOLなどで表示するための知識やメモ。
  - 特に`gmx trjconv`コマンド周辺を使った変換方法の説明。
  - 水分子・溶媒イオンのトラジェクトリを必要としない場合、**AmberToolsを使う方法の方が速い**。

## 前提知識

  - Gromacsのトラジェクトリファイルは`trr`形式と`xtc`形式がある。`xtc`形式はMDの設定（`.mdp`ファイル内）に指定していなければ自動的に生成されない。`trr`は各粒子の座標＋速度＋時間情報を保存できるのに対し、`xtc`ファイルは座標＋時間情報のみ。
  - Gromacsのトラジェクトリファイルは原子の種類や残基の種類などの情報を含まない。反対に、これらの情報は`pdb`、`gro`、`tpr`形式のファイルなら通常含んでいる。
      - したがって、VMDやPyMOLなどでタンパク質を見たい場合には、先に原子情報・残基情報を持っている`pdb`、`gro`を開いてタンパク質を構成する基本情報をロードし、その上で各タンパク質の座標の変化情報を持つトラジェクトリファイルをロードする必要がある。**このとき、`pdb`、`gro`を構成する粒子数と、トラジェクトリファイルを構成する粒子数が一致しないと、たいていエラーを吐かずにバグるので、常にこれらが一致しているかを確認する必要がある。**
  - Gromacsコマンドは基本的に、拡張子でファイルタイプを判別する。中身を見て自動判定してくれることはしないはず。
  - **Gromacsの周期境界条件（PBC）内でMDシミュレーションを行う場合、Gromacsトラジェクトリファイル内に保存される座標はそのPBCの外側に出ない。**
      - なので、何もしないまま`trr`ファイルをPyMOLで見ようとすると、こんな感じでタンパク質が途切れて見えることがほとんどである（内部処理的にはPBCでループしているため壊れていないが、人間には壊れているように見える）  
        ![1lkemovie0.gif](https://qiita-image-store.s3.amazonaws.com/0/224327/d51b73e9-75bd-e84e-e48b-a8ad02568f3f.gif)

  - 以上から、Gromacsで計算されたMDトラジェクトリを人間の目で理解しやすくするためには、以下の操作を行うことが多い。
      - タンパク質が途切れて表示されるのを修正する
      - タンパク質の重心を平行移動操作によって固定＆センタリング
      - さらに、回転操作を用いてあるリファレンス構造から大きく回転しないようフィッティングさせる
      - 場合によっては、水やイオンを見る必要はないのでこれらを適切に除去する

## 実際の操作

Gromacsで得られた`trr`ファイルを`gmx trjconv`で処理を施し、適切に見えやすくする。私が色々試した結果、トラジェクトリの変換操作は、

1.  PBCによる途切れを全体的に修正する（whole）
2.  さらにそのトラジェクトリをclusterごとに集める
3.  続いてタンパク質（溶質）の重心の位置をセンタリングする
4.  さらに回転操作を行う

の順に変換を行うと、最も普遍的に問題なく処理を行うことができるようだ。これらの変換操作を一度に行うことはできない様子。`gmx trjconv`によるトラジェクトリ変換にはリファレンス構造を必要とし、`gmx trjconv`の`-s`に続くオプションで指定する。通常`tpr`形式が最も都合が良い。このとき、**途切れた状態で保存されている`tpr`形式を`-s`入力に使うと表示がおかしくなるので十分に注意する**。したがって、`-s`に入れるリファレンス構造は、トポロジーを作ってエネルギー最小化を始める前の`tpr`ファイルが最も都合が良い。

### インデックスファイルの作成

`gmx trjconv`はトラジェクトリの変換を行ってくれるプログラムであり、引数`-n`にGROMACS形式のインデックスファイル（拡張子`.ndx`, 通常`index.ndx`とされることが多い）を指定することができる。インデックスファイルを指定しない場合でも最低限のグループ選択はできるような仕組みになっているが、インデックスファイルを作成しておくとグループ選択を自由に追加できるので、やっておいたほうが良い。

インデックスファイルを作るコマンドは`gmx make_ndx`である。さらに`-f`で何かしらの構造ファイルを与える。以下はここで用いるインデックスファイルの作成例である。タンパク質（PDB:1LKE）のMDシミュレーションで用いる座標ファイル`1lke.gro`が存在したとして、これに対し

    gmx make_ndx -f 1lke.gro

と入力する。すると続いて追加入力を求められる。

    There are:  8527      Water residues
    Analysing Protein...
    Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...
    Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...
    
      0 System              : 28284 atoms
      1 Protein             :  2605 atoms
      2 Protein-H           :  1338 atoms
      3 C-alpha             :   169 atoms
      4 Backbone            :   507 atoms
      5 MainChain           :   677 atoms
      6 MainChain+Cb        :   830 atoms
      7 MainChain+H         :   842 atoms
      8 SideChain           :  1763 atoms
      9 SideChain-H         :   661 atoms
     10 Prot-Masses         :  2605 atoms
     11 non-Protein         : 25679 atoms
     12 Other               :    62 atoms
     13 DOG                 :    62 atoms
     14 Na+                 :    20 atoms
     15 Cl-                 :    16 atoms
     16 Ion                 :    36 atoms
     17 DOG                 :    62 atoms
     18 Na+                 :    20 atoms
     19 Cl-                 :    16 atoms
     20 Water               : 25581 atoms
     21 SOL                 : 25581 atoms
     22 non-Water           :  2703 atoms
     23 Water_and_ions      : 25617 atoms
    
     nr : group      '!': not  'name' nr name   'splitch' nr    Enter: list groups
     'a': atom       '&': and  'del' nr         'splitres' nr   'l': list residues
     't': atom type  '|': or   'keep' nr        'splitat' nr    'h': help
     'r': residue              'res' nr         'chain' char
     "name": group             'case': case sensitive           'q': save and quit
     'ri': residue index
    
    >

ここでは`>`を求められるたびに、以下のように追加入力を行う。`q`を押すと終了。

``` 
 23 Water_and_ions      : 25617 atoms

 nr : group      '!': not  'name' nr name   'splitch' nr    Enter: list groups
 'a': atom       '&': and  'del' nr         'splitres' nr   'l': list residues
 't': atom type  '|': or   'keep' nr        'splitat' nr    'h': help
 'r': residue              'res' nr         'chain' char
 "name": group             'case': case sensitive           'q': save and quit
 'ri': residue index

> !23

Copied index group 23 'Water_and_ions'
Complemented group: 2667 atoms

 24 !Water_and_ions     :  2667 atoms

> name 24 nowation


> q
```

ここではnowationという名前で「タンパク質＋リガンドDOG（化合物名digoxigenin）」のグループを作成したいとする。これは「（系全体）-（水・イオン）」のグループと等価なので、この設定でグループを作成する。そこで追加入力の一度目は`!23`を入力する。これは`23: Water_and_ions`に含まれる原子の補集合を表す`!`をつけることを意味する。これを指定すると上記のようなComplemented group 2667 atomsが`24 !Water_and_ions : 2667 atoms`として出てくる（番号は系によって異なる）。この24番のグループを改名したい場合は`name 24 nowation`とする。最後に`q`を入力して終了する。

こうしてインデックスファイル`index.ndx`の24番に「タンパク質＋リガンドDOG（化合物名digoxigenin）」のグループが作成された。以降ではこうして作成した`index.ndx`がディレクトリ`../../top/index.ndx`にあるものとして操作を行う。

### 1段階目の変換操作

まず1段階目の変換は

    input=1lke
    gmx trjconv -f ${input}.trr -o ${input}nj.trr -s ../minimize/min1.tpr -n ../../top/index.ndx -b 0 -e 20000 -skip 100 -pbc whole

ここで重要なのは、全体のトラジェクトリの途切れを修正する`-pbc whole`の部分である。`-b`と`-e`のオプションではそれぞれ変換したいトラジェクトリの時間帯を指定できる。単位はpsなので、20000 ps = 20 ns。これは`-tu`オプションのデフォルトがpsになっているせい。詳細は`gmx trjconv -h`で確認できる。
これを実行すると、続けてこんな画面が出て追加入力を促される。

    Command line:
      gmx trjconv -f 1lke.trr -o 1lkenj.trr -s ../minimize/min1.tpr -n ../../top/index.ndx -b 0 -e 20000 -skip 100 -pbc whole
    
    Will write trr: Trajectory in portable xdr format
    Reading file ../minimize/min1.tpr, VERSION 2018 (single precision)
    Reading file ../minimize/min1.tpr, VERSION 2018 (single precision)
    Select group for output
    Group     0 (         System) has 28284 elements
    Group     1 (        Protein) has  2605 elements
    Group     2 (      Protein-H) has  1338 elements
    Group     3 (        C-alpha) has   169 elements
    Group     4 (       Backbone) has   507 elements
    Group     5 (      MainChain) has   677 elements
    Group     6 (   MainChain+Cb) has   830 elements
    Group     7 (    MainChain+H) has   842 elements
    Group     8 (      SideChain) has  1763 elements
    Group     9 (    SideChain-H) has   661 elements
    Group    10 (    Prot-Masses) has  2605 elements
    Group    11 (    non-Protein) has 25679 elements
    Group    12 (          Other) has    62 elements
    Group    13 (            DOG) has    62 elements
    Group    14 (            Na+) has    20 elements
    Group    15 (            Cl-) has    16 elements
    Group    16 (            Ion) has    36 elements
    Group    17 (            DOG) has    62 elements
    Group    18 (            Na+) has    20 elements
    Group    19 (            Cl-) has    16 elements
    Group    20 (          Water) has 25581 elements
    Group    21 (            SOL) has 25581 elements
    Group    22 (      non-Water) has  2703 elements
    Group    23 ( Water_and_ions) has 25617 elements
    Group    24 (       nowation) has  2667 elements
    Select a group:

ここの入力でアウトプット対象範囲を指定できる。番号または名前を入力するとその対象が出力される。

ここではまず水分子を含む全体のトラジェクトリのPBCによる分断を補正したいので、`0`または`System`を入力する（どちらでも良い）。

先行入力したい場合には

    input=1lke
    gmx trjconv -f ${input}.trr -o ${input}nj.trr -s ../minimize/min1.tpr -n ../../top/index.ndx -b 0 -e 20000 -skip 100 -pbc whole <<EOF
    0
    EOF

というヒアドキュメント文法が使える。

上記のGroup選択のときには`ndx`形式のインデックスファイルなるもの（デフォルトは`index.ndx`）を自分で作っておいて、`gmx trjconv`時にそのインデックスファイルを指定すれば良い。例えば水・イオン以外の全原子のGroupがほしいときには、このインデックスを作成しておいて、`-n`の引数で読み込ませて使う。`index.ndx`自体は単純なしくみになっているので簡単に後から手動で書き換えることもできると思う。この辺のGroup追加方法はマニュアルを見てください。

この時点で出てきた`1lkenj.trr`のトラジェクトリファイルをVMDやPyMOLなどで表示させてみると、こんな感じ  
![1lkemovie1.gif](https://qiita-image-store.s3.amazonaws.com/0/224327/b308e99e-7092-adf6-195c-b21cc0c53a97.gif)

周期境界ボックスの周辺での水やタンパク質の途切れは修正されているけれど、一緒に結合している基質Digoxigenin（水色）とタンパク質（緑色）が分離しているフレームがいくつか存在している。この補正を2段階目で行う。

### 2段階目の操作

2段階目の操作は、1段階目で出てきた`${input}nj.trr`を入力として利用して

    gmx trjconv -f ${input}nj.trr -o ${input}nj2.trr -s ../minimize/min1.tpr -n ../../top/index.ndx -pbc cluster

となる。すると、今度は1回目の追加入力が

    Select group for clustering
    Group     0 (         System) has 28284 elements
    Group     1 (        Protein) has  2605 elements
    Group     2 (      Protein-H) has  1338 elements
    Group     3 (        C-alpha) has   169 elements
    Group     4 (       Backbone) has   507 elements
    Group     5 (      MainChain) has   677 elements
    Group     6 (   MainChain+Cb) has   830 elements
    Group     7 (    MainChain+H) has   842 elements
    Group     8 (      SideChain) has  1763 elements
    Group     9 (    SideChain-H) has   661 elements
    Group    10 (    Prot-Masses) has  2605 elements
    Group    11 (    non-Protein) has 25679 elements
    Group    12 (          Other) has    62 elements
    Group    13 (            DOG) has    62 elements
    Group    14 (            Na+) has    20 elements
    Group    15 (            Cl-) has    16 elements
    Group    16 (            Ion) has    36 elements
    Group    17 (            DOG) has    62 elements
    Group    18 (            Na+) has    20 elements
    Group    19 (            Cl-) has    16 elements
    Group    20 (          Water) has 25581 elements
    Group    21 (            SOL) has 25581 elements
    Group    22 (      non-Water) has  2703 elements
    Group    23 ( Water_and_ions) has 25617 elements
    Group    24 (       nowation) has  2667 elements
    Select a group: 

となり、クラスタリングして表示したいグループの選択を要求される。ここで、1段階目ではうまくいっていなかった「タンパク質＋Digoxigenin」のグループ、つまり24のnowationを同じ場所にクラスタリングするよう、指定してあげれば良い。

nowation（または24）と入力すると、続いてさらに`Select group for output`と表示され、出力してほしいグループを要求される。ここでは再び水・イオンを含んで全原子の表示を確認したいので、`System`または`0`を選択する。

この出来上がった`1lkenj2.trr`をPyMOLの上で表示するとこんな感じになっている。  
![1lkemovie2.gif](https://qiita-image-store.s3.amazonaws.com/0/224327/b6081875-9fb4-27ee-f525-c3191614a341.gif)
先程はタンパク質とdigoxigeninが分離した状態のフレームが幾つか見られたが、ここでは常にこの2つが結合した状態で表示されているのがわかるだろう。周期境界の溶媒ボックスの外に飛び出してしまっているフレームの表示については次の段階で修正する。

### 3段階目の操作

3段階目は、タンパク質の重心が常に周期境界ボックスの中心にあるように移動させる変換を行う。

    gmx trjconv -f ${input}nj2.trr -o ${input}nj3.trr -s ../minimize/min1.tpr -n ../../top/index.ndx -pbc mol -ur compact -center

例によってここではセンタリング対象とアウトプット対象を尋ねられるので、`nowation`,`System`と順に入力する。変換後はこんな感じ。  
![1lkemovie3.gif](https://qiita-image-store.s3.amazonaws.com/0/224327/985e86fb-c5db-7237-0338-9e35e12598ac.gif)
きれいにボックスの中にセンタリングされたタンパク質+digoxigeninが存在するのがわかる。

### 4段階目の操作

4段階目は、タンパク質の向きを一定に揃えることである。これは`gmx trjconv`のオプションの1つ`-fit rot+trans`で揃えることができるのだが、これまでの`-pbc`オプションとタンパク質の向きを揃える`-fit rot+trans`は同時に使えないらしいので、最後にこれを行う。

    gmx trjconv -f ${input}nj3.trr  -o ${input}nj4.trr -s ../minimize/min1.tpr -n ../../top/index.ndx -center -fit rot+trans

least squares fitの対象, centeringの対象、アウトプットの対象を問われるので、順に`nowation`, `nowation`, `System`の順で入力する。変換後はこんな感じ。  
![1lkemovie4.gif](https://qiita-image-store.s3.amazonaws.com/0/224327/5713db1c-bf7c-cbaa-8bc3-f595de9b6d19.gif)
これで変換完了！！水分子を含むボックスがタンパク質の向きに合わせて回転されているのがわかります。今回のケースではN末端、C末端にやや長いループが存在することもあり、なんかタンパク質の重心がフレームごとに定まっていないように見えるが、ループ構造の少ないタンパク質ならばほとんど動かないように表示されるはず。

## 水・イオンのないトラジェクトリに変換＆出力して表示させたい場合

より現実的な使い方として、シミュレーションの中に非常に多く存在する水分子やイオンを表示させる必要がないケースは多々ある。そんなときは、これまで見てきたアウトプット対象の指定を`System`、つまり全原子から、`nowation`、つまりタンパク質＋digoxigeninだけにすれば、望み通りに水・イオンのないトラジェクトリファイルを作成することができる。
1段階目〜4段階目を水やイオンなしで順にまとめて処理するスクリプトはこちら。

``` bash
#!/bin/sh

input=1lke

gmx trjconv -f ${input}.trr -o ${input}nj.trr -s ../minimize/min1.tpr -n ../../top/index.ndx -b 0 -e 100000 -skip 100 -pbc whole <<EOF || exit 1
nowation
EOF

gmx trjconv -f ${input}nj.trr -o ${input}nj2.trr -s ../minimize/min1.tpr -n ../../top/index.ndx -pbc cluster <<EOF || exit 1
nowation
nowation
EOF

gmx trjconv -f ${input}nj2.trr -o ${input}nj3.trr -s ../minimize/min1.tpr -n ../../top/index.ndx -pbc mol -ur compact -center <<EOF || exit 1
nowation
nowation
EOF

gmx trjconv -f ${input}nj3.trr  -o ${input}nj4.trr -s ../minimize/min1.tpr -n ../../top/index.ndx -center -fit rot+trans <<EOF || exit 1
nowation
nowation
nowation
EOF
```

このような表示になる。  
![1lkemovie5.gif](https://qiita-image-store.s3.amazonaws.com/0/224327/879169cb-f1fa-16f9-498a-93784a77a832.gif)

## AmberToolsを使ってトラジェクトリを変換する方法
**AmberToolsがインストールされており、かつ水・イオンのトラジェクトリを必要としない場合は、AmberToolsの`cpptraj`コマンドを使った方法の方が簡単かつ正確かつ処理も速いです**。

AmberToolsには`cpptraj`というモジュールが同梱されており、主としてトラジェクトリの処理や解析に用います。詳しくはAmberToolsのマニュアルのcpptrajの章を読んでください。
ここで、`trajfix.in`というファイルを用意して中に以下のように書き込んでおきます。

```
# trajinで処理させたいトラジェクトリ名を指定、その後の1 last 10は、1フレーム目から最後のフレーム目まで、10フレームおきに読み込むという意味
trajin 1lke.trr 1 last 10
# 残基番号で1から170に対応する残基を、同じBoundary box内に位置するようにする（ボックスの外に行ってしまっても途中でループさせない）
unwrap :1-170
# 残基番号1-170のCα原子の重心を原点上にセンタリングするよう、各フレームのトラジェクトリを平行移動させる（Gromacsの3段階目の操作に対応）
center :1-170@CA mass origin
autoimage
# 1フレーム目（first）をリファレンスとして、rmsが最小になるよう回転操作を行う（上記4段階目の操作に対応）
rms first out rmsd.dat @CA
# stripコマンドで、SOD, WAT……などの残基に対応するトラジェクトリを出力対象から外す（水・イオンのないトラジェクトリを作成するとき用）。今はコメントアウトしてある。
#strip :SOD,CLA,WAT,TIP3,Na+,Cl-

#以上の処理を施したものを1lkenj4.trrという名前で出力する
trajout 1lkenj4.trr
go
```

この`trajfix.in`ファイルをMDトラジェクトリの`1lke.trr`と同じディレクトリに置き、`cpptraj -p /path/to/leap.pdb -i trajfix.in` と入力すれば（ここで`/path/to/leap.pdb`とは、対象の`leap.pdb`が置かれているディレクトリへのPATHを意味します）、上記の手法より圧倒的に早くかつ正確にトラジェクトリを変換してくれて、上の`1lkenj4.trr`と同じものが作成されます。

上の `gmx trjconv` と比べると、1発で変換でき、処理速度も5倍くらい速く、あまり深く考える必要がなく変換できます。

もし上記の方法でトラジェクトリが不審な動きをする（一部の原子が吹っ飛んでしまうなど）場合は、`trajfix.in`の１行目の`trajin 1lke.trr 1 last 10`の`10`の数字を`1`など小さい値に変えることで解決することがあります。

これでも解決しない場合や、水やイオンを入れたままのトラジェクトリ変換を行いたい場合は、上記の`gmx trjconv`コマンドの方で（時間がかかりますが）解決することがあります。

## （補足）水・イオンのないgroファイルの作成

上で挙げたように、トラジェクトリを載せる時用の`pdb`ファイルや`gro`ファイルを用意する必要があることもあるだろう。
このときの`pdb`ファイルは単純に水やイオンのある行を消去すればよいだけである。
一方で`gro`ファイルについては1つだけ気をつけることがある。
`gro`ファイルの仕様上、1行めはコメント、2行目に**系の原子数**が定義されているため、ここを修正しなければならない。例えば、溶媒和されている`min1.gro`ファイルのコピー`nowat.gro`を作成し、その中身を以下のように書き換える。

    leap
    2667
        1ASP      N    1   2.452   1.725   4.361
        1ASP     H1    2   2.355   1.752   4.347
        1ASP     H2    3   2.457   1.624   4.364
        1ASP     H3    4   2.506   1.758   4.283
    （中略）
      170DOG    H26 2659   2.969   2.969   3.226
      170DOG    H27 2660   2.846   2.911   3.019
      170DOG    H28 2661   2.961   3.008   2.927
      170DOG    H29 2662   2.772   2.974   2.790
      170DOG    H30 2663   2.566   3.081   2.876
      170DOG    H31 2664   2.693   3.211   2.751
      170DOG    H32 2665   2.868   3.193   2.759
      170DOG    H33 2666   2.788   3.391   2.893
      170DOG    H34 2667   2.682   3.286   2.983
       7.14920   7.29594   6.88980

こうすることで、2667原子からなるgroファイルに作り変えることが簡単に可能である。PyMOLやVMDでタンパク質を表示させるとき、必ずトラジェクトリに含まれる原子数と、ベースに載せるタンパク質の原子数は一致していなければならないことを覚えておく。

PyMOLを開く　→　原子数が一致するpdbまたはgroファイルをロードして表示（オブジェクト名をnowatとする） → コマンドラインに`load /path/to/1lkenj4.trr, nowat, state=1`と入力して実行すればうまく表示されるはずである（コマンドラインを使うときは実行ディレクトリやパスに注意）。