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