2022/09/13 追記、誤字修正
2024/01/16 タイトル修正、わかりにくい説明を修正
Current Protocols in BioinformaticsにBUSCOの使い方に関する論文が出ています。（引用1）。その論文のパートの１つ；”Support Protocol 3: BUILDING PHYLOGENOMIC TREES”では、BUSCOがシングルコピーの系統マーカー遺伝子を取り出すことを利用して、ゲノム配列から遺伝子セットを取し出し、複数の遺伝子から系統推定するプロトコルについて説明されています。BUSCOの標準のコマンドではサポートしていない方法のため、この論文では、snakemakeベースでいくつかの外部ツールとスクリプトを組み合わせてこれを行っています。手順を確認してみます。
インストール
依存
- snakemake
- conda
この系統推定パイプラインでは、busco_protocol/support_protocol3/envsにある３つのenvbiroment.yamlの環境がアクティベートされる（リンク）。これらの環境のツールは、初回ランで自動導入される。
git clone https://gitlab.com/ezlab/busco_protocol
cd busco_protocol/support_protocol3/
Snakefileとconfig.yamlが確認できる。
config.yamlファイル
config.yamlにはgenomeのfastaファイルのパス、GFF形式のアノテーションファイルのパス、出力、BUSCO5.2.2に対応したBUSCOデータセットなどを書く。３つの環境のYAMLファイルへの相対パスも書かれている。これらの環境は、BUSCOのラン、Rのスクリプトのラン、コア遺伝子の抽出とMSA＆トリミングなどで使い分けられる（リンク）。
genomesがgenomeのfastaファイルを含むディレクトリへのパス、GFFsがGFF形式のアノテーションファイルを含むディレクトリへのパス、resultsが出力ディレクトリ、geno_suff: ".fna" ならfastaファイルの拡張子は.fna、insecta_odb10がBUSCOデータセット（参考；v4リンク）。BUSCOデータセットは使うゲノムの分類によって変える必要がある。分類がよく分からないゲノムなら、マーカー遺伝子数は少ないがfungi_odb10のような大きな分類を使う手もある。percent_missing: 0だとすべてのゲノムに存在。keep_duplicates: "No"ならduplicates BUSCOに相当するマーカー遺伝子は使われない（デフォルトでNo）。
準備
上ではgenomeとなっているので、genome/に系統推定したいゲノムのfastaファイル（.fna）、GFFs/にアノテーションのgffファイルを配置する（対応するfastaとgffのファイル名は拡張子以外一致している必要がある）。
さらに、カレントにrulesディレクトリ、scriptsディレクトリ、envsディレクトリ、ルールの中にもenvディレクトリを配置する必要がある。本テストデータではこのままで問題ない。
Snakefileファイル
configfileではカレントパスにあるconfig.yamlが指定されている。Snakefileの中でincludeして呼び出されるSnakefileファイルは相対パスでrules下に含まれている。
最後に呼び出されてiqtreeを実行するSnakefile；rules/run_iqtree.smkが以下。
iqtree2のランで、アウトグループを-o ”Zootermopsis_nevadensis”としているが、これは論文中でgenomes/に配置されたfastaファイルの1つ；Zootermopsis_nevadensis.fnaをアウトグループにしているため。自分のデータでは修正するか、このオプション自体を消す必要がある。他のパラメータについても同様にチェックしておく。
ラン
準備ができたら実行する。30スレッド指定。マーカー遺伝子のトリミングやMSAで標準エラー出力が膨大に出てくるので、ファイルに保存するか捨てる。
snakemake --cores 30 --use-conda 2> error.log
それぞれのマーカータンパク質配列はdata/に、連結タンパク質のMSAや系統推定結果はresults/直下に保存される。
results/
phylogeny.treefileが最終出力。６つ見えるサブディレクトリはゲノム毎のBUSCO5.2.2出力。
（論文で使用されている６つの昆虫ゲノムの結果）。
data/extracted/にはゲノム毎にディレクトリができる。中には、それぞれのマーカータンパク質配列が保存されている。
smakemakeで実行されるフローはシンプルなため、解析手順に興味がある方は、rulesに含まれている８つのファイルを開いて確認してください。修正したいパラメータ設定があればこのファイルを書き換えるだけです。簡単ですね（その場合、コマンドだけでなく、どこを変えたかのも記録も残しておきましょう）。
この論文はオープンアクセスです。詳しくは論文を読んで下さい。
2024/01/15追記 ~tidyverseの導入失敗に伴うエラー~
SP3.yamlからcondaの仮想環境が作られるが、SP3-2.yamlのR3.6のtidyverseの導入に失敗するため、ラン途中でエラーが発生するようになってしまっている。対策として、SP3-2.yamlを開いてRのバージョンを3.6から4.2に書き換えて実行したがそれでもsnakemakeによる自動インストールで失敗していた。そこで環境をconda activateでアクティベート後（コマンド実行時のエラーログに長い数値からなる環境名をアクティベートしたというメッセージが残っている。おそらく最初に発生したエラーの少し前のlogあたり）、"mamba install r::r-tidyverse -y"でtidyverseを導入してdeactivate後、それから実行することで回避した。
引用
1
BUSCO: Assessing Genomic Data Quality and Beyond
Mosè Manni, Matthew R Berkeley, Mathieu Seppey, Evgeny M Zdobnov
Curr Protoc. 2021 Dec;1(12):e323
関連
参考
