GffCompare を使って GeMoMa と BRAKER の差分を取り、アノテーションを補完する

さて、[BRAKER2 による遺伝子アノテーション](https://qiita.com/drk0311/items/0263691a7cc4dcb0c33c/) と、[GeMoMa による近縁種のアノテーション移植](https://qiita.com/drk0311/items/32aa45b85e0ad3545067) が終わったら、それらを合わせしまう。ただ、単純に両方を合わせるだけだと同一ローカスに異なる 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　先生の解説ページ](https://kazumaxneo.hatenablog.com/entry/2019/08/26/073000) で。
<br>

# 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` だ。タブ区切りテクストなのでエクセルから開くこともできる。

![IMG-3026.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2419328/76213a93-f0d8-a145-271d-119231791399.jpeg)

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

# 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
```
<br>

# 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 に掛ける。
<br>

# 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 で移植したアノテーションファイルからそれぞれ差分を取って追加することも可能だ。どこまでやるかは本人次第だが、やり過ぎても却って良くない結果になるのでほどほどにしておこう。

以上で、ゲノムアノテーションは終了である。これで格好は付いたでしょ？
次回は、予測した遺伝子の機能アノテーションとなる。
<br>

# 関連
- [ゲノムアノテーション（目次）](https://qiita.com/drk0311/items/c024888962f840e3ce3e)
- [RepeatModeler/RepeatMaskerを使ってゲノムにSoftmaskを掛ける](https://qiita.com/drk0311/items/cba7958319737a1f12c0)
- [BRAKER2のインストール](https://qiita.com/drk0311/items/a3ac648f2780cfee57b1)
- [BRAKER2を使ってタンパク質コード領域を予測する](https://qiita.com/drk0311/items/0263691a7cc4dcb0c33c/)
- [RNA-seq のノイズを除去してゲノムにマッピングする](https://qiita.com/drk0311/items/56715d46dc7893c95f78)
- [GeMoMa で近縁種のアノテーションを移植する](https://qiita.com/drk0311/items/32aa45b85e0ad3545067) 
- [EnTAP のセットアップ](https://qiita.com/drk0311/items/9bbe56dc0669376eb98f)
- [EnTAP の実行と結果の解釈](https://t.co/kWNa01djMF)