概要

vmdで水素結合の解析をする。GUIとCUIでの方法を紹介する。

環境

macOS Mojave 10.14.5
vmd 1.9.3

GUIでの方法

まずは構造ファイルとトラジェクトリファイルを読み込む。読み込み方はこちら → VMDの使い方の基本 (Mac)

水素結合の数をプロット

VMD Main ==> Extensions ==> Analysis ==> Hydrogen Bonds
以下を参考に詳細を設定する。
- Selection 1 : protein
- Selection 2 : resname G39
- Update selections every frame? : チェックを外すと計算時間が短くなる。
- Only polar atoms (N, O, S, F)? : チェックする
- Calculate detailed info for : Residue pairsをチェックする
- Output directory : お好み
- Frame/bond data? : お好み
- Detailed hbond data? : お好み
- Find hydrogen bonds! : 設定が終わればここをクリックし、計算を始める
プロットが現れる。FileからExport to PostScriptを選択すると画像を保存できる。また、Export to ASCII matrix...を選択すると計算結果データを保存できる。

アウトプットファイルの説明

アウトプットファイルは2つある。Frame/bond data?で指定された結果ファイル(test.dat)とDetailed hbond data?で指定された詳細ファイル(test-detailed.dat)が出力される。

test.dat : 計算結果が書かれている。各フレームの水素結合の本数が記述されている。

#frame number
0 2
1 0
2 0

test-detailed.dat : 計算結果の詳細が書かれている。見つかった水素結合の組み合わせと、出現頻度が書かれている。Calculate detailed info forオプションの設定によって詳細ファイルの内容が変わる。以下、出現頻度について補足を加える。

Residue pairs open 下の例では、ARG290-SideとG39389-Sideでの水素結合出現頻度が0.30%となっている。今回は全部で1000フレームあったので1000*0.30/100=3フレームで水素結合があったことになる。

Found 16 hbonds.
donor        acceptor    occupancy
ARG290-Side      G39389-Side     0.30%
ARG212-Side      G39389-Side     0.10%

Unique hbondopen 下の例では、残基名は同じ（ARG290-SideとG39389）だが、NH1、NH2、O1B、O1Aという原子名の違いで出現頻度が計算されている。

ARG290-Side-NH1      G39389-Side-O1B     78.22%
ARG290-Side-NH2      G39389-Side-O1A     76.22%
ARG290-Side-NH1      G39389-Side-O1A     0.10%

All hbonds open 下の例では、ARG290-SideとG39389-Side間での水素結合出現頻度が154.55%となっており、100%を超えている。この理由は上記Unique hbondの例からわかる通り、78.22+76.22+0.10=154.54%となる。よってAll hbondsでは出現頻度は原子間で計算しているが、組み合わせが残基間で表示されていることがわかる。

ARG37-Side   G39389-Side     0.10%
ARG290-Side      G39389-Side     154.55%

Hydrogen Bondsメニューオプション

Input options
- Molecule：解析したい分子の指定
- Selection 1：解析したい原子の指定
- Selection 2：解析したい原子の指定
- Frame：解析したいフレームの指定
- Update selections every frame?：フレームごとにSelection 1,2で指定した原子の選択を毎回更新するかどうかの指定。原子を時間依存しないように指定（ex 原子名など）していればチェックを外すことで計算時間を節約できる。原子を時間依存するように指定（ex 座標での指定、またはprotein within 5 of (resname G39)のように、時間経過によって選択される原子が変化する）していれば毎回更新しなければならない。
- Only polar atoms：解析対象原子をN, O, S, Fに限定する
- Selection 1 is the：Selection 1で選んだ原子をドナー、アクセプター、またはその両方として解析するかどうか。Selection 1でドナーとすると、Selection 2で選んだ原子はアクセプターとなる。
- Donor-Acceptor distance：水素結合を定義するときの距離
- Angle cutoff：水素結合を定義するときの角度
- Calculate detailed info for：解析結果の詳細情報を定義する。ここが一番理解しにくい。具体例はtest-detailed.dat。
  - None：何もしない
  - All hbonds：原子間で水素結合があると１とカウントする。しかし、集計自体は残基間で行われる。よって、出現頻度が100%を超えるものもある。
  - Residue pairs：残基間で水素結合があると１とカウントする。残基間で水素結合が2本以上あったとしても、１とカウントする。
  - Unique hbond：原子間で水素結合があると１とカウントする。
Output options
- Output directory：出力ファイルのディレクトリ
- Log file：ログファイル名
- Frame/bond data：出力ファイル名
- Detailed hbond data：解析結果の詳細ファイル名

CUIでの方法

vmdのスクリプト作成が初めてという方はこちら→VMDのスクリプトの書き方。

measure hbonds

水素結合情報はmeasure hbondsで取得できる。このコマンドで取得できるのは水素結合を形成しているドナー、アクセプター、水素のインデックスを要素とするリスト。リストの長さは３つ（ドナー、アクセプター、水素）なので、ここから水素結合の本数、角度、距離など様々な情報を取得できる。ここでは、例としてGUIで出力した結果ファイル（test.dat）と同様の出力をする。

mol load pdb test.pdb dcd test.dcd

set out_file [open test.dat w]
set num_frames [molinfo top get numframes]

for {set i 1} {$i < $num_frames} {incr i} {
    set donor [atomselect top "protein" frame $i]
    set acceptor [atomselect top "resname G39" frame $i]
    set hbond_li [measure hbonds 3 20 $donor $acceptor]
    set num_hbonds [llength [lindex $hbond_li 0]]
    puts $out_file "$i\t$num_hbonds"
}
exit