#このページの目的
ゲノムからの遺伝子モデルはいつも悩ましい問題です。特に真核生物ではイントロンがあることなどからモデリングは複雑になり、完全無欠のゴールドスタンダードは存在しないと思います。そこで、ここでは、私が使用している遺伝子モデル構築法を晒すことで、初心者の方への参考や、達人の方からのブラッシュアップを受けたいと思っています。
#概要
複数の遺伝子モデリングソフトの結果をEVidenceModelerを使って融合させる形で予測を行います。使用ソフトは以下です
また以下のシステムを補助的に使います
#I.ExonerateでEST,RNAseqデータをゲノムにマッピング
トランスクリプトームデータをゲノムにマップすることで遺伝子領域を予測します。
RNAseqやESTデータは、対応領域が転写されていることの証拠となり、該当領域が遺伝子であることを示唆します。ゲノムへのcDNA配列のアライメントにはintron領域を考慮する必要がありますが、Exonerateはこれを行うのに適したソフトです。これ以外にPASAやGSNAPを用いることも出来ると思いますが、私はexonerateを使い慣れているのでこれを使いました。
###インストール法
省略します
###1.trinityからcomplete CDSだけを取り出す
Trintiy アセンブリングを行い、transdecoderでタンパク質コード領域を取り出します。その後、complte (開始点、終止点が含まれたCDS配列）を持っている物だけを取り出して、クエリセットとします。
complete配列を取り出すスクリプト
###2.高い類似性を示すベストヒットをgff形式で取り出す
以下の設定でexonerateをかける。
#!/bin/sh
#$ -S /bin/bash
#$ -cwd
#$ -v PATH
exonerate \
--model est2genome \
--bestn 1 \
--percent 90 \
--showtargetgff \
--query \
./fasta/Trinity.0.95.fasta.transdecoder.comp.cds \
--target \
/home/maedat/Eor_2016/over2000_pla_150327_out_gapClosed_subset.fa \
クエリ配列が多い場合は、fasta-splitterでクエリを分割した後
SGEのアレイジョブなどをつかっってかける
#!/bin/sh
fasta-splitter.pl \
--part-size 1000 \
--nopad \
./fasta/Trinity.0.95.fasta.transdecoder.comp.cds
#!/bin/sh
#$ -S /bin/bash
#$ -cwd
#$ -v PATH
#$ -t 1-1000
exonerate \
--model est2genome \
--bestn 1 \
--percent 90 \
--showtargetgff \
--query \
./fasta/Trinity.0.95.fasta.transdecoder.comp.part-$SGE_TASK_ID.cds \
--target \
/home/maedat/Eor_2016/over2000_pla_150327_out_gapClosed_subset.fa \
アレイを使った場合は、出力をcatで結合して
all.exonerate.out.txt
などとしておきます。
###3.マッチしたデータをEVM用に変換する
EVidenceModelerはGFF3形式を要求しますが、exonerateが吐くgff3はEVMの形式とマッチしないため、変換用のスクリプトを使って変換します
変換用スクリプトはEVidenceModelerにくっついてきます
/EVidenceModeler-1.1.1/EvmUtils/misc/Exonerate_to_evm_gff3.pl
#!/bin/sh
perl \
~/opt/EVidenceModeler-1.1.1/EvmUtils/misc/Exonerate_to_evm_gff3.pl \
all.exonerate.out.txt > all.exonerate.out.EVM.gff
#II.遺伝子モデル予測
Trinityのマッピングだけでは低発現遺伝子などを取り逃すので遺伝子モデリングソフトを使って領域を予想します。
###1, GeneMark-ESによる遺伝子予測
GeneMark-ESは、self-trainingによって自分で予測モデルの作成が出来るため、他の予測ソフトの項で出てくる手間や人為的バイアスの可能性を除去できます。(各ソフトの精度についてはここなどが便利かと)
コマンドはシンプルです。
#!/bin/sh
#$ -S /bin/bash
#$ -cwd
#$ -v PATH
perl ~/opt/gm_et_linux_64/gmes_petap/gmes_petap.pl \
--ES --cores 20 --sequence \
./genome_sequence.fa
出力はgtfなのでこれを
The Genome Annotation Library (gtf2gff3)
についてくるgtf2gff3スクリプトを使ってgff3に変換し、さらに以下の私自作のスクリプトでEVM用に変換します
#!/usr/bin/env python
# coding:utf-8
###各種遺伝子予想ソフトの吐くgff3をEVMの対応形式に変換するgenmark版
# scaffold22744_cov31 GeneMark.hmm gene 1341 2362 . + . ID=38108_g;Name=38108_g
# scaffold22744_cov31 GeneMark.hmm mRNA 1341 2362 . + . ID=38108_t;Parent=38108_g
# scaffold22744_cov31 GeneMark.hmm start_codon 1341 1343 . + 0 ID=start_codon:38108_t:1;Parent=38108_t
# scaffold22744_cov31 GeneMark.hmm exon 1341 1498 0 + . ID=exon:38108_t:1;Parent=38108_t
# scaffold22744_cov31 GeneMark.hmm exon 1722 1800 0 + . ID=exon:38108_t:2;Parent=38108_t
# scaffold22744_cov31 GeneMark.hmm exon 2114 2362 0 + . ID=exon:38108_t:3;Parent=38108_t
# scaffold22744_cov31 GeneMark.hmm CDS 1341 1498 . + 0 ID=CDS:38108_t:1;Parent=38108_t
# scaffold22744_cov31 GeneMark.hmm CDS 1722 1800 . + 1 ID=CDS:38108_t:2;Parent=38108_t
# scaffold22744_cov31 GeneMark.hmm CDS 2114 2362 . + 0 ID=CDS:38108_t:3;Parent=38108_t
# scaffold22744_cov31 GeneMark.hmm stop_codon 2360 2362 . + 0 ID=stop_codon:38108_t:1;Parent=38108_t
#あらかじめ、必要な行だけ取っておく
#grep -E "gene|CDS|mRNA" genemark.gff > genemark.edit.gff
import sys
import csv
import re
ori_gff_in = sys.argv[1] #１番目の引数には、変更したいgffファイルを指定する。
#ID項目を取り出すための正規表現をコンパイル
ID_get = re.compile("ID=(.+?);")
#Parent項目を取り出すための正規表現をコンパイル
Parent_get = re.compile("Parent=(.+)")
for q in open(ori_gff_in, "r"): #gffファイルを開く
ori_gff = q.split('\t') #タブで区切ってori_gffリストにいれる
#print(ori_gff)
if ori_gff[2] == "gene": #3番目のtype項目がgeneだったら
gene_line = '\t'.join([str(i) for i in ori_gff])
print(gene_line) #標準出力にそのまま出力
CDS_count=0 #CDSカウントを0に戻す
elif ori_gff[2] == "mRNA": #3番目のtype項目がtranscriptだったら
mRNA_gff = ori_gff[:] #mRNA項目用のリストを作成（コピー）
mRNA_line = '\t'.join([str(i) for i in mRNA_gff])
print(mRNA_line) #変更した内容を出力
elif ori_gff[2] == "CDS": #3番目のtype項目がCDSだったら
CDS_count += 1 #CDSカウントをインクリメント
exon_gff = ori_gff[:] #exon項目用のリストを作成（コピー）
exon_gff[2] = "exon" #exonにtype項目を変更
description = ori_gff[8] #9番目のdescriptionを取り出す
ID_mutch = ID_get.search(description) #IDを取り出す
NEW_ID = ID_mutch.group(1) + "_" + str(CDS_count) #IDにCDSの番号を振って一意にする
Parent_mutch = Parent_get.search(description) #Parentを取り出す
NEW_exon_description = "ID=" + NEW_ID + "_CDS;Parent=" + Parent_mutch.group(1) + "\n"
exon_gff[8] = NEW_exon_description
exon_line = '\t'.join([str(i) for i in exon_gff])
print(exon_line) #変更した内容を出力
CDS_gff = ori_gff[:] #CDS項目用のリストを作成（コピー）
description = exon_gff[8] #9番目のdescriptionを取り出す
ID_mutch = ID_get.search(description) #IDを取り出す（exonのIDになっているのでこれをparentにする）
CDS_description = "ID=" + ID_mutch.group(1) + "_CDS;Parent=" + ID_mutch.group(1) + "\n"
CDS_gff[8] = CDS_description
CDS_line = '\t'.join([str(i) for i in CDS_gff])
print(CDS_line) #変更した内容を出力
###2, glimmerHMMによる遺伝子予測
####マニュアルによる高信頼性遺伝子データの作成
近縁種のゲノム解読がなされていない場合、遺伝子予想モデル（パラメーター）から新規に作成する必要があります。augustus とglimmerHMMでは信頼出来る遺伝子情報を教師として、モデルをトレーニングする必要があります。そこでまず、信頼出来る遺伝子モデルを目視で作成します。最長contigなどを使って、100個以上の遺伝子を作ります。Trinityベースで作ったマッチングデータや、公共データベース（Refseq, uniprot)をクエリとしたtblaspやexonerate解析、RNAseqリードのtophatによるマッピングなどを駆使して過不足のない遺伝子モデルを作りましょう。WebApolloはここでの遺伝子モデル作成に有効なツールです。
####Glimmerのトレーニング
Webapolloなどで作成した高信頼性遺伝子のgffファイルをGlimmer用の形式に変換します。
以下の様に、コンティグ名とexon（CDS)領域がタブ区切りで並んだ形式で、遺伝子の変わり目には空欄行をいれます。
以下の例ならば、上の２行は2exonをもつ一つの遺伝子をしめしています。また逆鎖にコードされている場合は、3番目の遺伝子の様にCDSのendとstartを逆鎖に対応させることで表現します。
scaffold1408_cov40 1422843 1422916
scaffold1408_cov40 1423183 1423387
scaffold1408_cov40 1402083 1402148
scaffold1408_cov40 1402914 1403018
scaffold1408_cov40 1404643 1404764
scaffold1408_cov40 1406443 1406515
scaffold1408_cov40 1409733 1409810
scaffold1408_cov40 1410677 1410771
scaffold1408_cov40 1413277 1413389
scaffold1408_cov40 1413926 1414061
scaffold1408_cov40 1414601 1414706
scaffold1408_cov40 1416114 1416209
scaffold1408_cov40 894468 894363
scaffold1408_cov40 888957 888866
scaffold1408_cov40 883683 883515
scaffold1408_cov40 879683 879583
scaffold1408_cov40 877427 877353
scaffold1408_cov40 873904 873518
以上のアノテーションファイルと、対応するゲノム配列（高信頼性遺伝子の選抜に使った最長contigだけ）を使って以下のコマンドを実行し、予測モデルが入ったディレクトリを作成させます。
#!/bin/sh
#$ -S /bin/bash
#$ -cwd
#$ -v PATH
~/opt/GlimmerHMM/train/trainGlimmerHMM \
./seq.fa \
./seq_anno.txt
####GlimmerのRun
以上で作成した予測モデルディレクトリ(ここではTrainGlimmM2017-04-14D17:27:52)を使ってRUNを実行します。
#!/bin/sh
#$ -S /bin/bash
#$ -cwd
#$ -v PATH
~/opt/GlimmerHMM/bin/glimmerhmm \
./all.genome.fa \
./TrainGlimmM2017-04-14D17:27:52 \
-o gene.glimmer.out \
-g #gff形式で出力するときに指定
Glimmerはゲノム側がmutiple fastaであるときに、一番最初の配列しか認識しないため
複数のcontigからなるゲノムの場合は以下の様にしてsingle fastaにして回避します。
#!/bin/sh
CONT_NUM=`grep ">" $1 |wc -l`
echo $CONT_NUM
fasta-splitter.pl \
--n-parts $CONT_NUM \
--measure count \
--nopad \
$1
連番になっているsingle fastaファイルをSGEのarrayを使ってRUNすれば素早く処理が終わります。
得られた出力ファイルはcatで結合させて、all.glimmer.out.txt などと名前をつけましょう。
####EVM用のフォーマットへ変換
以下のツールでEVM用に形式を変更します。
/EVidenceModeler-1.1.1/EvmUtils/misc/glimmerHMM_to_GFF3.pl
使い方はExonerate_to_evm_gff3.plと同じ
###3,Augustusによる遺伝子予測
長いので別立てにします
#III.EvidenceModelerによるモデルの統合
複数の遺伝子モデルから、信頼性を評価し、手法毎に重み付けをして、各lociごとに統合された遺伝子モデルを構築します。
##gffファイルの準備
上記の手法で作成したgffファイルを準備し、遺伝子予測ソフトによるgffファイルは1つにまとめる。exonerateによるgffは別ファイルとしておく。
#!/bin/sh
cat augustus.perl.EVM.gff3 exonerate.perl.EVM.gff genemark.edit.EVM.gff > predicted.gff
重み付け情報をweightファイルに記述する
タブ区切りで、遺伝子予想ファイルの情報は1列目に ABINITIO_PREDICTION として2列目に使用ソフトウエアの名前を入れる。2列目の名前はgffファイルの2列目の記述と揃える。3列目には重みを数字でいれる。大きいほど重要視する。cDNAマッピングの結果に関する重み付けもこのファイルに書く。その場合、1列目はTRANSCRIPTとする.重み付けは、元につかった予想ソフトやcDNAの精度などを加味して決める。
ABINITIO_PREDICTION Augustus 2
ABINITIO_PREDICTION GlimmerHMM 1
ABINITIO_PREDICTION GeneMark.hmm 1
TRANSCRIPT match 10
##ゲノム領域を切り取って短い単位でEVM解析をする
ゲノムをのりしろ用の領域をつけながら100Kbp位に分断してそれぞれ別ディレクトリにする
#!/bin/sh
#$ -S /bin/bash
#$ -cwd
#$ -v PATH
perl ~/opt/EVidenceModeler-1.1.1/EvmUtils/partition_EVM_inputs.pl \
--genome /home/user/data/genome.fa \
--gene_predictions \
/home/user/data/predicted.gff \ #フルパスで遺伝子モデルソフトで予想したgffファイルを指定する
--transcript_alignments \
/home/user/data/exonerate.perl.EVM.gff \ #フルパスでcDNAマッピング結果のgffファイルを指定する
--segmentSize 100000 \
--overlapSize 10000 \
--partition_listing partitions_list.out ##分割したファイルの情報が入った出力ファイル
各ディレクトリ用の実行ファイルを作成する
#!/bin/sh
#$ -S /bin/bash
#$ -cwd
#$ -v PATH
perl ~/opt/EVidenceModeler-1.1.1/EvmUtils/write_EVM_commands.pl \
--genome \
/home/user/data/genome.fa \
--weights \
/home/user/data/weights.txt \
--gene_predictions \
predicted.gff \ #ファイル名だけをいれる
--transcript_alignments \
exonerate.perl.EVM.gff \ #ファイル名だけをいれる
--output_file_name evm.out \
--partitions partitions_list.out \ #上の処理で作った分割したファイルの情報が入ったファイルを指定する。
> commands.list
commands.listには各ディレクトリ毎に実行するためのコマンドが1行づつたくさん入っている。
SGEでアレイ処理するために100個ぐらいに分割する。
#!/bin/sh
split -d -l n commands.list
#!/bin/sh
#$ -S /bin/bash
#$ -cwd
#$ -v PATH
#$ -t 1-99
perl \
~/opt/EVidenceModeler-1.1.1/EvmUtils/execute_EVM_commands.pl \
x$SGE_TASK_ID \
| tee run.part.$SGE_TASK_ID.log
完成したファイルをのりしろを考慮してマージした後、gffに変換してさらにマージするスクリプト
#!/bin/sh
#$ -S /bin/bash
#$ -cwd
#$ -v PATH
perl ~/opt/EVidenceModeler-1.1.1/EvmUtils/recombine_EVM_partial_outputs.pl \
--partitions partitions_list.out --output_file_name evm.out
perl ~opt/EVidenceModeler-1.1.1/EvmUtils/convert_EVM_outputs_to_GFF3.pl \
--partitions partitions_list.out \
--output evm.out \
--genome \
/home/user/data/genome.fa \
find . -regex ".*evm.out.gff3" -exec cat {} \; > EVM.all.gff3