真核生物の遺伝子モデルの作り方素案

#このページの目的
ゲノムからの遺伝子モデルはいつも悩ましい問題です。特に真核生物ではイントロンがあることなどからモデリングは複雑になり、完全無欠のゴールドスタンダードは存在しないと思います。そこで、ここでは、私が使用している遺伝子モデル構築法を晒すことで、初心者の方への参考や、達人の方からのブラッシュアップを受けたいと思っています。

#概要
複数の遺伝子モデリングソフトの結果をEVidenceModelerを使って融合させる形で予測を行います。使用ソフトは以下です

+ [augustus-3.2.1](http://bioinf.uni-greifswald.de/augustus/)

+ [GeneMark-ES v.4.33](http://exon.gatech.edu/GeneMark/)
+ [glimmerHMM 3.0.4](https://ccb.jhu.edu/software/glimmerhmm/)
+ [exonerate version 2.4.7](http://www.ebi.ac.uk/about/vertebrate-genomics/software/exonerate-manual)
+ [EVidenceModeler v1.1.1](https://evidencemodeler.github.io)

また以下のシステムを補助的に使います

+ [WebApollo 2.0.6](http://genomearchitect.github.io)
+ [fasta-splitter](http://kirill-kryukov.com/study/tools/fasta-splitter/)
+ [The Genome Annotation Library (gtf2gff3)](http://www.sequenceontology.org/software/GAL.html)

#**I.ExonerateでEST,RNAseqデータをゲノムにマッピング**
トランスクリプトームデータをゲノムにマップすることで遺伝子領域を予測します。
RNAseqやESTデータは、対応領域が転写されていることの証拠となり、該当領域が遺伝子であることを示唆します。ゲノムへのcDNA配列のアライメントにはintron領域を考慮する必要がありますが、Exonerateはこれを行うのに適したソフトです。これ以外にPASAやGSNAPを用いることも出来ると思いますが、私はexonerateを使い慣れているのでこれを使いました。

###インストール法
省略します

###1.trinityからcomplete CDSだけを取り出す

Trintiy アセンブリングを行い、transdecoderでタンパク質コード領域を取り出します。その後、complte (開始点、終止点が含まれたCDS配列）を持っている物だけを取り出して、クエリセットとします。
[complete配列を取り出すスクリプト](http://qiita.com/MaedaTaro_Umiushi/items/be1811a040ab17a6bd39)

###2.高い類似性を示すベストヒットをgff形式で取り出す

以下の設定でexonerateをかける。


```sh:exonerate.sh
#!/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のアレイジョブなどをつかっってかける

```sh:split.sh
#!/bin/sh

fasta-splitter.pl \
--part-size 1000 \
--nopad \
./fasta/Trinity.0.95.fasta.transdecoder.comp.cds

```

```sh:exonerate.array.sh
#!/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

```sh:exone_EVM_comv.sh
#!/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によって自分で予測モデルの作成が出来るため、他の予測ソフトの項で出てくる手間や人為的バイアスの可能性を除去できます。(各ソフトの精度については[ここ](http://augustus.gobics.de/accuracy)などが便利かと)

コマンドはシンプルです。

```sh:GeneMark_run.sh
#!/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用に変換します


```py:EVM_gff_converter.genmark.py

#!/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はここでの遺伝子モデル作成に有効なツールです。

[webapolloの導入方法](http://qiita.com/MaedaTaro_Umiushi/items/3adf70e7264cd13d229d)

####Glimmerのトレーニング
Webapolloなどで作成した高信頼性遺伝子のgffファイルをGlimmer用の形式に変換します。

以下の様に、コンティグ名とexon（CDS)領域がタブ区切りで並んだ形式で、遺伝子の変わり目には空欄行をいれます。
以下の例ならば、上の２行は2exonをもつ一つの遺伝子をしめしています。また逆鎖にコードされている場合は、3番目の遺伝子の様にCDSのendとstartを逆鎖に対応させることで表現します。

```txt.cread_gene.txt

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だけ）を使って以下のコマンドを実行し、予測モデルが入ったディレクトリを作成させます。

```sh:training.sh

#!/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を実行します。


```sh:glimmer_Run.sh

#!/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にして回避します。

```sh:single_fasta.sh

#!/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による遺伝子予測

長いので[別立て](http://qiita.com/MaedaTaro_Umiushi/items/407dc5b283aada13b64d)にします



#**III.EvidenceModelerによるモデルの統合**

複数の遺伝子モデルから、信頼性を評価し、手法毎に重み付けをして、各lociごとに統合された遺伝子モデルを構築します。

##gffファイルの準備

上記の手法で作成したgffファイルを準備し、遺伝子予測ソフトによるgffファイルは1つにまとめる。exonerateによるgffは別ファイルとしておく。

```sh:gff_merge.sh

#!/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の精度などを加味して決める。

```txt:weight.txt
ABINITIO_PREDICTION	Augustus	2
ABINITIO_PREDICTION	GlimmerHMM	1
ABINITIO_PREDICTION	GeneMark.hmm	1
TRANSCRIPT	match	10

```


##ゲノム領域を切り取って短い単位でEVM解析をする


ゲノムをのりしろ用の領域をつけながら100Kbp位に分断してそれぞれ別ディレクトリにする

```sh:EVM_split.sh
#!/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  ##分割したファイルの情報が入った出力ファイル

```

各ディレクトリ用の実行ファイルを作成する


```sh:EVM_split_command.sh
#!/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個ぐらいに分割する。

```sh:split.sh
#!/bin/sh

split -d -l n commands.list

```

```sh:EVM_run.sh

#!/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に変換してさらにマージするスクリプト

```sh:EVM_merge.sh

#!/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

```