#2022/7/1 追記：hintsファイルによるサポートのないアノテーションを除去
前回の記事 BRAKER2のインストール によってBRAKER2を実行する環境ができていることを前提とする。特に GeneMark 関連の設定ミスが大抵のエラーの原因なので、ハマッてしまったと思ったらよく見直そう。
インプットは RepeatModeler/RepeatMasker で Softmasking を実施した GENOME.softmasked.fasta と、OrthoDB からダウンロード・マージした proteins.fasta だ。
コマンドはシンプル
conda activate braker
braker.pl --genome=GENOME.softmasked.fasta --prot_seq=proteins.fasta --softmasking --cores=24
--cores は並列計算で使うCPUの数。途中のプロセスで最大8コアしか使わないところがあるが、結局沢山使った方が早く終る。が、48を超えると挙動が不安定になりがちなようだ。
サイズが 1 Gbp 以下の比較的小さいゲノムなら２～３日待てば終わる。結果はカレントディレクトリの中に braker というディレクトリが生成されてそこに書き出される。最終産物は braker.gtf。
余計なアノテーションを除去
braker の遺伝子予測はかなり優秀で、brakerだけでだいたい遺伝子の95%くらいは予測できてしまう。ただ問題は、予測される遺伝子数がべらぼうに多いことだ。2倍体の植物ではだいたい3万遺伝子前後が妥当な数字だが、braker で予測される遺伝子数はたいてい5万前後になる。ORFが組めるというだけで余計な遺伝子まで予測してしまうからだ。そして、gtf ファイルの中身をを見てもどれが Protein.fasta をもとに予測した遺伝子で、どれが純粋に予測しただけの遺伝子なのか区別がつかない。
だが、braker では解析の途中で Protein.fasta のマッピングで遺伝子の位置を予測したhintsfile.gff というものが生成されるので、このヒントとオーバーラップしてる遺伝子は、少なくとも既存の遺伝子コード配列に基づいた予測遺伝子ということになる。
やり方だが、https://github.com/Gaius-Augustus/BRAKER/tree/report/scripts/predictionAnalysis のページにある、predictionAnalysis.py と selectSupportSubsets.py を使う。
ということで、↓のコマンドを braker ディレクトリの中で実行して、スクリプトをダウンロードし、実行権限も付与する。
wget https://raw.githubusercontent.com/Gaius-Augustus/BRAKER/report/scripts/predictionAnalysis/predictionAnalysis.py
wget https://raw.githubusercontent.com/Gaius-Augustus/BRAKER/report/scripts/predictionAnalysis/selectSupportedSubsets.py
chmod +x predictionAnalysis.py selectSupportedSubsets.py
次に、↓を打つ。
./selectSupportedSubsets.py --fullSupport braker.FULL.gtf --anySupport braker.ANY.gtf --noSupport braker.NO.gtf braker.gtf hintsfile.gff
これで、hintsファイルと完全一致している予測遺伝子は braker.FULL.gtf に、少しでも hints とオーバーラップがある予測遺伝子は braker.ANY.gtf に、hints に基づかない、純粋な予測は braker.NO.gtf に分類して吐き出される。当たり前だが、ANY には FULL も含まれるよ。fully-supported と partially-supported を敢えて分けたい場合は自分で差分を取ろう。