9
6

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

IQ-TREEで系統樹の最尤推定

Last updated at Posted at 2021-09-02

塩基配列やアミノ酸配列を用いた最尤法による樹形探索プログラムはいくつもあるが、とりわけ高機能で速いのがIQ-TREEである。ここでは、1遺伝子座を用いた推定および複数遺伝子座を用いた推定の2パターンを例にIQ-TREEの使い方を説明する。なお、公式のマニュアル(英語)が超充実しているので、基本的にはそちらを参考にしたほうがいい。

IQ-TREEを使うメリットは以下の4つである。

  • 類似プログラムと比較して探索能力と効率が高い
  • 置換モデル選択も同時にやってくれる
  • 各枝の信頼度をUltrafast BootstrapやSH-aLRT testなど、高速なアルゴリズムで評価できる
  • Webサーバーでも使える

間違えがあればコメント等でご指摘下さい

環境

macOS 11.4
IQ-TREE2 2.0.7

インストール

Homebrewを使っている人はbrew install iqtree2でインストール可能(バイオインフォ系のパッケージをHomebrewに加えていない人はbrew tap brewsci/bioを先に実行しよう)。

直接ダウンロードする場合は公式サイトから。
IQ-TREE Webサーバーを利用する場合はインストール不要。128コア、メモリ1TBのLinux環境で解析を回せるぞ!

IQ-TREEはバージョン1と2が提供されているが、基本的にはIQ-TREE 2の方がオススメ。
Webサーバー版では現在バージョン1.6.12が提供されている。

1遺伝子を用いた樹形推定

インプットデータ

Phylip, Fasta, Nexus, Clustal形式に対応している。
アライメントは事前にMAFTT等で済ませておくこと。

解析

mtDNAの1遺伝子や部分配列、パーティションを設定しない連結配列を用いた場合は以下のようにします。

ターミナル
iqtree2 -s input-data.fasta -m MFP -bb 1000 -alrt 1000 -nt AUTO

# -s でインプットデータを指定
# -m で置換モデルを指定。MFP(ModelFinder Plus)を指定するとModelFinderを使ったモデル選択を実施
# -b は標準的なBootstrap法の繰り返し数を指定。100以上を指定すること。遅いので基本使わない
# -bb はUltrafast Bootstrap法の繰り返し数を指定。1000以上を指定すること。速い
# -alrt はSH-aLRT法の繰り返し数を指定。各枝の信頼度を評価する。速い
# -nt でスレッド数を指定。AUTOで自動設定

-m MFPを指定したとき、ModelFinderによってベイズ情報量規準(BIC)を基準に置換モデルが選択される。赤池情報量規準(AIC)と補正赤池情報量規準(AICc)の計算も行われ、こちらを採用したい場合は-AIC-AICcオプションを付ける。モデル選択のオプションについては公式HP参照。置換モデルの一覧はこちら(公式HP)

高速な枝の信頼度の評価は上記のUFbootとSH-aLRT以外にも、近似ベイズ法-abayesを指定)によるものと、高速ローカル・ブートストラップ法(-lbp 1000を指定)によっても可能である。これらの手法はすべて両立可能で、系統樹上でそれぞれの支持率(近似ベイズでは事後確率)が/で区切られて表示される。

Ultrafast Bootstrapの信頼度は一般的なブートストラップの値とは解釈が異なり、95%以上のNodeのみ信頼できるといえる。SH補正近似尤度比検定(SH-aLRT)の値は80%以上の場合に信頼度が高いと言える。UFbootとSH-aLRTの両方を算出することがオススメされている。

複数遺伝子座を用いた推定

インプットデータ

配列ファイルはPhylip, Fasta, Nexus, Clustal形式に対応している。
アライメントは事前にMAFTT等で済ませておくこと。

配列は連結(concatenate)していても、領域(パーティション)ごとに個別のファイルになっていてもOK。
ただし、1サンプルごとにファイル分けをしていると解析はできないので、ひとつのファイルに解析したい個体をすべてまとめる必要がある。

パーティションファイルはRAxML形式とNEXUS形式に対応している。公式は柔軟な設定が可能なNEXUS形式をオススメしている。

NEXUS形式のパーティションファイル(モデルを指定する場合)

Runする際に-mは使わない。

partition.nexus
#nexus
begin sets;
    charset COI = 1-658;
    charset cytb = 659-1700;
    charpartition mine = HKY+G:COI, GTR+G:cytb;
end;

パーティションファイルでモデルを指定しない場合は charpartition mine = の行を消して、-m オプションで特定のモデルもしくはモデル選択を指定する。

NEXUS形式のパーティションファイル(複数の配列ファイルを使う場合)

配列ファイルが複数に分割されている場合、同じディレクトリに配列ファイルを置いた上で、パーティションファイルにはcharset 領域名=ファイル名: 位置と記述する。

partition.nexus
#nexus
begin sets;
    charset COI = PCG.phy: 1-658;
    charset cytb = PCG.phy: 659-1700;
    charset 12S = rRNA.phy: 1-1000;
    charset 16S = rRNA.phy: 1001-2600;
    charpartition mine = HKY+G:COI, GTR+G:cytb, JC69:12S, TN93:16S;
end;

配列ファイル内にパーティションが含まれない(=ファイル内のすべての配列をそのまま使う)場合は、位置を指定する代わりにカンマで区切り、*を指定する。
charset D-loop = Control-region.phy: *;

NEXUS形式のパーティションファイル(コドン位置ごとに解析する場合)

1遺伝子座のみの解析の場合でも、コドン位置ごとに分割して解析する場合は位置の後に\3と書いたパーティションファイルを準備する。

partition.nexus
#nexus
begin sets;
    charset COI-1st = 1-657\3;   #1stコドン(1, 4, 7...655塩基)のみの解析
    charset COI-2nd = 2-657\3;   #2ndコドン(2, 5, 8...656塩基)のみの解析
    charset COI-3rd = 3-657\3;   #3rdコドン(3, 6, 9...657塩基)のみの解析
    charpartition mine = HKY+G:COI-1st, GTR+G:COI-2nd, GTR+G:COI-3rd;
end;
RAxML形式のパーティションファイル(モデルを指定しない場合)

すべての遺伝子座に-mで指定した同じモデルが適用される。ただし、-m MFPとした際は領域ごとにベストフィットモデルが適用される。

partition.raxml.txt
DNA, COI = 1-658
DNA, cytb = 659-1700
RAxML形式のパーティションファイル(モデルを指定する場合)
partition.raxml.txt
GTR+I, COI = 1-658
HKY+I+G4, cytb = 659-1700

マニュアルにはRAxML形式ではすべての領域に同じ置換モデルが適用されると書いてあったが、上記のような記述をすると領域ごとに異なるモデルが適用された。正しい挙動をするか不明なので、おとなしくNEXUS形式で書くことをオススメする。

解析

パーティションを設定した複数遺伝子座を用いた解析は以下のようにします。

ターミナル
iqtree2 -s input-data.fasta -p partition.nexus -bb 1000 -alrt 1000 -nt AUTO

# -s でインプットデータを指定。配列ファイルが複数の場合は指定しない
# -p でパーティションファイルを指定
# -m はパーティションファイルで置換モデルを指定していない場合のみ使う

その他の解析

いくつかは時間ができたら詳しく解説したい。とりあえず、羅列しておく。

SNPデータの解析

SNPデータを扱う場合、-m GTR+ASCのように置換モデルにascertainment bias correctionを含めて解析する。モデル選択する場合は-m MFP+ASCとする。なお、ミッシングが多すぎると解析できない。

Time treeの推定

バージョン2.0.3からサンプルや祖先分岐の時間情報を用いて、最小二乗法による分岐年代推定をする機能が実装された。化石較正を使った古い分岐年代の推定というよりはウイルスの分岐のような短いタイムスケールの解析を想定しているようにみえる。詳細は公式のマニュアルを参照のこと。

Concordance factorsの計算

複数の遺伝子座を使って樹形推定をしたとき、すべての領域を連結した時の系統樹と各遺伝子の系統樹および各サイトの系統樹がどれだけ一致しているかを調べる手法。遺伝子のconcordance factorsをgCF、サイトのconcordance factorsをsCFと呼ぶ。コンカテ樹の推定それぞれの遺伝子樹の推定gCFとsCFの算出という3ステップで解析する。詳細は公式のマニュアルRob Lanfear博士の解説を参照のこと。

アウトプットファイル

解析が終わると、ログの最後の方に以下のメッセージが表示され、いくつかのファイルが書き出される。

ターミナル
...
Analysis results written to: 
  IQ-TREE report:                example.phy.iqtree
  Maximum-likelihood tree:       example.phy.treefile
  Likelihood distances:          example.phy.mldist

Ultrafast bootstrap approximation results written to:
  Split support values:          example.phy.splits.nex
  Consensus tree:                example.phy.contree
  Screen log file:               example.phy.log

Date and Time: Thu Sep  2 21:18:30 2021


% ls #ディレクトリ内のファイル一覧

example.phy		example.phy.iqtree	example.phy.splits.nex
example.phy.bionj	example.phy.log		example.phy.treefile
example.phy.ckp.gz	example.phy.mldist
example.phy.contree	example.phy.model.gz

このうち、filename.iqtreeには結果の概要やコンセンサス樹が記述されている。filename.treefileはnewick形式の系統樹で、系統樹を可視化する際はこのファイルをFigtreeなどで開く。

出力ファイルに特定の接頭辞を付けたい場合は、解析時に--prefixで指定する。上書きを防ぐ時に便利。

参考文献

9
6
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
9
6

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?