全ゲノム配列を使って系統樹を「素早く」「簡単に」作りたい。

概要

系統関係の推定にシングルコピーオルソログを用いることは信頼性の高い系統樹を得るための一つの手段です。これまではシングルコピーの遺伝子を得るために、各サンプルのゲノム配列の用意ーアノテーションの実行ー得られたタンパク質の配列の相動性検索が必要で、系統解析を始める前の準備に手間と時間がかかっていました。

ここでは、これらの課題を解決するための CUSCO パイプラインについて、特徴と最小限の手順で使い始める方法を紹介します。インストールから基本的な実行例までを通して、バイオインフォマティクス初学者にとってもすぐ試せる内容になっています。これから系統解析を始めたい方や、既存の方法に煩雑さを感じている方にとって、本記事が第一歩の助けとなれば幸いです。

特徴

「簡単」シングルコピーの取得から系統樹構築までをコマンド1行で実行してくれます。
「早い」シングルコピーオルソログ取得を従来の相動性検索から、アノテーションに切り替えたワークフローを実装しました。
「汎用性」ゲノム配列のfastaさえあれば、データベースから取ってきたものでも実行できます（アノテーションされてなくてもOK）。倍数体モードを実装したことで、倍数体を含んだ系統関係の推定も可能になりました。
「お手軽」手持ちのデータセットに新しいサンプルを追加して系統関係の推定をしたい時に、種の系統関係を反映するシングルコピーオルソログ（マーカー遺伝子）をリストアップしてくれます。サンプルのゲノムをわざわざ読まなくても、サンガーシークエンスでその遺伝子の配列だけ分かれば、全ゲノムを使った系統樹と同じ樹形の遺伝子系統樹が得られます。

インストール(25/08/20時点)

各ソフトのバージョン管理などの解析環境を自動で調整してくれるminiconda (anaconda) を使ってインストールします。minicondaのインストールは各自でお願いします。

順番としては、①新しいconda環境を作る→②Polyphestをインストール→③CUSCO環境のインストール→④Astral, pangeneをインストールしてすべてのソフトにパスを通す、となります。やってることは、うまくcondaパッケージに入らなかったソフト３つを手動で入れているだけです。

CUSCOのgithubページ（ https://github.com/seikot345/CUSCO ）に飛んで、Astralのリンクを押し、インストールされたAstral.5.7.8.zipだけカレントディレクトリにおいてください。以下の15行を1行づつポチポチすればインストールは完了です。

Install

# Create new environment
conda create -n cusco python==3.10.16 pip
conda activate cusco

# Install polyphest
git clone https://github.com/NakhlehLab/Polyphest.git
cd Polyphest
pip install -r requirements.txt

# update environment
cd ../
git clone https://github.com/seikot345/CUSCO.git
cd CUSCO
conda env update --file environment.yml
mv cusco.py ../
mv multree_builder.py utils.py ../Polyphest/polyphest
cd ../

# Install pangene
curl -L https://github.com/lh3/pangene/releases/download/v1.1/pangene-1.1-bin.tar.bz2|tar jxf -

# Unzip Astral
unzip Astral.5.7.8.zip

# Add PATH of both softwares
export PATH=$PATH:$(pwd)/pangene-1.1-bin/bin_x64-linux:$(pwd)/pangene-1.1-bin/scripts:$(pwd)/Astral:$(pwd)/Polyphest:$(pwd)/Polyphest/polyphest

※ 途中 multree_builder.py と utils.py をmvしてますが、これは投稿時点でPolyphest公式のスクリプト中のコマンドが古いため、Polyphestのエラーを防ぐために修正したスクリプトを置き換えています。

実行

カレントディレクトリが次のようになってたら準備OKです。

./

Astral  cusco.py  pangene-1.1-bin  Polyphest    

inputディレクトリとoutputディレクトリを作り、さらにinputの中にspeciesディレクトリを作ります。

mkdir input
mkdir output
cd input
mkdir species

inputディレクトリの中にinputファイルを用意します。CUSCOに必要なものは系統樹を作りたいデータセットのゲノム配列のfastaとその内１種のタンパク質データ（アミノ酸配列）のfastaです。今年はヘビ年なのでコブラ科のヘビのデータを使った例を見せます。

NCBIのAssemmblyデータベースにセットし、Elapidae（コブラ科）と入力し、適当にアセンブルファイルをダウンロードします。

今回は、Hydrophis major, Hydrophis melanocephalus, Laticauda colubrina, Naja naja, Notechis scutatus の５種のゲノムfastaをダウンロードして ./input/species に移します。また Notechis scutatus のタンパク質fastaもダウンロードして ./input に移します。すべてが終わるとinputディレクトリ内は次のようになります。

./input/

protein.faa  species

./input/species/

GCA_004320005.1_hydMel_1.0_genomic.fna  GCA_009733165.1_Nana_v5_genomic　　GCA_015471245.1_latCor_2.0_genomic.fna  GCA_033807585.1_HMAJ_1.0_genomic.fna  GCA_900518725.1_TS10Xv2-PRI_genomic.fna

分かりやすいインストラクションのため、ゲノムファイルの名前も変えたいですが、種数が多いと面倒くさいですね。
そこでCUSCOのedid機能を使いましょう。

cd ../
python cusco.py edid -species

これで ./input/species/ 内のファイルが次のようになりました。

./input/species/

Hmajor.fna  Hmelanocephalus.fna  Lcolubrina.fna  Nnaja.fna  Nscutatus.fna

この直したファイル名でシングルコピーオルソログの対応表などが出力されるので、自分の分かりやすいファイル名に変更することをお勧めします。

ゲノムファイル、リファレンスタンパク質ファイルが準備できたら、次のコマンドを実行します。

実行

python cusco.py -t 8 -r proteins.faa

-tで使用するthread数を指定しないと使用可能なCPUを全て使用するので注意。PCのスペックにもよりますが、実行時間は１日あれば終わっていると思います。実行終了後にはinputディレクトリにいろいろファイルができますが、今回は説明を省きます。outputディレクトリをみると次の４つが出力されています。

実行結果

marker.list  single_copy_genes  single_copy_genes.tsv  species_tree.tree

・single_copy_genesディレクトリにはシングルコピーオルソログのfastaが出力されます。今回のデータセットでは13427個のシングルコピーオルソログが検出されたかと思います。ファイル名はリファレンスのproteins.faaのIDと一致します。
・single_copy_genes.tsvにはリファレンスと対応する各サンプルの遺伝子IDが表になっています。各サンプルの遺伝子IDは./input/fasta/にあるfastaファイルに由来します。
・species_tree.treeはnewick形式の系統樹ファイルです。このファイルが全シングルコピーオルソログを使って構築された系統樹になります。Figtreeなどのフリーソフトを使って描画してみてください。
・marker.listについては後述します。

マーカー遺伝子

CUSCOの機能の一つにマーカー遺伝子の抽出があります。ここで言うマーカー遺伝子とは”その遺伝子の系統樹のトポロジーが全ゲノムを使って構築した系統樹と同じである遺伝子”を指します。この抽出には、遺伝子系統樹と種の系統樹間の Robinson-Foulds distance (ファイル中のnRF) を指標にしています。CUSCOの実行が終わったら、outputディレクトリに marker.list というファイルが出来上がります。ファイルの中身はシングルコピーオルソログとその nRF がタブ区切り形式で出力されています。nRFは０～１の間を取り、０に近いほど比較した二つの系統樹のトポロジーが同じ（０で完全一致）になります。

marker.list

gene    nRF
XP_026532018.1  0.00
XP_026533831.1  0.00
XP_026537500.1  0.00
XP_026534022.1  0.00
XP_026534834.1  0.00
...

このリストが活躍するのは次のような状況です。

”自分はゲノムの読まれていないコブラ科ヘビのDNAサンプルを持っている。そのサンプルを６種目として、いま作った系統樹に含めて系統関係が知りたい。だがしかし、ゲノム解読のための研究費なんてない。。”

なんてあきらめてた人、朗報です。このマーカー遺伝子で系統樹を構築すると、全ゲノムデータで系統樹を描いたときと同じのトポロジーの系統樹になる可能性の高い、つまり、ゲノムを読まなくても信頼性の高い系統関係の推定が可能になるのです！

などと説明しましたが、実際は何個の遺伝子を追加で読めばいいのか、どの遺伝子がベストなマーカーなのか、など検証しきれていない部分は多数あり、必ずしも”正解”の系統樹を持つ遺伝子をリストアップしているわけではないこと、注意してください。詳しい検証と議論は元論文を参照いただければと思います。あくまで全ゲノムシークエンスが想像を絶するくらい低コストになった世界線を迎えるまでの、ニッチな需要に答えた機能、だと認識してもらえたらと思います。

もしマーカー遺伝子戦略を実行する場合、マーカー遺伝子用のプライマーを設計して、サンガーシークエンスで読むことになると思いますが、そのプライマー設計をサポートする機能もCUSCOには実装しています。マーカー遺伝子としてひとつの遺伝子と２～３適当な種を指定してcusco.py primeを実行してください（今回はmarker.listの一番上の遺伝子と H. major, N. naja の2種）。

python cusco.py prime -g XP_026532018.1 -s Hmajor Nnaja

すると XP_026532018.1_primerdesign.fasta というfastaが ./output/primer_design に出力されます(下図)。①はマーカー遺伝子の各exon、②は-sで指定した２種の①を包括したゲノム領域を示してくれます。

倍数体モード

CUSCOの隠れた長所。CUSCOは倍数体を含むデータセットでも実行可能です。ここでは６倍体の小麦を例にCUSCOを実行してみたいと思います。
今回は、Aegilops speltoides, Aegilops tauschii, Triticum aestivum, Triticum　monococcu, Triticum urartu の５種のゲノムとTriticum urartu のタンパク質を使用します。

倍数体モードの場合は ./input/ ディレクトリ内に次のようなcopy_number.tsvを用意する必要があります。

copy_number.tsv

Aspeltoides	1
Atauschii	1
Taestivum	3
Tmonococcum	1
Turartu	1

小麦（T. aestivum）は Aegilops speltoides(２倍体), Aegilops tauschii(２倍体), Triticum urartu(２倍体) から成る異質６倍体です。そのためデータセット内では他の種に比べ３倍のサブゲノムごとのシングルコピーオルソログが存在すると考えられます。倍数体モードにおいてcopy_number.tsvで各サンプルごとに保有するであろう遺伝子数を指定することで、データセットごとに柔軟な”シングルコピーオルソログの抽出”が可能になります。

実行は引数 -p をつけるだけです。

（倍数体モード）実行

python cusco.py -t 8 -r proteins.faa -p

実行後の ./output/ ディレクトリ内には次のような出力があります。

（倍数体モード）実行結果

marker.list  single_copy_genes  single_copy_genes.tsv  species_tree.enewick  species_tree.tree

シングルコピーオルソログの数は9789個になったかと思います。ここではあまり見かけない二つの系統樹ファイルについて説明します。

１つ目がspecies_tree.enewickです。このファイルを見てみると、

species_tree.enewick

(Aspeltoides,(((Tmonococcum,(Turartu,(Taestivum)#H0)),(Atauschii,#H0)),#H0));

というよく見るtreeファイルと少し違うことが分かります。このファイルは Extended Newick と言われる形式でtreeの一部にネットワークが含まれる形式です。このファイルを視覚化したい場合、Dendroscopeを使うことをお勧めします。実際に開いてみるとこんな感じです。

T. aestivum のサブゲノムの由来が示せました。

２つ目がspecies_tree.treeです。このファイルを見てみると、

species_tree.tree

(Aspeltoides,(Taestivum,((Tmonococcum,(Taestivum,Turartu)),(Atauschii,Taestivum))));

となります。普通のtreeファイルじゃないかって？よく見てください。copy_number.tsvで3と指定した T. aestivum が３回登場しています。この Multi-label tree は Extended Newick と同様に、倍数体のサブゲノムがどの系統に由来するものなのかが示されています。

終わりに

以上、CUSCOの簡単な使い方を紹介しました。「あれ、思ったより簡単じゃん」と感じてもらえたでしょうか。
まだまだ発展途上のツールなので、「もっとこういう機能が欲しい」「こう使うと便利だった」みたいなアイデアがあれば、ぜひコメントやGitHubのIssueで教えてください。系統解析という生物学研究の第一歩を効率化する一助になれば嬉しいです。