LoginSignup
6
7

More than 5 years have passed since last update.

@MaedaTaro_Umiushi

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

Last updated at Posted at 2017-04-20

このページの目的

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

概要

複数の遺伝子モデリングソフトの結果を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をかける。

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のアレイジョブなどをつかっってかける

split.sh
#!/bin/sh

fasta-splitter.pl \
--part-size 1000 \
--nopad \
./fasta/Trinity.0.95.fasta.transdecoder.comp.cds
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

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

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

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用に変換します

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の導入方法

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

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

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

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

長いので別立てにします

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

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

gffファイルの準備

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

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

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

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

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

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  ##分割したファイルの情報が入った出力ファイル

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

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個ぐらいに分割する。

split.sh
#!/bin/sh

split -d -l n commands.list
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に変換してさらにマージするスクリプト

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
6
7
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
6
7