- 追加された行はこのように表示されます。
- 削除された行は
このように表示されます。
{{category 系統解析}}
*License
This script is distributed under GNU GPL.
本スクリプトは、GNU GPLに基づいて配布されています。
*決まり文句
例によって例のごとく内容は無保証なので要注意。このページの情報に従って何らかの損害が出ても補償できませんのであしからず。要するに自己の責任において実行してちょうだいねということです。
*概要
近年、系統解析に用いるデータの増加によって複数の遺伝子領域を用いた系統解析が頻繁に行われるようになってきていますが、モデル選択と解析ソフト側の設定が煩雑なため、全領域を連結しただけのデータを用いて塩基置換モデル選択と系統解析が行われていることがしばしばあります。しかし、実際には異なる遺伝子領域は異なる塩基置換の様式を持っていることが多く、それぞれの領域毎にモデル選択を行って選ばれたモデルを適用すべきことが多いと思います。これはそのためのPerlスクリプトです。内部でReadSeqを呼び出してデータ変換を行うことで、ReadSeqで読み込み可能なデータ形式ならほぼ全て扱うことができますので事前のデータ形式変換は必要ありません。とりあえずFASTA・PHYLIP・GenBank・NEXUS形式では動作確認をしています。尤度計算には当初GPLなPHYMLを採用予定でしたが、できたスクリプトで色々なデータを解析すると頻繁に落ちたり算出する尤度が複雑なモデルの方が悪いなど様々な問題があったため、結局PAUP*を使うことになりました。
*注意
これは複数遺伝子領域データを扱う際の一般的な注意点なのですが、各遺伝子領域単独での解析結果が明らかに異なる場合には、組み換えや遺伝子浸透などによりそれぞれの遺伝子が異なる系統に由来する可能性が高いと思います。そのような場合にはどれか、もしくは全体がそもそも系統解析には適していないと考えて下さい。また、単一の領域しか使わない場合は問題無いのかというとそうではなく、単一の領域しか使わない場合には検討のしようがないのでただ分からないだけです。複数領域分のデータがある場合にはこのスクリプトを使って「各遺伝子領域が明らかに異なる系統樹を支持するのか」の検定を補助することは可能です。
*必要なもの
-PAUP*
モデル選択に必要な尤度を計算するために用いる。スクリプトからは「paup」で呼び出せる必要がある(PATHの通っている場所かスクリプトと同じ場所に置いておく)。スクリプトの冒頭で呼び出しコマンドを定義していますのでテキストエディタで開いてその部分を書き換えれば必ずしもPATHの通っている場所に置かなくても構いません。
-ReadSeq
入力データ形式の変換のために用いる。スクリプトからは「readseq」で呼び出せる必要がある(PATHの通っている場所かスクリプトと同じ場所に置いておく)。スクリプトの冒頭で呼び出しコマンドを定義していますのでテキストエディタで開いてその部分を書き換えれば必ずしもPATHの通っている場所に置かなくても構いません。Java版でも使えます。Jarファイルをスクリプトと同じ場所に置いて下さい。
-Perl
スクリプトを実行するために必要。Windowsでは[ActivePerl|http://www.activestate.com/Products/ActivePerl/]でOK。MacOS Xでは標準で含まれているはずです。
*処理内容
+コマンドラインオプションが指定されていれば領域毎に配列を分割する
+タンパクコード領域はコドン位置毎にさらに分割する
+未分割・分割・コドン位置分割のそれぞれでモデル選択
+結果を出力
+MrBayesとPAUP*用のコマンドを出力
MrBayesとPAUP*以外のソフト用の出力は今後実装していくつもりです。また、尤度計算は一度に複数のPAUP*を起動しての疑似マルチスレッド計算に対応していますので、マルチCPUマシンで高速化を図ることが可能です。
*スクリプト
-http://www.fifthdimension.jp/products/kakusan/
からダウンロードして下さい。
*使用方法
使い方は
perl kakusan.pl --help
でメッセージが出ますのでまずは読んで下さい。領域を指定しない場合は
perl kakusan.pl datafile
だけでデフォルトの設定で解析と結果の出力が行われます。
**領域の指定方法
perl kakusan.pl --partition=hoge1:1-300,hoge2:301-600 datafile
みたいな感じで指定して下さい。どこの領域にも含まれない座位が見つかったら止まります。また、複数の領域に重複して含まれる座位が見つかった場合もエラーになります。また、C版のReadSeqを使っている場合は仕様上OTU名にはアンダースコアやハイフンを用いることができません(ReadSeqで正常に変換できない)。どうしてもアンダースコアやハイフンを使う必要がある場合にはJava版のReadSeqを使ってみて下さい。
タンパクコード領域を含む場合は領域名に_Pを付けます。
perl kakusan.pl --partition=hoge1:1-300,protein_P:301-600,hoge2:601-900
このように指定すると、コドン位置毎に分割してさらにモデル選択を行います(分割しない場合も検討します)。ただし、出力オプションにPAUPを指定している場合は、PAUPがProportional modelやSeparate modelに対応していないため、コドン位置毎に分割したモデルは検討しません。
**マルチCPUマシンでの高速化
perl kakusan.pl --numthreads=2 datafile
などとコマンドラインオプションを指定して起動すると、尤度計算時にPAUP*を複数起動して一度に計算することでマルチCPUを生かした高速化が可能です。上記の例では最大起動数を2つまでにしています。Athlon 64 X2やCore Duo、Core 2 DuoなどのデュアルコアCPUでも高速化の恩恵を受けることができます。
*結果を利用した解析方法
実行すると、datafile.kakusanというディレクトリが生成されてその中に色々ファイルができます。MrBayes、PAUP*用の設定ファイルはそれぞれMrBayes、PAUPという名前のディレクトリの中に出力されますのでまずはカレントディレクトリを移動します。
cd datafile.kakusan
cd MrBayes
このディレクトリにはwhole_AIC_hogehoge.nexなどのファイルがあります。配列全体に一つの塩基置換モデルを当てはめる場合(Concatenate model)・各領域毎に異なる塩基置換モデルを当てはめた上で領域間の置換速度比を導入(枝長は領域間で共有)したモデルを当てはめる場合(Proportional model)・各領域毎に異なる塩基置換モデルを当てはめた上で領域間で枝長も共有しないモデルを当てはめる場合(Separate model)のそれぞれのためのファイルがありますので適当に選んで使って下さい。
mb
MrBayes > Execute whole_BIC1_concatenate.nex
AICc/BICの1〜6というのは、サンプルサイズに使う値が異なるため分けています。1が固定した樹形における最小置換数(最節約スコア)、2が各座位の最小置換数の和、3が各座位の形質状態の最小値の和、4が座位数、5が可変座位数、6が座位数×OTU数です。
このスクリプトでは、Concatenate・Proportional・Separate modelの比較は行っていませんのでご注意下さい。Concatenate・Proportional・Separate modelのそれぞれの意味については[Pupko et al. (2002)|http://mbe.oxfordjournals.org/cgi/content/abstract/19/12/2294]のIntroductionをご参照下さい(無料で全文読めます)。ただし、ここで私が言っているConcatenate modelとPupkoらのそれは少し意味が違います。PupkoらのConcatenate modelはProportional modelにおいて領域間で置換速度が等しい場合(つまり塩基置換モデルは領域間で異なる)も含むモデルを指していますが、ここで言っているのは領域間の塩基置換モデルは共通のモデルです。PupkoらがConcatenate modelと呼んでいるモデルはMrBayesではProportional modelの設定ファイルで
PrSet ApplyTo=(all) RatePr=Fixed;
とすることで取り扱うことができます。
現段階では特に出力される設定ファイルのテストが不十分なため、ファイルに問題がある可能性があります。ファイルの内容をよく読んで確認してご使用下さい。
*履歴
-2006/10/29 RC2リリース
-2006/10/26 RC1リリース
*バグ報告・要望
トップページのメールアドレスまでお願いします。
*複数遺伝子領域を用いた系統解析のための塩基置換モデル選択
[分子系統学演習 - データセットの作成から仮説検定まで|http://www.fifthdimension.jp/documents/molphytextbook/]にまとめ直しました。リンク先をご覧下さい。