さて、BRAKER2 による遺伝子アノテーションと、GeMoMa による近縁種のアノテーション移植が終わったら、それらを合わせしまう。ただ、単純に両方を合わせるだけだと同一ローカスに異なる CDS を予測してしまうことも多いので、基本的には GeMoMa のアノテーションを採用し、GeMoMa では全く遺伝子が予測されなかった領域に BRAKER が予測した遺伝子があればそれを採用する、という形にしたい。もちろんこの辺はもう各ユーザーの好みになってくる話だが。
何にせよ、複数のアノテーションツールでそれぞれ得られた GFF ファイルや GTF ファイルから差分を取りたいときに有用なのが gffcompare (https://ccb.jhu.edu/software/stringtie/gffcompare.shtml) である。2つの GFF/GTF ファイルの被り具合を、各遺伝子座について評価してくれるので非常にわかりやすい。詳しくは上記リンク先の Transcript classification codes を参照。分かりやすいイメージが使われているので直感的に理解できるだろう。

1. インストール

mamba install gffcompare

ついでにもう一つの便利ツールである gffread もインストールしておくとよいかも知れない。詳しくは kazumax　先生の解説ページで。

2. 実行

GeMoMa でできた final_annotation.gff と BRAKER でできた braker.gtf をインプットに使う。GeMoMa の OUTPUT フォルダに braker.gtf をコピーしておくと操作が楽。

gffcompare braker.gtf -r final_annotation.gff -o gffcmp

gff/gtf ファイルはオプションなしのスペース区切りでいくつでも指定可能。
-r では基準とする gff/gtf を指定。
-o でアウトプットファイルの接頭文字 (prefix) を指定。

結果としていくつかのファイルが出力されるが、必要なのは gffcmp.braker.gtf.tmap だ。タブ区切りテクストなのでエクセルから開くこともできる。

3 列目の class_code というところが２つのアノテーションの被り具合を表している。それぞれの遺伝子ID 対応関係も含めたテーブルになっているので大変扱いやすい。

3. 結果のフィルタリング

class_code が "u" と "i" と "y" あたりに限定するのが妥当だと思う。保守的に "u" だけでも十分かも知れない。

cat gffcmp.braker.gtf.tmap | awk '$3=="u" || $3=="i" || $3=="y" { print $4 }' > braker.keep.list

とすれば、BRAKER でだけ予測された遺伝子をリスト化できる。もちろん grep -w や cut -f を使っても構わない。list ができたら、braker.gtf から list に載った遺伝子だけをフィルタリングする。

cat braker.gtf | grep -wf braker.keep.list > braker.keep.gtf

4. GTF ファイルを GeMoMa 用に変換する

次に、final_annotation.gff に　braker.keep.gtf をあらためて GeMoMa に放り込むのだが、GeMoMa は各遺伝子座ごとに "gene" というアノテーションが付いた gff ファイルしか受け付けてくれない。だが安心しよう。そのためのコマンドは GeMoMa が用意してくれている。ワーキングディレクトリは前回と同じ /mnt/c/アカウント名/GEMOMA/ を想定。

conda activate gemoma
GeMoMa AnnotationEvidence a=OUTPUT/braker.keep.gtf g=GENOME.fasta outdir=OUTPUT

結果として、annotation_with_attributes.gff というファイルができる。これと、final_annotatin.gff を一緒に Annotation Fileter に掛ける。

5. GeMoMa Annotation Filter

GeMoMa GAF g=OUTPUT/final_annotation.gff g=OUTPUT/annotation_with_attributes.gff tf=true outdir=OUTPUT

filtered_predictions_1.gff というファイルができる。

6. Annotation Finalizer

ここで遺伝子名を付ける。

GeMoMa AnnotationFinalizer g=GENOME.fasta \
 a=OUTPUT/filtered_predictions_1.gff \
 u=NO \
 tf=true \
 rename=COMPOSED \
 p=遺伝子の接頭文字 \
 infix=G \
 s=00 \
 crp="chr,contig,scf" \
 d=遺伝子ID番号の桁数  \
 outdir=OUTPUT

p= で遺伝子名の Prefix を指定する。イネの **O++s とかダイズの Glyma みたいなやつ。
infix= は Glyma01G の G　にあたるやつ。
s= は遺伝IDの最後に付ける文字。00 が普通じゃない？
crp=　は、コンティグ名・スキャフォールド名から削りたい文字を指定。何も指定しない場合には Glyma.chromosome01G0000010 みたいになる。多分多くの人は遺伝子IDに chromosome はない方を好むだろう。コンマ区切りで幾つでも指定可能。あとから sed やテキストエディタの置換機能を使って一括処理してもいいけどね。
d= は遺伝子ID番号の桁数。 p=Glyma. infix=G s=00 d=5 だと Glyma.Chr01G0000100 みたいな感じになる。

なお、無論 'u=NO tr=true' の代わりに u=YES i=OUPUT/denoised_introns.fasta c=UNSTRANDED coverage_unstranded=OUTPUT/coverage.bedgraph を指定すれば、あらためて UTR を付け直すことができる。
また、BRAKER 以外のアノテーション結果や、別の近縁種を GeMoMa で移植したアノテーションファイルからそれぞれ差分を取って追加することも可能だ。どこまでやるかは本人次第だが、やり過ぎても却って良くない結果になるのでほどほどにしておこう。