バイオインフォマティクス体験的入門（分子系統樹の統計推測）
本演習ではバイオインフォマティクスの一例として「分子系統樹の推定」を実際に体験する．公共のデータベースから自分の興味のある動物群（３０種程度）のミトコンドリアDNA配列をダウンロードし，このデータに統計手法を適用して進化の系統樹を推定する．「確率モデルで探る生物の進化」をまず読んで，イメージをつかんで欲しい．
参加者は数人づつのチーム毎に課題Ａ，Ｂのいずれかにとりくむ．チームのメンバーで作業を分担する．さまざまな側面（確率と統計，数値解析，プログラミング，英語読解，ウェブからの情報取得）があるので，それぞれの側面に得意なメンバーがいればスムーズに作業が進むだろう．課題Ｂに関しては難易度が高いので，必ずしもすべてを行う必要はない．出来るところまでを，楽しんでやってもらえればよい．
- 課題Ａ： 系統樹推定のために開発された様々なソフトウエアを利用して系統樹を推定する．最先端の生物学で現在用いられている手法を体験する．具体的には，ＤＮＡ配列のアライメント（位置あわせ）を行うためのclustalw，系統樹推定を行うためのmolphy，得られた系統樹の信頼度を計算するためのconselである．高度な手法を体験するが，すでに開発されているソフトウエアを実行するだけであり，演習としては比較的容易である．（ある学生からの意見： なんだかよく分からないけど，指示通りに手順を実行したら，答えが出ました！）
- 課題Ｂ： 系統樹推定のためのアルゴリズムを学習してソフトウエアを自分で書き，それをＤＮＡ配列データに適用する．データは課題Ａの前半（Ａ１−Ａ４）までを実行して得る．作成するプログラムは，確率過程のモデルからＤＮＡ配列の距離行列を計算するプログラム（Ｂ１）と，それから近隣結合法によって系統樹を再構築するプログラム（Ｂ２），さらにブートストラップ法によって信頼度を計算するプログラム（Ｂ３）である．課題Ａに比べれば若干古い手法であるが，アルゴリズムを自ら実装するため，演習としての難易度は高い．もしすべての実装が困難であれば，一部を課題Ａから流用してもよいから，チャレンジしてみることを薦める．系統樹推定のために開発されたソフトウエアのリストをみれば，この問題がいかに生物学で重要視されているかが分かるだろう．もし自信作が出来れば，その仲間入りができるかもしれない．
アルゴリズムの説明
課題Ａの手順
- 課題の具体的手順は「２００３年度のページ」を参照すること．ただし以下にあげるような変更点がある．
- 「Ａ１．作業ディレクトリの準備」は旧計算機システムを前提に書かれている．２００５年度からの新システムではcmd001は実行しないこと．代わりに次の手順を行う．
- 自分のホームにある ".bash_profile" や".bashrc"に PATH=.:~shimo/class/bin:$PATH という行を加え，ログインしなおす．
- 自分のホームに適当な名前の作業ディレクトリを作成しそこにcdして作業を行う．
- cp ~shimo/class/enshu200209/bin/* . を実行して，ファイル(cmdXXX等）を作業ディレクトリにコピーする．これらのファイルはシェルスクリプトで書かれた実行ファイルである．
- 「Ａ２．進化系統樹」，「Ａ３．データベースへのアクセス」はリンク先の構成が多少変更になっているので注意．データベースのファイルを保存するときは，"save to file"を選択すると良い．またデータベースから得られるデータファイルの形式も若干変更になっているが，分析結果には影響しないはず．
- 「Ａ４．DNA配列のアライメント」以降でcmd101等のシェルスクリプトを実行するときは，sh cmd101[enter]のようにする．もしくは，各ファイルの最初の行にある変更して，#!/usr/local/bin/bash を #!/bin/bash などとしても良い．
- 「Ａ５．系統樹の推定」の最後でcmd371を実行して作成される「work-conc-user-name.eps」は使わないこと．これは与えたリストの１番目のトポロジの図になってしまっている．work-conc-user.mlrとwork-conc-user.treのほうは，正しい結果（尤度を最大にするトポロジ）を与えるので，そのあとの信頼度（確率値）の計算は正しく行える．もしwork-conc-user-name.epsに相当するEPSファイルを作成したい場合は，「TreeView」を用いるか，もしくはwork-conc-user.mlrをみて最尤トポロジがどれかを確認した後で，その系統樹をひとつだけ含むような候補リストをwork-conc-user.tplとして新たに作成し，cmd371を再実行すれば絵がかける．（そうするまえに，最初のcmd371で得られたファイルをすべて保存するなどしておかないと，引き続く確率値の計算がおこなえなくなることには注意してください．）
- 「Ａ６．系統樹の信頼集合」は変更なし．
課題Ｂの手順
- 課題Ａの説明を参照して，「Ａ３．データベースへのアクセス」まで実行すること．これで各生物ごとにNC_XXXXXX.txtという形式のデータファイル（サンプル： NC_001807.txt）の準備ができる．
- 引き続き課題Ａの「Ａ４．DNA配列のアライメント」を実行する．ただし，課題ＡではDNA配列をアミノ酸配列に変換してからアライメントしていたが，課題ＢではDNA配列のままアライメントする．具体的には次の手順を行う．
- 手順1から３までは，課題Ａと同様に行う．これでnsqファイルが作成される（サンプル： work-ND1.nsq)．
- 手順４のcmd141は実行しない．
- 手順５のcmd201の代わりにcmd211を実行する．これでnsqファイルから直接alnファイルが作成される（サンプル： work-ND1.aln）．
- 手順６のcmd221の代わりにcmd231を実行する．これでalnファイルからptnファイルの代わりにnucファイルが作成される（サンプル： work-ND1.nuc）．
- 手順７のcmd241の代わりにcmd251を実行する．これでnucファイルをつなぎ合わせ，work-conc.nucが出力される．このファイルを自分で作成した系統樹推定プログラムに入力する．なお生物種間のDNA配列の違いを分かりやすく見せるために，work-conc-nuc.txtも同時に出力される．
- このようにして準備したＤＮＡ配列データから系統樹推定を行う．「アルゴリズムの説明」を参考にして，自分で以下のＢ１，Ｂ２，Ｂ３に相当するプログラムを作成すること（一般的に入手できるものであれば，Java, C, R等いずれを用いても良い）．すべてを作成することが困難な場合は，部分的に既存ソフトウエアを利用してもよい．また課題Ａで利用するソフトウエアのソースコードを参考にしてもよい．
- プログラムＢ１： アライメントされたＤＮＡ配列データ（サンプル： work-conc.nuc）を入力し，生物種間の距離行列（サンプル： work-conc.dis）を出力するプログラムを作成し，実行する．
- プログラムＢ２： 距離行列（サンプル： work-conc.dis）を入力し，系統樹（サンプル： work-conc-nj.tre, work-conc-nj.tpl）を出力するプログラムを作成し，実行する．系統樹を図（サンプル： work-conc-nj.pdf, work-conc-nj-name.pdf）にするには「TreeView」等を用いるか，自分でプログラムを作成し描画する．
- 上記Ｂ１，Ｂ２に関しては，課題Ａの「Ａ５．系統樹の推定」を参考にしてよい．具体的には次の手順を行う．
- 手順１は実行しない．（すでにcmd251によってwork-conc-nuc.txtが出力されている．）
- 手順２のcmd311の代わりにcmd321を実行すると，距離行列と系統樹が出力される（Ｂ１，Ｂ２のサンプルファイル）．
- 手順３のcmd331の代わりにcmd341を実行すると，近隣結合法とは別の推定方式である最尤法によって系統樹が推定される（サンプル： work-conc.mlr, work-conc-ml.tpl, work-conc-ml.pdf, work-conc-ml-name.pdf）．この出力には，系統樹の各枝に信頼度が記入されている．これは以下のＢ３で計算するものと類似のため参考のためにあげてあるが，実際にこの手順３で用いた方法（ＬＢＰ）は，Ｂ３で用いる方法とは異なるので注意．
- プログラムＢ３： ブートストラップ法によって系統樹の信頼度を計算する．各枝の信頼度を記入した図（サンプル： work-conc-ml-name.pdf）を作成する．
- 上記Ｂ３に関しては，基本的にはデータに乱数を適用してブートストラップデータを多数（１００個から１００００個程度）生成し，各データにＢ１＋Ｂ２を繰り返し適用する．各データからひとつ系統樹が得られるので，それを並べたファイルを作成する（サンプル： work-conc-user.tpl ただしこれはブートストラップで得られるものではなく，課題Ａ５の手順５で得られる105個の系統樹のものである．あくまでファイル形式の参考のためにあげてあるだけ）．これをconselパッケージのtreeassプログラムに入力すれば，各系統樹と各枝（クラスタ）の頻度が計算される．treeassはＵＮＩＸでいえばuniq+sortに相当するプログラムで，系統樹を並べたファイルを読み込んで，そのなかで同じ系統樹や同じ枝（クラスタ）が何回出てきたかを数える．もちろん，この部分も自分たちで作成すればなおよい．