近年、ベイズ系統推論手法の精緻化が急ピッチで進められており、進化データの共同モデリングのほぼすべての側面において、新たな大きな進歩が見られる。ゲノム配列、サンプリング日、表現型データ、放射性炭素年代、発掘された化石、生物地理学的範囲情報など、複数の独立したデータ源からの証拠を組み合わせることでしか適切に答えられない進化上の問題があることが、ますます理解されるようになっている。関連するすべてのデータを単一のジョイントモデルに含めることは、概念的にも計算的にも非常に困難である。完全なモデル階層に構成できる互換性のある(サブ)モデルのロバストな開発を可能にする高度な計算ソフトウェアパッケージが、こうした開発において重要な役割を果たしている。このようなソフトウェアフレームワークの開発は、それ自体がますます主要な科学的活動となっており、実用的なソフトウェア設計、開発、工学的課題から統計的・概念的モデリングの課題まで、特有の課題を伴っている。BEAST 2はそのような計算ソフトウェアプラットフォームの一つであり、4年以上前に初めて発表された。ここでは、BEAST 2の最初のリリース以来、最近の2.5リリースに至るまでの、BEAST 2のコアプラットフォームとモデル階層における一連の主要な新開発について述べる。
HPより
BEAST 2は、分子配列のベイズ系統解析のためのクロスプラットフォームプログラムです。厳密分子時計モデルまたは緩和分子時計モデルを用いて、rooted(根付き)の時間測定された系統を推定する。系統樹を再構築する方法として使用できるが、単一のツリートポロジーに条件付けすることなく進化仮説を検証するためのフレームワークでもある。BEAST 2はマルコフ連鎖モンテカルロ法(MCMC)を用いて樹木空間を平均化するため、各ツリーはその事後確率に比例して重み付けされる。BEAST 2には、標準的な解析を設定するためのグラフィカル・ユーザー・インターフェースと、解析結果を分析するためのプログラムが含まれている。
BEAST 2はBEAST 1.xを書き直したもので、モジュール性に重点を置いている。これにより、BEAST 2のパッケージシステムによる拡張が容易になっている。現在BEAST 2で利用可能な機能/モデルとBEAST 1で実装されている機能/モデルの比較については、機能表に示されている。
BEAST 2 は非常に大規模で複雑なアプリケーションであるため、初めて使う人にとっては、 プログラムとその機能を理解するのは大変な作業かもしれない。そのため、以下のことを行うことを勧める。
BEAST 2のユーザーや開発者が使用する用語に慣れる。
利用可能なチュートリアル、特にご自分のデータに適用可能なチュートリアルに目を通す。
beast-usersグループにアクセスするか、beast-users+subscribe@googlegroups.com にメールを送って、ユーザーメーリングリストに参加する。このメーリングリストは新バージョンの発表や、バグやその他の問題についてのアドバイスに使われる。
Taming the BEASTワークショップに参加する。
- Model selection(link)
インストール
2つの異なる開発チームのバージョンがある点に注意する*1(どちらもOSS)。
BEAST2系
https://github.com/CompEvol/beast2/releases
https://beast.community/2024-07-09_BEAST_X_released.html
Beast 2のWindowsとMac版をダウンロードしてテストした(ダウンロードリンク)。より古いバージョンもReleaseで公開されている。
#conda install (link:beast2は少し古いバージョン)
mamba install beast -c bioconda -y
mamba install beast2 -c bioconda -y
#launch beauti v1.10
> beauti
#launch beast v1.10
> beast
実行方法
BEASTによる系統推定の手順だけ簡単に確認します。基本的な流れは、
1)入力ファイルの準備(配列データ、あれば地理やコレクション日、信頼できる化石データなど)
2)多重整列
3)BEASTuiでXML作成
4)BEASTのラン
5)(任意) Tracerで確認(Github)
6)Tree Annotatorのラン
7)系統樹ビューアで可視化(Figtree、iTOL、MEGA,TreeViewerなど)
すでに配列の整列まで終わっているとして、ここでは3のBEASTuiからの流れを確認します。
最新版のバージョン2.7.7のwindows版をダウンロードして開いた。BEASTのバージョンに注意する。1.1ならおそらくBEAST X(BEAST 1.1)をダウンロードしている。

BEAUtiとBEASTのバージョンが違うとXMLファイルの互換性がなくなるので注意する(利用可能なモデルの違いなどに起因)。
3)BEASTuiでXML作成
BEAUtiを立ち上げると初期はpartitionsパネルになっている。
Partitionsパネル:パーティションパネルでは、新しいアラインメントやその他のデータ、例えば地理データなどを追加する。
https://beast2-dev.github.io/hmc/hmc//Standard/Partitions/index.html

File => Import Alignmentか左下の+ボタンをクリックしてMSAファイル(NEXUS、BEAST、FASTA形式)を読み込ませる。anolis.nexを使う(アノールトカゲ科アノール属の29種のNADH dehydrogenase遺伝子配列)。ファイルを画面へドラッグしても読み込める。

読み込んだ(パーティション)データを選択後、下のSplitボタンをクリックすると、コドン位置でアラインメントをスプリットすることができる(コーディングヌクレオチドアラインメントを分割し、それぞれに異なるサイトモデルを使用することで、データへの適合を改善する)。下の写真だと、コドンの1番目、2番目、3番目で分割する。

注;partitionとは配列データを分析目的に応じて意味のある部分に分割して扱う事を指す(例;複数の遺伝子、コドン、構造領域)。異なる歴史を持つパーティション(例えば、核とミトコンドリア遺伝子)をリンク(共有)させるとバイアスが起こる可能性があるので注意する。
Splitすると、ツリーは新しいパーティション間で共有される。

分割することで、例えばPartition1にはXXモデルを使い、Partition2にはYYモデルを使うといったことができる。
分割すると、メニューにLinkとUnlinkボタンが出現する。Linkは、ミトコンドリアDNAの異なる遺伝子座や、DNAパーティションと地理データパーティションなど、共通の歴史を共有している場合にのみ、ツリーをリンクさせる。異なる歴史を持つパーティション(例えば、核とミトコンドリア遺伝子)などをリンクさせると、偏った過信的な推定につながる可能性がある(Ogilvie et al.)(マニュアルより)。

LInkは、SiteとClockそれぞれにボタンが用意されている。時計モデルのみ共有して置換モデルは共有しないなど柔軟に設定出来る(異なるツリーを持つパーティションのクロックはリンクしない)。具体的には、パーティションがツリーを共有していて、平均進化速度が異なるが枝ごとの速度変動は同じと考えられる場合(例:ミトコンドリアDNAの異なる遺伝子座)、クロックをリンクさせる。パーティションがツリーを共有していても、枝ごとの速度変動が同じとは期待できない場合(例:遺伝子配列パーティションと地理的位置パーティション)、クロックはリンク解除する。相対速度はSiteパネルの置換速度を用いて推定できる(マニュアルより)。リンクされているかどうかは、パネルの文字が同じかどうかで判断できる。
下の置換ボタン (r) をクリックすると、ファイル選択ダイアログが開いてパーティション内のアライメントの代わりに使用するアライメントを選択できる。また(-)ボタンをクリックするとパーティションデータを削除できる。

パーティションの右端にはAmbiguities(DNAならATGC以外の文字)のチェックボックスがある。あいまいな文字を欠損データとして解釈するか、考慮するかを決める。デフォルトでは欠損データとして扱う(チェック無し)。
データのごく一部のみがあいまいな場合は欠損データとして扱っても問題ないが、多くの部位があいまいな場合は解析結果に影響を及ぼす可能性があるため、あいまいさを考慮する必要がある(マニュアルより)。
Tip datesパネル:個々のサンプルのコレクション日の指定。左上のUse tip datesボタンを有効にすると、taxaごとに過去のコレクション日を指定できる。
https://beast2-dev.github.io/hmc/hmc//Standard/Tip_Dates/index.html

Siteパネル: 置換モデルを選択するパネル。サイトモデルは、アライメントのタイプによって変わる。

ここではパッケージマネージャーでインストールしたOBAMAを指定した(パッケージマネージャーは一番下で説明)。OBAMAは、アミノ酸アラインメントの解析で、モデルの不確実性を考慮して適した置換モデルをベイズモデルで最適なモデルを推定する(Github)。

ヌクレオチドデータにはbModelTest(Bouckaert & Drummond, 2017)などで適したモデルを推定できる(パッケージマネージャーでインストール可能)。
Clockパネル:時計モデルを選択するパネル。
https://beast2-dev.github.io/hmc/hmc//Standard/Clock_Model/index.html

デフォルトではStrict Clock(枝ごとに変えない)モデルが選ばれている。データの種類に関する経験がない初回の解析では、Strict Clockと単純な置換モデル、樹形事前分布を用いて実行することがマニュアルで提案されている。その後、データに十分な数の置換と適切な数の分類群がある場合は、relaxed clock modelを指定する(変動係数の分布がゼロに近い場合は、枝間の速度変動がほとんどないことを示しており、データは厳密時計を支持している。この分布がゼロに触れていない場合は、速度変動があることを示しており、緩和時計が支持される。判断がつかない場合は、パスサンプリングやネストサンプリング解析を実行してベイズ仮説検定を行う(マニュアルより))。
Priorsパネル;ベイズ統計での事前分布(Prior distribution)を指定するパネル。
https://beast2-dev.github.io/hmc/hmc//Standard/Priors/index.html

解析に使うパラメータそれぞれに事前分布を設定する。事前分布に基づいてデータの尤度を計算し、事後分布(Posterior distribution)を得るので、樹形、置換率や時計それぞれについてパラメータの初期想定をする。
(ノード校正が化石に基づく場合は、化石発見過程に関する情報に基づいてノード事前分布を指定できるCladeAgeパッケージ(Matschiner et al., 2017)を検討する。もし化石発見過程が解析内のすべての化石で均一と仮定できるなら、サンプル祖先パッケージ(Gavryushkina et al., 2014)も考慮する(マニュアルより))。
化石から推定可能な祖先ノードの分岐年代情報がある場合、priorsパネルで指定する。
まずAdd Priorを押してMRCA prior(最も近い共通祖先の事前分布)を選ぶ。

そのクレードに含まれるtaxaを選択して右のリストに移し、Taxon set labelに名前をつけたあとOKをクリック。

それから、分布を指定し、Mean、High, Lowなどの値で指定していく。

Nomalだと、Myaなどでノードの年代が正規分布に従ってxx Mya を中心に、yy Mya 程度の範囲を主に取り得ると仮定させる。

(lowerは化石が存在する地層の最も新しい年代、Upperは化石が存在する地層の最も古い年代)
Yule ModelのClockモデルに対しても事前分布の仮定を適用することが可能(未テスト)

MCMCパネル:ベイズ推定の中核であるMCMC(Markov Chain Monte Carlo)アルゴリズムの設定を行うパネル
https://beast2-dev.github.io/hmc/hmc//Standard/MCMC/index.html

MCMCのステップ数(=何回回すか)、何ステップごとに結果の系統樹をファイルに書き出すかなどを指定する。規定単位ごとに系統樹は繰り返し書き出されるので、1000万など大きめの値を指定しておけば、ESS(Effective Sample Size)やトレースの収束具合を確認し、途中で十分な ESS(通常は >200)になっていれば、途中で収束したと判断してそれ以降の出力は無視し、BEAST の出力として使う事もできる(小さい chainLengthだと当然計算は早く終るが、収束しないと無意味なデータになってしまう)。
burnIn(焼き付き:電子機器の品質チェックの用語から来ている。初期テスト・動作確認の目的)は最初の何%のサンプルを無視するかの設定で、初期状態のバイアスを取り除くために指定する。通常は10%〜25%程度を指定する。burinのデータはトレースやツリーログには表示されず、スクリーン出力もない。不要なデータを出さないことで、ログファイルのサイズを小さくし、ディスクスペースとログの転送と処理の時間を削減するのに役立つ(マニュアルより)1000万なら100-250万の値を設定。

最後にXMLファイルを書き出す。v1.1ではgenerateボタンがあったが、2.7.7ではSaveでXMLのパスと名前を指定して保存する。

4)BEASTのラン
Beastを立ち上げる。

保存したXMLを指定する。以前の結果があると書き込み禁止でエラーになるので、その場合はOverwrite: 上書きを選択する。

最後にCPUと命令セット、スレっド数、あるいはGPUを指定し実行する。
ループ回数次第だが、大きなデータではかなりの時間がかかる。

プリントされるのは、MCMC サンプルのiteration数、posteriorがサンプル時点での事後確率(posterior probability)値、likelihoodがデータがモデルとどれだけ一致しているかを示す尤度となる。ループ回数が増えると収束し、負の値の絶対値が小さいほど良くフィットしていることを示している。右端のpriorは事前分布(prior probability)の対数値。収束しない場合は公式の解説を参照。
BEAST 2 の公式ドキュメントでは、デフォルトの 10,000,000 ステップでは不十分な場合があり、データやモデルに応じて調整する必要があると書かれている。指定するチェーン値が分からない場合、まずはデフォルトの Chain Length(10,000,000)で解析を実行し、解析後、Tracer などのツールを使用してESS 値やトレースプロットを確認するのが良いと思われる。ESS 値が十分でない場合やトレースプロットに収束の兆候が見られない場合、Chain Lengthを増やして再解析する。*2, 3
指定したMCMCサンプルのiteration数ごとにプリントされ、指定したiteration数ごとに系統樹ファイルが同じファイルに追記されていく(プリント単位と保存のiteration単位は個別に設定可能)。

繰り返しになるが、十分収束すれば停止してもOK。操作パネルのQuitボタンで終了、Stopボタンで停止できる。ただしStopで停止するとそのまま操作パネルだけ閉じるバグがあり、画面からは直接Resumeできないので注意(JAVAバージョンにより異なる可能性あり)。
XMLのあるパスに出力は書き出される(v1.1では出力パスを指定できたが、2.7.7では出来ない?)。

5)Tracerで確認
Tracerを起動する。TracerはBEAST 2には含まれないのでレポジトリからダウンロードする。
Release Tracer v1.7.2 · beast-dev/tracer · GitHub
Tracerを起動してlogファイルを読み込んだ。

左側のstatisticsから分布を見る統計量を選択する。
このように、Tracerでは、パラメータの推定値:平均、中央値、95% HPD 区間、MCMC の収束状況を視覚的に確認するためのトレースプロット、ESS(有効サンプルサイズ):各パラメータのサンプルの独立性を評価などを確認できる。各統計量でESSが200以上がマニュアルでは推奨されている。
6) Tree Annotatorのラン
TreeAnnotatorのラン。MCMC によってサンプリングされた多数のツリーの中から、代表的な系統樹を1本選び出し、そのツリーに計算した統計情報や進化速度を付与して1本のツリーファイルとして書き出す。
TreeAnnotatorを立ち上げたら、Beastで計算したツリーファイルを指定する。

buninの設定がある(デフォルトは初期10%をpreburninとして扱う)。出力名も指定して実行する。

node heights(縦向きの系統樹でrootからの時間的距離なのでheightと思われる)では、Keep heights を選択した。kepp heightsは選択されたターゲットツリーのノードの高さをそのまま保持する。準備が出来たらRunする。
計算時間は短め

LogCombinerのラン
今回は不要だが、Beast2に含まれるLogCombinerアプリを使うと、複数回のBEAST実行から得られたツリーファイルやログファイルを適切なburninを除去した後に1つのファイルにまとめることができる(TreeAnnotatorと違って複数ツリーからの統計量はファイルにツリーファイルに含まないため、誤差範囲を表示したツリープロットには使えない)。

ファイルを画面にドラッグして指定、出力名を指定する。チェックボックスに1つ以上チェックを付けないと何も実行されない。
7)系統樹ビューアで可視化(FigTreeなど)
最後に系統樹ビューアに読み込んで確認する。
FigTree:FigTreeはBEASTで作成された要約および注釈付きの系統樹を表示することを目的として設計されているので、Figtreeは最も適している(HP)。こちらからダウンロードできる。
https://github.com/rambaut/figtree/releases
Macで起動できない場合など、.jarファイルをターミナル上で指定する(”java -jar Figtree.jar”)。
Icytree:ブラウザ上で使える系統樹ビューア
TreeViewer:MAC, Windows,Linux向けのモダンな系統樹ビューア。表示方法を非常に柔軟に設定することが可能(Pubmed)。
TreeViewerの特徴として、複数のサンプリングされたツリーファイルからコンセンサスツリーやその統計量をアプリ上で計算できる点がある。したがって、Beastからのサンプリングされた多数のツリーファイルを直接読み込んで使用することも出来る。ファイルを読み込み、Module => Transformaer Module => Consensusをクリック、枝長をMeanかMedianで計算してPlotする。burninも指定出来る。

Plot例

(注;time treeはtipの位置が現代で揃う。しかしTreeViewerアプリだとTreeAnnotatorでkeep height設定でベストツリーを選択していてもtip先端がズレて見える可能性がある。=> FIgtreeを使うか、Figtreeでフォーマット変換してから読み込む(ただしnewlick変換だと一部の統計情報が消えてしまう))
その他
BEAUtiアプリのFile => Manage Packagesでは、多くのモデルや解析手法をパッケージとしてインストールできるようになっている。
https://www.beast2.org/managing-packages/index.html
BEASTで使えるMODEL_SELECTIONをインストールした。

インストールしたパッケージを使用するには、一度プログラムを閉じて再度開く必要がある点に注意。
また、BEAUtiのFile => Lauch Appsではアプリを立ち上げられる。ただし開いたアプリを閉じるとBEAUtiも閉じるので注意。

引用
BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis
Remco Bouckaert ,Timothy G. Vaughan,Joëlle Barido-Sottani,Sebastián Duchêne,Mathieu Fourment,Alexandra Gavryushkina,Joseph Heled,Graham Jones,Denise Kühnert,Nicola De Maio,Michael Matschiner,Fábio K. Mendes,Nicola F. Müller,Huw A. Ogilvie,Louis du Plessis,Alex Popinga,Andrew Rambaut,David Rasmussen,Igor Siveroni,Marc A. Suchard,Chieh-Hsi Wu,Dong Xie,Chi Zhang,Tanja Stadler,Alexei J. Drummond
PLoS Comput Biol. 2019 Apr 8;15(4):e1006650.
参考
*1
Why the two different versions of BEAST?
https://groups.google.com/g/beast-users/c/RP-wB_BVxpY?utm_source=chatgpt.com&pli=1
*2
*3
How to beat bad mixing
https://www.beast2.org/2023/11/01/bad-mixing.html?utm_source=chatgpt.com
*4
Figtreeについては、こちらの動画でわかりやすく説明されています。
BEAST cannot find Beagle library
https://github.com/beagle-dev/beagle-lib/issues/146#issuecomment-834010453
https://github.com/beagle-dev/beagle-lib/issues/152
関連
