GeMoMa で近縁種のアノテーションを移植する

GeMoMa (http://www.jstacs.de/index.php/GeMoMa) は同種あるいは近縁種など、すでに存在する Reference ゲノムから遺伝子アノテーションを移すためのツールである。コードされるタンパク質の配列をマッピングすることによって遺伝子を予測するが、そのときに RNA-seq のデータを使ってイントロン領域を予測したりするので精度・感度ともに高い。単なる転写産物のマッピングに基づくGeneLift 等はシンプルで使いやすいが、私としては圧倒的に GeMoMa を勧める。GeMoMa の操作で打ち込むコマンドは結構沢山あり、また一つ一つが長くて打ち込むのが面倒ではあるが、計算自体はどれも短時間で終わる。一日頑張れば全てのプロセスを完遂できるだろう。

# 用意するもの
以下の全てが単一のディレクトリに保存されているものとする。
- Reference として使う近縁種のゲノム `REFERENCE.fasta` 
- Reference のアノテーションファイル `REFERENCE.gtf` または `REFERENCE.gff` 
- 自分でアセンブルしたゲノム `GENOME.fasta` 
- [RNA-seq (ノイズ除去) をマッピングした](https://qiita.com/drk0311/items/56715d46dc7893c95f78) `CLEAN.bam` 

なお、この解説ではWSL環境の `/mnt/c/アカウント名/GEMOMA/` が作業ディレクトリであると想定しておく。

# 1. インストールと準備諸々
要件として [Java-1.8 以降](https://openjdk.java.net/install/) と [tblastn](https://blast.ncbi.nlm.nih.gov/Blast.cgi?PAGE_TYPE=BlastDocs&DOC_TYPE=Download) または [mmseqs2](https://github.com/soedinglab/MMseqs2) が必要。
```
mamba create -n gemoma java-1.8.0-openjdk-devel-cos7-s390x mmseqs2
conda activate gemoma
```
私は anaconda 版の Java は試していないので、ダメそうだったら https://openjdk.java.net/install/ を参照して自分でインストールしよう。 Ubuntu なら `sudo apt-get install openjdk-8-jre` でイケるはず。
<br>

GeMoMa は Anaconda パッケージも存在するが、V1.7（現時点での最新は 1.8 ）なので http://www.jstacs.de/index.php/GeMoMa から最新版を入手。ダウンロードして解凍すればすぐ使える。この解説中では conda 以外の方法でインストールするツールは `~/local` 中に置くこととする。 GeMoMa をダウンロードして解凍した後、元ファイルは削除して作業ディレクトリに移動、また結果の出力用に `OUTPUT` ディレクトリを作る。
```
cd ~/local
mkdir gemoma
cd gemoma

wget http://www.jstacs.de/download.php?which=GeMoMa
unzip GeMoMa-1.8.zip
rm GeMoMa-1.8.zip

cd /mnt/c/アカウント名/GEMOMA
mkdir OUTPUT
```
unzip コマンドが入っていなかった場合は `sudo apt install unzip` などでインストール。
以下、一応すべてのコマンドをパイプライン化したスクリプトもあるのだが、作者の想定した通りにインプットファイルを用意したりオプションを設定したりするのはなかなか難しく、途中でエラーを吐かれてしまいがちなので一つずつ進める。
<br>


# 2. Extract RNA-seq Evidence
要するにRNA-seqをマッピングしたBAMファイルからリードの厚みとイントロンを抽出する。
```
GeMoMa ERE s=FR_UNSTRANDED m=CLEAN.bam outdir=OUTPUT
```
stranded RNA-seq を実施していた場合は、

- 1st strand library なら `s=FR_FIRST_STRAND`
- 2nd strand library なら `s=FR_SECOND_STRAND` 

とすること。 分からない場合は業者さんに聞こう（いまどきは 2nd strand library が多いと思うけど）。
結果として `introns.gff` と `coverage.bedgraph` ができる。 
<br>


# 3. Check & Denoise Introns
**exon-intron junction の配列をチェック**
```
GeMoMa CheckIntrons t=GENOME.fasta i=OUTPUT/introns.gff outdir=OUTPUT
```
<br>

**怪しいイントロンの除去**
```
GeMoMa DenoiseIntrons i=OUTPUT/introns.gff c=UNSTRANDED coverage_unstranded=OUTPUT/coverage.bedgraph outdir=OUTPUT
```
stranded RNA-seq を実施していた場合は、`c=STRANDED`　として、`coverage_unstranded` のところを `coverage_forward=OUTPUT/coverage_forward.bedgraph coverage_reverse=OUTPUT/coverage_reverse.bedgraph` とすること（注意：bedgrph のファイル名は違っているかも）。
結果、`denoised_introns.gff` ができる。以降、イントロンファイルを求められるコマンドではこれを使う。
<br>

# 4. Extract CDS Sequences from Reference
遺伝子座からタンパク質をコードした領域だけを抜き出す。
```
GeMoMa Extractor a=REFERENCE.gtf g=REFERENCE.fasta outdir=OUTPUT
```
結果、`cds-parts.fasta` と `assignment.tabular` ができる。


# 5. REFERENCE の CDS を GENOME.fasta にマップする
速いので mmseqs2 を使う。便利なツールなので意識の高い者は自習してモノにせよ。
**mmseqs2 データベースの作成**
```
mkdir TARGET
mmseqs createdb GENOME.fasta TARGET/GenomeDB -v 2
mkdir QUERY
mmseqs createdb OUTPUT/cds-parts.fasta QUERY/cdsDB -v 2
```
<br>

**相同性検索**
```
mmseqs search QUERY/cdsDB TARGET/GenomeDB OUTPUT/mmseqs.out OUTPUT/mmseqs_tmp -e 100.0 --threads 24 -s 8.5 -a --comp-bias-corr 0 --max-seqs 500 --mask 0 --orf-start-mode 1 -v 2
```
GeMoMa が推奨するパラメータ設定だけど、かなり緩い検索条件だなという印象はある。だが結果としてうまくいくので信じておく。この時点で出力される結果は人には読めないので次のコマンドで変換する。
<br>

**結果の変換**
```
mmseqs convertalis QUERY/cdsDB TARGET/GenomeDB OUTPUT/mmseqs.out OUTPUT/search.txt --threads 24 --format-output "query,target,pident,alnlen,mismatch,gapopen,qstart,qend,tstart,tend,evalue,bits,empty,raw,nident,empty,empty,empty,qframe,tframe,qaln,taln,qlen,tlen" -v 2
```
`--format-output` ですごく長いパラメータを指定しているが、結果を開けてみれば何が起きているかは分かると思う。blast 検索の `-outmft 6` でタブ区切り出力を経験したことがあればなおよく分かるはず。 「？？？」となった人は https://www.metagenomics.wiki/tools/blast/blastn-output-format-6 を参照。
これによって得られた `search.txt` が次のインプットとなる。

# 6. Gene Model Mapper
ここでついに GeMoMa の核心である。相同性検索の結果をもとに、GENOME.fasta に遺伝子アノテーションを入れる (この時点では CDS の予測のみ)
```
GeMoMa GeMoMa s=OUTPUT/search.txt \
  c=OUTPUT/cds-parts.fasta a=OUTPUT/assignment.tabular \
  t=GENOME.fasta sort=true Score=ReAlign i=OUTPUT/denoised_introns.gff \
  coverage=UNSTRANDED coverage_unstranded=OUTPUT/coverage.bedgraph \
  outdir=OUTPUT
```
stranded RNA-seq を使っていた場合は、
`coverage=STRANDED converage_forward=OUTPUT/coverage_forward.bedgraph coverage_reverse=OUTPUT/coverage_reverse.bedgraph` とすること。

結果、`predicted_annotation.gff` が得られる。
<br>

# 7. GeMoMa Annotation Filter
```
GeMoMa GAF g=OUTPUT/predicted_annotation.gff outdir=OUTPUT
```
信頼性の低いアノテーションはここで除去。結果として `filtered_predictions.gff` が得られる。
<br>

# 8. Annotation Finalizer
RNA-seq の結果を使って UTR を付加する。遺伝子座のネーミングもここで出来る（ここではまだやらない）。
```
GeMoMa AnnotationFinalizer g=GENOME.fasta a=OUTPUT/filtered_predictions.gff 
u=YES i=OUTPUT/denoised_introns.gff c=UNSTRANDED coverage_unstranded=OUTPUT/coverage.bedgraph rename=NO outdir=OUTPUT
```
stranded RNA-seq の場合は（以下略
これで、UTR 付きで近縁種のアノテーションを移植した `final_annotation.gff` が出来る。UTR があることで、BRAKER2 が変な場所に予測してしまった遺伝子を弾くことができるのだが、それは次のステップで行う。この回はここまで。
<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)
- [Gffcompare を使ってGeMoMaとBrakerの差分を取り、アノテーションを補完する](https://qiita.com/drk0311/items/7c9910f12a2a959c710d)
- [EnTAP のセットアップ](https://qiita.com/drk0311/items/9bbe56dc0669376eb98f)
- [EnTAP の実行と結果の解釈](https://t.co/kWNa01djMF)