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 (ノイズ除去) をマッピングした CLEAN.bam

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

1. インストールと準備諸々

要件として Java-1.8 以降と tblastn または 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 でイケるはず。

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 などでインストール。
以下、一応すべてのコマンドをパイプライン化したスクリプトもあるのだが、作者の想定した通りにインプットファイルを用意したりオプションを設定したりするのはなかなか難しく、途中でエラーを吐かれてしまいがちなので一つずつ進める。

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 ができる。

3. Check & Denoise Introns

exon-intron junction の配列をチェック

GeMoMa CheckIntrons t=GENOME.fasta i=OUTPUT/introns.gff outdir=OUTPUT

怪しいイントロンの除去

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 ができる。以降、イントロンファイルを求められるコマンドではこれを使う。

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

相同性検索

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 が推奨するパラメータ設定だけど、かなり緩い検索条件だなという印象はある。だが結果としてうまくいくので信じておく。この時点で出力される結果は人には読めないので次のコマンドで変換する。

結果の変換

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 が得られる。

7. GeMoMa Annotation Filter

GeMoMa GAF g=OUTPUT/predicted_annotation.gff outdir=OUTPUT

信頼性の低いアノテーションはここで除去。結果として filtered_predictions.gff が得られる。

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 が変な場所に予測してしまった遺伝子を弾くことができるのだが、それは次のステップで行う。この回はここまで。

関連

ゲノムアノテーション（目次）
RepeatModeler/RepeatMaskerを使ってゲノムにSoftmaskを掛ける
BRAKER2のインストール
BRAKER2を使ってタンパク質コード領域を予測する
RNA-seq のノイズを除去してゲノムにマッピングする
Gffcompare を使ってGeMoMaとBrakerの差分を取り、アノテーションを補完する
EnTAP のセットアップ
EnTAP の実行と結果の解釈