生物情報工学・実習

分子動力学法と構造ゲノム科学について

00/06/23

  1. UNIX環境へのログイン
  2. 分子動力学法について
  3. 分子運動の可視的な観測
  4. 構造ゲノム科学
  5. 最後に
  6. 付録

UNIX環境へのログイン



  1. はじめに


    前回の実習と同様に、今回の実習でも教育用計算機システムを利用して実習を行います。今回は分子動力学法の復習をした後、実際にプログラムをコンパイルし実行してみます。そのために教育用計算機システムに用意されているUNIXの環境を使用します。 (教育用計算機システムのWindows NTにはプログラムを開発する環境がありません。)
     

  2. UNIXサーバへのログイン


    今回の実習ではUNIXの環境を利用します。 NCデスクトップ画面で左側にあるアイコンメニューの中の「Unix」アイコンを左クリックすると、 UNIXのログインウィンドウが現れます。NCを使い始めるときにユーザー名とパスワードを入力して、ユーザー認証を行いましたが、UNIXを使うときは更に別の計算機資源を利用するので、再度ユーザー名とパスワードを入力してください。

    ログインに成功するとUNIX画面(正確にはXウィンドウ)が現れます。

分子動力学法



  1. 分子動力学法とは


    液体や溶液のように多数の原子や分子が動き回る系において、それらの粒子の振る舞いを観測する手法として分子動力学法(Molecular Dynamics = MD)があります。
    MDは個々の粒子にNewtonの運動方程式を適用してそれを積分することで、個々の粒子座標の時間経過の観測や系全体のエネルギーの算出ができます。
    しかし非常に多くの粒子が相互作用しあっているので、膨大な計算が必要になります。このために計算機を利用します。

    今回の実習では、アルゴン原子の動きを観測しその系のエネルギーを算出するために、実際にプログラムを使用して実行結果を考察します。

    実際にプログラムを使用する前に、MDのアルゴリズムの復習をします。
     

  2. Verlet法による運動方程式の解法


    N個の分子からなる系において、2分子間のポテンシャルがφ(r)で与えられたとすると、系全体のポテンシャルエネルギーΦは2分子間のポテンシャルの和になります。

        -(1)

    ここでのrijは分子iとjの距離でri=(xi,yi,zi)と表すと、

        -(2)

    となります。
    こうして与えられたポテンシャルエネルギーから、分子iの座標でこのポテンシャルエネルギーを微分することで、分子iに働く力Fi=(Fxi,Fyi,Fzi)が分かり、例えばFxiについて次のようになります。

         -(3)

    ｙ、ｚ成分においても同様です。

    次に運動方程式を近似的に解くために、ごく短い時間刻みΔtごとに運動方程式の解を求めることにします。
    時刻tにおける分子iの座標をri(t)とすると、時刻t±Δtにおける座標はri(t±Δt)となります。これを時刻tのまわりでテーラー展開を2次まですると次のようになります。

         -(4)

    この両辺の和と差をとると

         -(5)
         -(6)

    ここで時刻tでの分子iの速度vi(t)と加速度ai(t)は

          -(7)
         -(8)

    とあたえられます。よって(6)と(7)、(5)と(8)から

         -(9)
         -(10)

    という、速度と座標がtの漸化式として得られます。ここでは速度は座標を計算した後に計算します。

    またこの式は初期条件ri(0+Δt)では使えないので、次の式を利用します。

         -(11)

    このようにVerlet法によって導かれた漸化式は、運動方程式の近似解を計算機で求める際に都合がよい式になります。

    また一般に分子間のポテンシャルを表す関数として、レナード-ジョーンズポテンシャルが用いられます。

          -(12)

    以上のことを踏まえて、今回の実習ではC言語のプログラムを利用します。


  3. ソースファイルのダウンロード


    実習時間中に１からＣ言語のプログラミングをしてもらう時間がないので、用意したソースファイル(プログラミングのコードが書いてあるファイル)を使います。

    下の２つのファイルをダウンロードしてください。

    md.c - MDによる分子シュミレーションのプログラム
    input.dat - アルゴン原子の初期座標と初速度のデータファイル



  4. 分子動力学法の実行


    プログラムを実際に計算機上で走らせるために、まずコンパイルします。
    次のようにTerminalに入力してください。

    % cc md.c -o md -lm

    これでコンパイルをする事ができました。次に今コンパイルしたプログラムを実行します。

    % ./md
    ITERATE STEP 50
    ..
    ITERATE STEP 100
    ..
    ITERATE STEP 150
    ..
    ..
    .

    実行開始後しばらくの間標準出力にステップの途中経過が表示されます。
    計算が終わるとTerminalにプロンプトが戻り、４つのファイル「output.dat」「kinetic.dat」「energy.dat」「temperature.dat」ができているはずです。

    % ls
    energy.dat    input.dat    kinetic.dat    output.dat    temperature.dat    ....    .....    ...
    ....    ...
    ...
     

    と、入力する事で確認する事ができます。


  5. 計算結果のプロット


    計算結果のファイルの中身はすべて数値なので、プロットしてグラフとして計算結果を考察します。

    教育用計算機システムのUNIXにはプロット用のツールとして「gnuplot」という UNIXの標準的なプロットツールがインストールされているのでこれを利用してプロットします。
    次のようにTerminalに入力してください。

    % gnuplot

            G N U P L O T
            unix version 3.5
            ..    ..    ..    ..
            ,,    ,,
            ..

            ..    ..    ..
            Send bugs, suggestions and mods to
    gnuplot >

    gnuplotがこれで起動しました。出力ファイルをプロットするには次のように入力してください。

    gnuplot > plot "FILE NAME" with linespoints

    'FILE NAME'には、ここでは「kinetic.dat」、「energy.dat」、「temperature.dat」のどれか１つ見たいものを入力してください。
    すると新しいwindowが開いて、そこに計算結果がプロットされます。更に複数のグラフを同時に比較したい時は、コンマで区切って並べます。

    gnuplot > plot 'FILE NAME 1' with linespoints, 'FILE NAME 2' with linespoints

    3つのグラフを比較すると次のことが分かります。



分子運動の可視的な観測


    MDの計算によってごく短い時間後との分子の座標も計算結果として得られる。この結果を３次元処理をして、分子の動きを可視的に観測することができます。

    先ほど計算した、アルゴン原子の動きを見てみます。
    画面の左端にあるデスクトップメニューからWindows NTを起動して下さい。 (注：Unixシステムから起動したNetscapeから実行する事はできません。)

    Windows NTのNetscapeから

    アルゴン原子の分子運動

    にアクセスして下さい。

    Windows NTはのちほど使うのでログアウトしないで下さい。画面右上の最小化のボタンを押すことでWindows NTの画面を一時的に見えなくすることができます。

    データのダウンロードが終わると、再生画面があらわれます。再生ボタンをクリックすると動画が動きます。

    それ以外にもいくつかの動画が用意してありますので見て下さい。

    tRNA分子の基準振動モード
    (糖-リン酸骨格の二面角およびグリコシル結合の二面角をパラメータとして計算。)
    rasp21の分子動力学計算結果
    (水をいれて計算していますが、蛋白質だけ表示しています。)
    へリックス構造を持つペプチドのMDによる計算結果
    (Poisson-Boltzmann を組み込んだ計算)



構造ゲノム科学 (Structural Genomics)


  1. 構造ゲノム科学とは

    ゲノム時代…ヒトをはじめとして、全ゲノムの塩基配列を解読しようとするプロジェクトが世界各地で進められています。すでにさまざまな生物種で全ゲノムが決定されています。例: 微生物、 ヒトなどの日本でのプロジェクト情報

    遺伝子の塩基配列が決定されると、それぞれの遺伝子の機能を知ることが重要になります。塩基配列だけがわかっていて、機能が未知の遺伝子の機能を予測するためにはどのような方法をとればよいでしょうか。(下をクリックする前に少し自分で考えてみてください。)

    遺伝子の機能を理解するためには、アミノ酸配列からもう一歩進めて、その遺伝子がタンパク質として発現したときの立体構造を知ることも重要です。このように、構造を手がかりに遺伝子の機能を理解しようとする試みは「構造ゲノム科学」と呼ばれています。 タンパク質のアミノ酸配列だけでなく、立体構造がわかると、何ができるようになるかを考えてみましょう。

    配列は似ていないけれど構造と機能が似ている、というタンパク質の例をPDBで検索して、Rasmolを用いて表示してみましょう。
    まず、「レグヘモグロビン」を検索します。

    手順：

    1. PDBのホームページにアクセスする。
    2. 右の「Search」の中の「SearchLite」をクリックする。
    3. 中央のテキストボックスに「Leghemoglobin」と入力して「Search」ボタンを押す。
    4. 一番上の「1BIN」の欄の右端の「{EXPLORE}」をクリックする。1BINに関する情報が表示される。
    5. 左の「Structural Neighbors」をクリックする。
    6. 様々な構造類似性検索ソフトが表示される。上から２番目の「CE」をクリックする。
    7. 検索に用いるパラメータ群を指定するウィンドウが表示される。「Specify Selection Criteria:」の中の「Length Difference:」を「<30%」に、「Sort by:」を「RMSD」に、それぞれ変更する。
    8. 「Search Database」ボタンを押す。
    9. 1BINのA鎖に類似した構造をもつタンパク質の一覧が、RMSD(各原子が平均何オングストロームずれているか)の順に表示される。その中からミオグロビン「1EBC:A」(上から２０数個目あたりにある)を見つけ、左端のボタンをクリックして選択状態(ボタンが押し込まれた状態)にする。配列相同性は12.3%と低いが(配列は異なる)、RMSDは2.6オングストロームと小さい(構造は似ている)ことを確認する。
    10. ウィンドウの上方の「GET ALIGNMENT」ボタンをクリックする。
    11. 配列アラインメント結果が表示される。配列相同性が低いことを実際の配列の上で確認すること。
    12. 下方の「View Results」の中の一番上「Download alignment as a PDB file」をShiftキーを押しながら左クリックし、「1BINa_1EBCa.txt」という名前をつけて保存する。
    13. そのままではRasmolで２つ同時に表示することができないので、kterm上で以下のコマンドを実行する。

      % sed -e "s/ENDMDL/TER/" 1BINa_1EBCa.txt > 1BINa_1EBCa.pdb

    14. これで２つのタンパク質の立体構造を同時に見ることができるPDBファイルを作成することができました。のちほどNT環境の方で前回使用したRasMolを使って見ることにします。
    15. 表示する際にはRasmolで1BINa_1EBCa.pdb を表示させる。1BINと1EBCに別の色をつけて表示するときは、「Colours」メニューの「Chain」を選択。さらに「Display」メニューの「Ribbons」を選択すると、見やすくなる。

    1BIN(レグヘモグロビン)と1EBC(ミオグロビン)は、配列相同性が低いのに、構造が非常によく似ていることを確かめて下さい。いずれもヘムと呼ばれる鉄錯体を結合して酸素運搬・貯蔵に関わるなど、機能的にも類似しています。

    ※PDBにうまくアクセスできないときは、このページからPDBファイルを取得してください。
    ※UNIX上で作成したファイルをWindows NTで使うためには、Windows NTの「スタート」から「Unix Home Directory」を開くとUNIXのホームディレクトリを見ることができます。そこからWindows NTのホームディレクトリにドラッグ＆ドロップすることで移すことができます。

  2. タンパク質の立体構造を予測する

    タンパク質の立体構造を実験的に求める手法には、 が広く行われていますが、結晶化が困難なタンパク質やNMRで原子座標の決定が困難なタンパク質も存在します。このようなタンパク質に対しては、コンピュータを用いた立体構造予測が有効になります。

    背景：

    代表的な構造予測手法と特徴：

    また、タンパク質の立体構造をもとに、複数のタンパク質がどのような配置で結合するかを予測する手法(ドッキング予測)も、さまざまなものが考えられており、Web上でも公開されています。

    構造予測コンテスト(Critical Assessment of Techniques for Protein Structure Prediction, CASP)と呼ばれる会議が、２年に１度開催されています。これは、立体構造がまだ公表されていないタンパク質のアミノ酸配列だけが先に公表され、各グループがそのアミノ酸配列をもとに予測構造をsubmitし、各グループのsubmitが終了したあとに正解の立体構造が公表されて、どのグループの手法が成績がよかったか、などが議論されるというものです。2000年に行われたCASP4では、ab initio法の台頭が注目されました。

    (余力があれば、以下にチャレンジしてみてください)

    CASP4において、ab initio法でよい成績をおさめたBakerらのグループのターゲットNo.102の予測結果とターゲットの正解の構造をダウンロードし、どれだけ似ているかをRasmolで確認せよ。

    ダウンロード手順はこちら。

最後に



    以上で、実習を終わります。
    今日の実習の出欠の確認のために前回と同様に

    report@bi.a.u-tokyo.ac.jp

    へ、自分の名前、所属と実習の感想を一言書いたメールを送ってください。

付録：ソースコードについて


    授業で扱ったプログラムについて解説します。実際のプログラムと照らし合わせながら見て下さい。

  1. データファイルの定義

    分子の初期状態がかかれているデータファイルである「input.dat」と、計算結果を出力するためのファイル「output.dat」「kinetic.dat」「energy.dat」「temperature.dat」を定義しています。

  2. 分子座標データの読み込み>

    「input.dat」にかかれている座標データを、このプログラムで割り当てたメモリにロードして、座標データをプログラムで利用できるようにします。

  3. 相互作用する相手の分子の決定

    分子の空間を無限に設定すると、計算機で扱える程度の分子の数ではエネルギーが発散してしまいます。
    逆に空間を制限すると分子と壁との相互作用が大きくなってしまいます。
    そのために周期境界条件を用います。これはある一定の基本空間は設定するものの、その境界は物理的な壁ではなく分子は移動できます。しかし系全体は基本空間の連続したものとして計算するとこで、分子数の減少やエネルギーの発散も起きないという仕組みです。
    つまりここでは周期境界条件が適用されてるので、相互作用の相手の分子が基本空間から飛び出ている分子に対して、その飛び出た分子が周期境界条件を適用することで、基本空間中のどの分子かを指示します。(図を参照)
     

  4. ポテンシャルエネルギーと分子に働く力の算出

    このプログラムでは「レナード-ジョーンズ ポテンシャル」と呼ばれる最も広く使われているポテンシャル関数型を利用しています。
    また分子に働く力はポテンシャルをその分子の座標で微分することにより得られます。
    (式(12)と式(3)を参照)

  5. 加速度の算出

    すでに分子に働く力がわかっているので、運動方程式から加速度が求まります。
    (F=ma ⇔ a=F/m)

  6. Verlet法による分子座標と速度の更新

    運動方程式をVerlet法で近似して解くので、ごく短い時間であるΔtだけ時間を進めてそのときの分子の速度を計算して、そこからその時の分子の座標を求めることを繰り返します。
    (式(9)と式(10)と式(11))

  7. 運動エネルギーの算出

    運動エネルギーは
        -(13)
    で与えられるので。この式を計算します。
  8. 分子座標の更新

    Verlet法によって求めた分子の座標を周期境界条件と照らし合わせて、基本空間内における分子座標を決定します。

  9. 計算結果の書き出し

    「kinetic.dat」に運動エネルギー、「energy.dat」に全エネルギーの総和、「temperature.dat」に温度、「output.dat」に座標を書き出しています。
TOPページ
生物情報工学研究室