意識低い系の系統解析・Modeltest-NGとRAxML-NGで系統解析 発動篇

説明するより読んでもらった方が早いようなこと。

詳細なModeltest-NGの使い方はgithubのWikiを、RAxML-NGの使い方もやはりgithubのWikiなりpdfのマニュアルを観れば良い。 ここでは解析の流れをだけ示し、とりあえず動かしてみる事を目的とする。 なんとなく使い方を把握したら公式の情報を読んで各オプションについて理解を深め、使いこなせる様になって欲しい。

201811にやったもの。

1. 単独の遺伝子座を用い、パーティションを指定しない場合

まずは簡単な解析の事例を示す。

準備

てきとうなアライメントファイルを作ってFastaかPhylipで書き出す。

Modeltest-NGでの解析

とりあえずフォルダを作成し、中にアライメントファイルを配置する。 ここではphilip形式の「10OTU.phy」とするが、Fastaファイルも解析可能。 Terminalで作成したフォルダにcdで移動、↓のコマンドを実行する。

modeltest-ng -i 10OTU.phy -o testestes

補足:
-i アライメントファイルの指定。
-o 出力ファイル名の指定。 無指定の場合はアライメントファイルと同じ名前が付く。
-p マルチスレッドの指定。 今回は指定していないがCPUに合わせて指定することで計算が速くなる。 「-p 4」等、物理コア数に合わせるのが良いと思われるのでシステムレポートなどで確認のこと。

コマンドを実行すると、いくつかのファイルがフォルダ内に作られる。「.out」が付くファイルがモデルテストの結果で最重要! ここでは出力ファイル名の指定をしているので「testestes.out」として作られる。

testestes.outをテキストエディタで開くと3種の情報基準BIC、AIC、AICcに基づくモデルテストの結果とコマンド例が記されている。 どれに基づくかはお好みで。 注目するのはコマンド例。

Commands:
> phyml -i 10OTU.phy -m 010020 -f m -v 0 -a e -c 4 -o tlr
> raxmlHPC-SSE3 -s 10OTU.phy -m GTRGAMMAX -n EXEC_NAME -p PARSIMONY_SEED
> raxml-ng --msa 10OTU.phy --model TrN+G4
> paup -s 10OTU.phy
> iqtree -s 10OTU.phy -m TrN+G4

↑の様に、各種ソフトウェア向けのコマンド例があるが、続く解析でRAxML-NGを使うので「raxml-ng --msa 10OTU.phy --model TrN+G4」に注目。 この例では進化モデル「TrN+G4」が選ばれていることがわかる。

RAxML-NGでの解析

次にRAxML-NGでの解析に移る。 コマンドを自分の目的に合わせて編集する。

raxml-ng --all --msa 10OTU.phy --model TrN+G4 --prefix testesphy -bs-trees 1000

補足:
--all ブートストラップ値と最尤系統樹の推定を行う。
--msa マルチプルシークエンスアライメントファイル、要はアライメントファイルの指定。 Modeltest-NGで使ったものを指定。
--model 進化モデルの指定。 上で出力されたものを入力。
--prefix 出力ファイル名の接頭語指定。
-bs-trees ブートストラップ試行数指定。
--threads マルチスレッドの指定。 Modeltest-NGの「-p」と同様。

最初に作った、配列ファイルがあるフォルダをTerminalで開き、上のコマンドを実行する。

フォルダ内にいくつかのファイルが作られる。 その中の「testesphy.raxml.support」にブートストラップ値つき最尤系統樹が記されている。

樹形とブートストラップ値の確認

「testesphy.raxml.support」をFigTreeで開き、Branch Labelsにチェックを入れ、Displayオプションでlabelを選択して表示させる。

これだけの作業でとりあえずブートストラップ値を含む最尤系統樹が作成できる。

2. 複数の遺伝子座を用いる場合、あるいはパーティションを用いる場合

複数の遺伝子座や、タンパク質コード領域をコドン別に用いたい場合、配列をパーティション分けして解析する必要がある。 Modeltest-NGでは、事前にパーティションファイルを用意しておく事でそういった解析を行うことができる。

準備

アライメントファイルseq.txtとパーティションファイルpartition.txtを使うものとする。 アライメントファイルはFastaまたはPhylip形式。 どこからどこまでがそれぞれの遺伝子なのかメモっておく。 パーティションファイルはテキストエディタで↓の形式に編集する。 例ではCOI領域の各コドンをパーティションとして分けている。 RAxMLと同じような形式っぽい。 2つのファイルをフォルダにまとめておく。

DNA, COI_1 = 1-1551\3
DNA, COI_2 = 2-1551\3
DNA, COI_3 = 3-1551\3

コドン分けではなく、単に複数遺伝子座を連結しただけなら↓の様な感じ。 当然、複数遺伝子座かつコドン分け、というのも可能。

DNA, tekitou = 1-1551
DNA, sorenari = 1551-2500
DNA, nantoka = 2501-4000

Modeltest-NGでの解析

↓の様なコマンドをまとめ、Terminalでアライメントファイルとパーティションファイルをおいたフォルダに移動し、実行する。

modeltest-ng -i seq.txt -q partition.txt -t ml

補足:
-i アライメントファイルの指定。
-q パーティションファイルの指定。 テキストエディタで作成。
-o 出力ファイル名の指定。 無指定の場合はアライメントファイルと同じ名前が付く。
-t 開始樹形の指定。 デフォルトはmpらしいので今回はMLにしてみた。
-p マルチスレッドの指定。 今回は指定していないが「-p 4」とか指定できると速くなる。
-T 今回は不要だが、解析ソフトに合わせて選択されるモデルを制限することができる。 MrBayesあたりを使う時に必要。

パーティション指定した場合出力されるファイルが増える。 「.out」が付くファイルが結果である事に代わりはないが、続く解析の為のファイルが別途出力される。 出力ファイルのうち、seq.txt.part.bic、seq.txt.part.aic、seq.txt.part.aiccが系統解析用のモデルファイル。 bic,aic,aiccに基づくものがそれぞれ出力される。 

出力されたモデルファイルはRAxML-NGではそのまま使えるが、MrBayesだとコマンドブロックを手動で編集する必要がある様だ。

今回はmodeltest-ngが出力したBICに基づくモデル選択ファイルseq.txt.part.bicを使用するものとする。 ファイルの中身は↓の様になっており、各パーティションに対して選択された進化モデルが記されている。

TIM2ef+I+G4, COI_1 = 1-1551\3
F81+G4, COI_2 = 2-1551\3
TIM3+I+G4, COI_3 = 3-1551\3

RAxML-NGでの解析

次にRAxML-NGでの解析に移る。 コマンドを自分の目的に合わせて編集する。

raxml-ng --all --msa seq.txt --model seq.txt.part.bic --prefix testrun2 --bs-trees 1000

補足:
--all ブートストラップ値と最尤系統樹の推定
--msa マルチプルシークエンスアライメントファイル、要は配列ファイルの指定。 Modelstest-NGでの解析に使ったファイルを指定する。
--model 進化モデルの選択、進化モデルを記したパーティションファイルの指定。Modeltest-NGが吐いたファイルをそのまま指定できる。
--prefix 出力ファイル名の頭につく文字列の指定。 整理の得意な人は指定しなくても困らないかもしれない。
--bs-trees ブートストラップの計算回数指定。 
--threads マルチスレッドの指定。 Modeltest-NGの「-p」と同様。

TerminalでModeltest-NGの解析を行ったフォルダに移り、コマンドを実行する。 しばらく、あるいは長時間待つとフォルダ内に解析結果ファイルが出力される。

出力ファイルの内の「testrun2.raxml.support」がブートストラップ値つきでスコアの良い最尤系統樹のデータ。 とりあえずこれが重要。

樹形とブートストラップ値の確認

「testrun2.raxml.support」をFigTreeで開き、Branch Labelsにチェックを入れ、Displayオプションでlabelを選択して表示させる。

これだけの作業でとりあえずブートストラップ値を含む最尤系統樹が作成できる。