macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

複数の遺伝子モデルを統合する EVidenceModeler

2023/02/15 追記

 

 EVidenceModeler (EVM) は、真核生物の遺伝子構造を自動アノテーションするツールであり、真核生物の遺伝子構造を、利用可能なすべての証拠の重み付きコンセンサスとして報告するものである。EVMは、Program to Assemble Spliced Alignments (PASA)と組み合わせることにより、タンパク質コード遺伝子と交互にスプライスされたアイソフォームを予測する、包括的で設定可能なアノテーションシステムを実現することができる。イネとヒトのゲノム配列を用いた実験により、EVMは手動キュレーションに匹敵する品質の自動遺伝子構造アノテーションを生成することが実証された。

 

wiki

https://github.com/EVidenceModeler/EVidenceModeler/wiki

 

コメント

EVMはさまざまな遺伝子アノテーション結果(structural anotation)を統合するために使用されるツールです。最良のアノテーションは、十分に訓練された研究者が目でエビデンスを確認したり既知タンパク質と比較しながら進めるマニュアルアノテーションだと思いますが、それでは数万の遺伝子があるゲノムではスケールしません。自動で再現性のある計算機的手法が必要ということですね。発表は古いですが、調べた限り、最近のゲノムアセンブリ論文でも多くがEVMを使用しています。性能的に代替がほとんどないことを示唆しています。論文の第一著者はTrinityやPASA、Traansdecorderなどこの界隈で非常によく使われるソフトウェアを発表されているBrian J Haasさんです。忙しいと思いますがissuesでも積極的に応対されており、本当に凄い方だと思います。

インストール

Github

#conda(link)
mamba install -c bioconda evidencemodeler -y

> evidence_modeler.pl

$ evidence_modeler.pl 

 

################# Evidence Modeler ##############################

#

#  parameters:

#

#  Required:

# 

#  --genome              genome sequence in fasta format

#  --weights              weights for evidence types file

#  --gene_predictions     gene predictions gff3 file

#

#  Optional but recommended:

#  --protein_alignments   protein alignments gff3 file

#  --transcript_alignments       transcript alignments gff3 file

#

#  Optional and miscellaneous

#

#  --repeats               gff3 file with repeats masked from genome file

#

#  

#  --terminalExons         supplementary file of additional terminal exons to consider (from PASA long-orfs)

#

#  --stop_codons             list of stop codons (default: TAA,TGA,TAG)

#                               *for Tetrahymena, set --stop_codons TGA

#  --min_intron_length       minimum length for an intron (default 20 bp)

#  --exec_dir                directory that EVM cd's to before running.

#

# flags:

#

#  --forwardStrandOnly   runs only on the forward strand

#  --reverseStrandOnly   runs only on the reverse strand

#

#  -S                    verbose flag

#  --debug               debug mode, writes lots of extra files.

#  --report_ELM          report the eliminated EVM preds too.

#

#  --search_long_introns  <int>  when set, reexamines long introns (can find nested genes, but also can result in FPs) (default: 0 (off))

#

#

#  --re_search_intergenic <int>  when set, reexamines intergenic regions of minimum length (can add FPs) (default: 0  (off))

#  --terminal_intergenic_re_search <int>   reexamines genomic regions outside of the span of all predicted genes (default: 10000)

#

#################################################################################################################################

> ./EVM_to_GFF3.pl -h

usage: ./EVM_to_GFF3.pl evm_output_file contig_accession 

#こちらはEvmUtils/に含まれているユーティリティスクリプト。condaでは導入されない。

 

 

テストラン

git clone https://github.com/EVidenceModeler/EVidenceModeler.git
cd EVidenceModeler/simple_example/
bash runMe.sh

出力

以下の2つのコマンドが実行されている。

evidence_modeler.pl --genome genome.fasta \
--weights ./weights.txt \
--gene_predictions gene_predictions.gff3 \
--protein_alignments protein_alignments.gff3 \
--transcript_alignments transcript_alignments.gff3 > evm.out 

../EvmUtils/EVM_to_GFF3.pl evm.out.orig Contig1 > evm.out.gff3

weights.txt

1列目には予測の種類;外部エビデンスのPROTEINアラインメントかRNA seqアラインメント、もしくはab initioかどうか記載する。2列目にはソフトウエアの種類(指定するGFF3ファイル中の名前と一致している必要がある)、3列目には重みの数値を記載する(*1)。 --gene_predictions、 --protein_alignments、 --transcript_alignmentsそれぞれで指定するGFF3ファイルが複数ある場合、あらかじめ結合して指定する。EVM に提供する各エビデンスの重みは、2008年の論文のTable1にある4 つのクラスに分けられている(PMC)。重み値のデフォルト値(一般に推奨される値)は、wikiに書かれてる。それによると、ABINITIO_PREDICTIONの予測はデフォルト1、PROTEINのタンパク質アラインメントはデフォルト1だがGenewiseなどを使った高品質な遺伝子アラインメントがある時は2〜5、TRANSCRIPTの関連する生物種のESTのアラインメントは1、PASAアラインメントアセンブル10、信頼できる完全長遺伝子モデルが提供できる時はOTHER_PREDICTIONクラスで任意の数値を使うとなっている。

 

著者のBrian J Haasさん(TrinityやTransDecoderの開発者でもある)はEVM用のGFF3ファイルにするスクリプトをリクエストに応じて作成されており、Evmtilsにはそのためのスクリプトが配置されている(link)。

#EVMがパース可能なGFF3かどうか検証するスクリプト
perl EvmUtils/gff3_gene_prediction_file_validator.pl input.gff3 2> error

変換用スクリプト集はmisc(その他)に多く配置されている。

https://github.com/EVidenceModeler/EVidenceModeler/tree/master/EvmUtils/misc

 

より実践的には、EVMのwikiやこのレポジトリで説明されているように、複数に分けて並列ランし、結果を統合、gff3に再変換する。そのためのスクリプトもEvmUtils/に含まれている。

https://github.com/galikbence/genome_annotation

追記

2023年1月にアナウンスされたEVM2.0では1コマンドで分割・EVMの計算、再結合までできるようになっています。

 

issueより

  1. EVMがアライメントデータしか持っていないのであれば、遺伝子をモデル化する能力は非常に限られてしまう(PROTEIN および TRANSCRIPT クラスは開始コドンから停止コドンまでの完全な遺伝子構造を提供するのではなく、内部エクソンなどの遺伝子構造の構成要素を提供することが期待されると論文に書かれている。アライメントデータは遺伝子モデルを作っているのではなく、exonを予測してそれが遺伝子モデルぽく見えているに過ぎない。そのため、開始/停止コドンなども書かれていない。これらのエビデンスをベースとしては進めないという事)。十分に訓練されたab.iniito遺伝子予測結果を準備し(例;braker2は標準の設定でAUGUSTUSを5回訓練して使用する)、これをアラインメントベースのエビデンスと共にEVMに提供する(#39)。コメント;逆に言えば、ab.iniito遺伝子予測結果だけでも遺伝子は効率的にモデル化されない。
  2. 各サブフォルダに含まれるevm.out.logファイルを調べると、読み込めないエビデンスファイルがあるのかどうか確認できる。例えば"-WARNING -, skipping predType: GeneMark"と書いてあれば、そのgff3ファイルは読み込みエラーになっている(例えばweightファイルに書いてない、フォーマットがおかしいなど)。
  3. EVMはUTRをモデル化せず、コーディングエクソンだけを遺伝子モデルとして扱う(#46)。UTRを追加したい場合は、EVMのデータを入力アノテーションとして、後でPASAを実行する。
  4. EVMは遺伝子がコンティグの末端に位置しない限り完全なORFを必要とする
  5. #46 )。
  6. EVMは、入力されるすべてのエビデンスが最高クラスの品質であることを期待する。必ずしもそうでないエビデンスは排除することも考える(#15)。
  7. weightファイルで提供する、アライメントと予測のエビデンスタイプは扱い方が大きく異なる。ABINITIO_PREDICTIONには、完全なエクソン構造コンポーネントがある場合のみ使用することが推奨される(#22)。
  8. 1と重複するが、タンパク質アライメントデータは、内部エクソンを補完するためにのみ使用される。EVMはタンパク質アライメントの証拠だけでは遺伝子をモデル化しない。ただし、信頼度が高く完全なものを「OTHER_PREDICTION」のエビデンスタイプとして定義した場合は別である(#35)。
  9.  Transdecoderで得られたgffにはpartialの翻訳産物とcompleteの翻訳産物両方が含まれている。これをweightファイルのOTHER_PREDICTIONのエビデンスとして提供する時に、completeのみ提供すれば最高の品質のエビデンスを提供できる
  10. #22)。
  11. 最終出力の遺伝子モデル数が少なすぎる場合=> 入力フォーマットのエラー(#22)(No.2参照)。
  12. 最終的なアノテーションを得た後、低品質の遺伝子、例えば、少ない証拠でサポートされている遺伝子をどのようにフィルタリングするか=> 時間がかかるだろうが、BLAST結果と発現量の組み合わせである程度フィルタリングできる(#22)。
  13. エビデンスの間に多くの不一致がある場合 =>  最良の結果を得るために、EVMに最良のエビデンスのみを与えることを考える(#22)。(6と同じ)
  14. Stringtie2のGTF2.1は変換スクリプトが用意されていないが、cufflinks_gtf_to_alignment_gff3.plでEVM互換GFF3に変換できる(#38)。
  15. PacbioのIso-seqのカバレッジが非常に高い場合、それがスコアリングを支配している可能性がある。CD-HIt などでpacbioのリードの冗長性を無くし、代表的なpacbioのアラインメントををevmの入力として使用する。これがより良い結果に繋がるかどうかは1つの可能性でしかなく、十分に検討する必要がある(#39)。
  16. EVMで開始、ストップコドンがない遺伝子が報告されることがあるが、それはcontig末端にあるために部分的な遺伝子しかモデル化されないなどの現実的な問題に対処するためのEVMの仕様である。それでも削除するならEVM以外の方法か手動で行う必要がある(#47)。(4と重複)

 

コメント

  • エビデンスのGFF3ファイルを色々調整していると、weightファイルに記載漏れがでてくる可能性があります。統合した後のGFF3ファイルついて、weightファイルの2列目に書いた名前と一致しているか確認してからEVMをランしましょう。

#GFF3の2列目に書いてあるツールや由来生物種名を確認するなら

$ cut -f 2 gene.ab.initio.gff3 |sort - |uniq - 

  • どのような重み値を使えば良いのか悩ましいところですが、使用した限り、wikiで書かれている通りab.initioの重みは1, 2などと小さく設定し(ゲノムブラウザで見て高品質に見えたとしても)、近縁種のタンパク質アラインメントも1などの控えめな値にし、自身のcDNAアラインメントの重みは4~10とある程度高く設定したほうがBUSCOの完全性指標ではより完成度の高い遺伝子モデルが作れる印象です。つまり、ab.initioの遺伝子モデルが基準としてあって、そのモデルをより信頼性の高い(スコアの高い)cDNAやタンパク質アラインメント情報を使って塗り替えることでexon-intron境界を調整する、というイメージで重み値を決めるのが良いと思っています。どのエビデンスの重みをより高く設定するかは事前にゲノムブラウザや客観的評価手法で入念にチェックして決めます(それでも主観性は消えませんが)。
  • EVMのランは一定のチャンクサイズに分割して行われます(各コンティグごとにディレクトリができ、その中に指定したサイズと指定したサイズ分の重複領域があるサブディレクトリができる)。それでも全ての領域の計算が完了するまでには長い時間がかかります。重みやデータを変えて試行錯誤を繰り返すと、特にゲノムサイズの大きな生物ではかなりの日数を要する可能性があります(分からないけどopenMPIに対応したMAKER内でEVMをランするなら環境によっては早い?)。試行錯誤にはかなり時間がかかってしまうので、意図しないエラーがないかどうか、ランしてしばらくしたら、初期に計算される領域チャンクのサブディレクトリにある.logファイルを見ておきましょう(例えばchr1_1-100000/とか)。他の領域の計算が進行中でも、その領域の計算が終わっていればエラーが発生していないか確認できます。例えば、weightファイルにないエビデンスがあるなどのエラーが確認できるでしょう。一方でcDNAとタンパク質のエビデンスがそのチャンクで使われたかどうかは.logには書いてありません。サブディレクトリの.outファイルも開いて、cDNAとタンパク質のエビデンスが使われているかどうかも確認しましょう。
  • 重みやデータを変えて試行錯誤するとかなり時間を要します。通常、EVMのラン前にゲノムブラウザでGFF3を開いて入念にチェックしているはずです。この領域にある重みを高くしたcDNAエビデンスが使われてないのはおかしいのでは?、という感じで、計算が終わったサブディレクトリを早めにチェックしておきましょう。技術的問題が発生していれば、ジョブをkillし、修正後にやり直します。
  • 全体的にEVMの結果は良くても、ゲノムブラウザで信頼性の高いRNA seqからのエビデンスと照合していると、納得いかない遺伝子モデルも稀に見つかってくると思います(もしくはEVMで出力されてない遺伝子)。その時は該当する領域のEVMサブディレクトリにアクセスして、どんな理由でその遺伝子モデルが出力されたか、もしくは棄却されたのかを確認してみましょう。自分が遭遇したケースの1つは、CDSエビデンスだけ存在する/もしくは近縁種タンパク質のエビデンスだけ存在し、ab.initio遺伝子モデルはいずれのツールも存在しないか整合性のないモデルを報告しているという領域でした(この領域からは、EVMによって如何なる遺伝子モデルも報告されていなかった)。EVMはab.initio遺伝子モデルがベースのため、最初に正確なab.initio遺伝子モデルが構築できなかったか、CDSやタンパク質のエビデンスでもアップデート出来なかったことが原因と考えられます。
  • ab.initioの遺伝子モデルには十分に訓練されて予測された遺伝子モデルを複数使う事とが推奨されていますが、複数あると、品質が相対的に他よりに劣る、もしくは優れているモデルが見つかると思います。その時は、ab.initioの遺伝子モデル内で相対的に重み値を1くらい変えることができますが、RNA seqやESTなどのアラインメントエビデンスの重み値は超えないようにしましょう。
  • RNA seqやESTなどのアラインメントエビデンスでab.initioの遺伝子モデルを更新するならab.initioの遺伝子モデルは不要と思うかもしれませんが、そうすると発現していない条件の遺伝子は全く予測できなくなります。高品質なab.initioの遺伝子モデルを複数の手法で構築し、最高品質のアラインメントエビデンスで修正することが重要です(要確認)。

 

2023/02/15追記

1月にEVM2.0が出たとのことです。EVMのGoogle forumで著者のBrian Haasさんが報告しています。セグメント分割からコンセンサス遺伝子モデル構築、統合まで1コマンドでやってくれるようになったようです。カレントに大量のディレクトリを作ることもなくなりました。早速テストしましたが、問題なく動作しています(注;今までのバージョンとの違いは見ていません)。EVMとDockerも公式にサポートするようになったとのことです(link)。

https://groups.google.com/g/evidencemodeler-users/c/pH8NzKnOolI

(リンクされているマニュアルでは--weightが抜けてる)

 

引用

Automated eukaryotic gene structure annotation using EVidenceModeler and the Program to Assemble Spliced Alignments

Brian J Haas, Steven L Salzberg, Wei Zhu, Mihaela Pertea, Jonathan E Allen, Joshua Orvis, Owen White, C Robin Buell, Jennifer R Wortman

Genome Biol. 2008 Jan 11;9(1):R7

 

関連


参考

 *1

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

https://qiita.com/MaedaTaro_Umiushi/items/683d094144c7147f0be7

 

 

 

 

 

 

Ab initio遺伝子予測器 SNAP

 

計算機による遺伝子予測は、特に実験データの少ないゲノムに対して重要な問題であり続けている。様々なゲノムに容易に適応できるように設計されたSNAP遺伝子検出器を紹介する。また、SNAP遺伝子検出器のパラメータは、系統的に最も近いゲノムのパラメータとは限らないことを示した。また、ブートストラップパラメータ推定では、外来遺伝子探索がより有効であり、その結果得られるパラメータは非常に正確であることがわかった。

遺伝子予測は種特異的なパラメータに敏感であるため、すべてのゲノムに専用のgene finderが必要である。

 

HP

https://korflab.github.io/

Documentation

http://korflab.ucdavis.edu/software.html

 

Githubより

SNAPは、真核生物および原核生物ゲノムの両方に適した汎用遺伝子探索プログラムである。SNAPはSemi-HMM-based Nucleic Acid Parserの頭文字をとったものである。

インストール

Github

https://github.com/KorfLab/SNAP

git clone https://github.com/KorfLab/SNAP.git
cd SNAP/
make

#conda(link)
mamba install -c bioconda snap

./snap

$ ./snap 

 

SNAP - Semi-HMM-based Nucleic Acid Parser (version 2006-07-28)

 

usage: snap [options] <HMM file> <FASTA file> [options]

options:

  -help           report useful information

  -lcmask         treat lowercase as N

  -plus           predict on plus strand only

  -minus          predict on minus strand only

  -gff            output annotation as GFF

  -ace            output annotation as ACED

  -quiet          do not send progress to STDERR

  -aa <file>      create FASTA file of proteins

  -tx <file>      create FASTA file of transcripts

  -xdef <file>    external definitions

  -name <string>  name for the gene [default snap]

> ./fathom 

 

FATHOM - sequence and annotation tool (version 2006-07-28)

 

usage: fathom <ann> <dna> <commands>

commands:

  -help           report useful information

  -validate [-quiet]

  -gene-stats [-errors-ok -nucleotide -dinucleotide]

  -categorize <int>

  -export <int> [-plus -errors-ok]

  -extract <feature> -length <int> -offset <int>

  -exon-intron

  -split <-number <int> | -training <float> | -GC <float> | -repeat <float>>

  -ace-format <-gene-method <string> [-dna -extra <string>]>

  -compare-genes <predictions> [-details]

  -score-genes <hmm> [-errors-ok]

  -filter-genes <hmm> -min-score <float> -min-length <int>

 

 

テストラン

./snap HMM/thale DNA/thale.dna.gz
./snap HMM/worm DNA/worm.dna.gz

> ./snap HMM/thale DNA/thale.dna.gz

 

HMM/には代表的な生物種のHMMモデルファイルが含まれている。それでは遺伝子予測精度に問題があるならレポジトリのチュートリアルの手順を踏んで訓練しないといけない(gff3_to_zff.plを使ってZFFという形式に変換し、続いてfathomの一連のコマンドとhmm-assembler.plを走らせて3 種類のゲノムの SNAP HMM ファイルを作成する)*1。

 

 

実行方法

モデルのHMMファイルとゲノムのfastaファイルを指定する。

./snap -gff -aa output_protein.faa -tx output_cds.fna HMM/At.HMM input_genome.fna > output.gff
  • -plus         predict on plus strand only 
  • -minus     predict on minus strand only
  • -gff           output annotation as GFF
  • -ace          output annotation as ACED
  • -quiet        do not send progress to STDERR
  • -aa <file>  create FASTA file of proteins
  • -tx <file>   create FASTA file of transcripts 

 

makerのチュートリアルで紹介されているように複数回訓練させることで遺伝子予測性能は大きく向上する。

引用

Gene finding in novel genomes

Ian Korf 
BMC Bioinformatics volume 5, Article number: 59 (2004) 

 

*1

もしくはこのHPのような手順を踏む必要がある。

 

イントロン位置の保存性とRNA-seqを活用したホモロジーに基づく遺伝子予測を行う GeMoMa

明けましておめでとうございます。今年もよろしくお願いいたします。

今年も忙しくなりそうなので、更新できるタイミングがあれば積極的に更新していきます。

 

 

 GeMoMaは、進化的に関連するリファレンス種の遺伝子モデルを基に、対象種の遺伝子モデルを予測する、相同性に基づく遺伝子予測プログラムである。GeMoMaは、アミノ酸配列保存、イントロン位置保存、およびRNA-seqデータを利用して、タンパク質をコードする転写産物を正確に予測することができる。さらに、GeMoMaは複数のリファレンス種に基づく予測の組み合わせをサポートし、異なるリファレンス種の高品質アノテーションを対象種に移植することができる。ここでは、GeMoMa モジュールと GeMoMa パイプラインの詳細な説明と、特定の生物学的問題に対処するためのコマンドラインでの使用方法について紹介する。

 

Documentation

http://www.jstacs.de/index.php/GeMoMa-Docs

 

インストール

condaを使って導入した(ubuntu18)。

依存

For running the GeMoMa, you need the following software on your computer

  • Java v1.8 or later
  • blast or mmseqs

GeMoMa

http://www.jstacs.de/index.php/GeMoMa

#conda(link)
mamba create -n gemoma -y
conda activate gemoma
mamba install -c conda-forge -c bioconda gemoma=1.9 -y

GeMoMa -h

Searching for the new GeMoMa updates ...

You are using the latest GeMoMa version.

 

This jar allows to run all parts of GeneModelMapper (GeMoMa) except the external search algorithm (e.g. tblastn).

 

 

For more information please visit http://www.jstacs.de/index.php/GeMoMa

If you have any questions, comments or bugs, please check FAQs on our homepage, our github page https://github.com/Jstacs/Jstacs/labels/GeMoMa or contact jens.keilwagen@julius-kuehn.de

 

If you use this tool, please cite

 

@article{Keilwagen:2016:GeMoMa,

 author = {Keilwagen, Jens and Wenk, Michael and Erickson, Jessica L. and Schattat, Martin H. and Grau, Jan and Hartung, Frank},

 title = {{Using intron position conservation for homology-based gene prediction}},

 journal = {Nucleic Acids Research},

 volume = {44},

 number = {9},

 pages = {e89-e89},

 year = {2016},

 month = {02},

 issn = {0305-1048},

 doi = {10.1093/nar/gkw092}

}

 

@article{Keilwagen:2018:GeMoMa_RNAseq,

 author = {Keilwagen, Jens and Hartung, Frank and Paulini, Michael and Twardziok, Sven O. and Grau, Jan},

 title = {Combining RNA-seq data and homology-based gene prediction for plants, animals and fungi},

 journal = {BMC Bioinformatics},

 year = {2018},

 month = {May},

 day = {30},

 volume = {19},

 number = {1},

 pages = {189},

 issn = {1471-2105},

 doi = {10.1186/s12859-018-2203-5}

}

 

 

Available tools:

 

GeMoMaPipeline - GeMoMa pipeline

ERE - Extract RNA-seq Evidence

CheckIntrons - CheckIntrons

DenoiseIntrons - DenoiseIntrons

NRR - NCBI Reference Retriever

Extractor - Extractor

GeMoMa - GeneModelMapper

GAF - GeMoMa Annotation Filter

AnnotationFinalizer - AnnotationFinalizer

AnnotationEvidence - Annotation evidence

Attribute2Table - Attribute2Table

SyntenyChecker - Synteny checker

AddAttribute - AddAttribute

GAFComparison - GAFComparison

Analyzer - Analyzer

BUSCORecomputer - BUSCORecomputer

GFFAttributes - GFFAttributes

TranscribedCluster - Transcribed Cluster

 

Syntax: GeMoMa <toolname> [<parameter=value> ...]

 

Further info about the tools is given with

GeMoMa <toolname> info

 

For tests of individual tools:

GeMoMa <toolname> test [<verbose>]

 

Tool parameters are listed with

GeMoMa <toolname>

 

GeMoMa ERE

$ GeMoMa ERE 

Searching for the new GeMoMa updates ...

You are using the latest GeMoMa version.

 

Parameters of tool "Extract RNA-seq Evidence" (ERE, version: 1.9):

s - Stranded (Defines whether the reads are stranded. In case of FR_FIRST_STRAND, the first read of a read pair or the only read in case of single-end data is assumed to be located on forward strand of the cDNA, i.e., reverse to the mRNA orientation. If you are using Illumina TruSeq you should use FR_FIRST_STRAND., range={FR_UNSTRANDED, FR_FIRST_STRAND, FR_SECOND_STRAND}, default = FR_UNSTRANDED) = FR_UNSTRANDED

The following parameter(s) can be used multiple times:

m - mapped reads file (BAM/SAM files containing the mapped reads, type = bam,sam) = null

v - ValidationStringency (Defines how strict to be when reading a SAM or BAM, beyond bare minimum validation., range={STRICT, LENIENT, SILENT}, default = LENIENT) = LENIENT

u - use secondary alignments (allows to filter flags in the SAM or BAM, default = true) = true

c - coverage (allows to output the coverage, default = true) = true

mmq - minimum mapping quality (reads with a mapping quality that is lower than this value will be ignored, valid range = [0, 255], default = 40) = 40

mc - minimum context (only introns that have evidence of at least one split read with a minimal M (=(mis)match) stretch in the cigar string larger than or equal to this value will be used, valid range = [1, 1000000], default = 1) = 1

maximumcoverage - maximum coverage (optional parameter to reduce the size of coverage output files, coverage higher than this value will be reported as this value, valid range = [1, 10000], OPTIONAL) = null

f - filter by intron mismatches (filter reads by the number of mismatches around splits, range={NO, YES}, default = NO) = NO

    No parameters for selection "NO"

    Parameters for selection "YES":

    r - region around introns (test region of this size around introns/splits for mismatches to the genome, valid range = [0, 100], default = 10) = 10

    n - number of mismatches (number of mismatches allowed in regions around introns/splits, valid range = [0, 100], default = 3) = 3

    t - target genome (The target genome file (FASTA). Should be in IUPAC code, type = fasta,fas,fa,fna,fasta.gz,fas.gz,fa.gz,fna.gz) = null

e - evidence long splits (require introns to have at least this number of times the supporting reads as their length deviates from the mean split length, valid range = [0.0, 100.0], default = 0.0) = 0.0

mil - minimum intron length (introns shorter than the minimum length are discarded and considered as contiguous, valid range = [0, 1000], default = 0) = 0

repositioning - repositioning (due to limitations in BAM/SAM format huge chromosomes need to be split before mapping. This parameter allows to undo the split mapping to real chromosomes and coordinates. The repositioning file has 3 columns: split_chr_name, original_chr_name, offset_in_original_chr, type = tabular, OPTIONAL) = null

outdir - The output directory, defaults to the current working directory (.) = .

 

 

実行方法

v1.9では以下のコマンドを利用できる。

  • GeMoMaPipeline - GeMoMa pipeline
  • ERE - Extract RNA-seq Evidence
  • CheckIntrons - CheckIntrons
  • DenoiseIntrons - DenoiseIntrons
  • NRR - NCBI Reference Retriever
  • Extractor - Extractor
  • GeMoMa - GeneModelMapper
  • GAF - GeMoMa Annotation Filter
  • AnnotationFinalizer - AnnotationFinalizer
  • AnnotationEvidence - Annotation evidence
  • Attribute2Table - Attribute2Table
  • SyntenyChecker - Synteny checker
  • AddAttribute - AddAttribute
  • GAFComparison - GAFComparison
  • Analyzer - Analyzer
  • BUSCORecomputer - BUSCORecomputer
  • GFFAttributes - GFFAttributes
  • TranscribedCluster - Transcribed Cluster

 

 

GeMoMaPipelineコマンド

このコマンドはGeMoMaパイプライン全体を実行する。基本的には以下のように実行される;Extract RNA-seq evidence (ERE), DenoiseIntrons, Extractor, external search (tblastn or mmseqs), Gene Model Mapper (GeMoMa), GeMoMa Annotation Filter (GAF), AnnnotationFinalizerが基本的に動作する。マルチスレッドで、1台のマシンのすべての計算コアを利用することができるが、計算クラスタのように分散させることはできない。

アノテーションをつける自身のゲノム(t)、リファレンス種1の名前(i)とアノテーション(a)と配列(g)、外部エビデンスであるRNAseqのbamファイル(ERE.m)(r=MAPPEDの時のみ)、スレッド数(threads)、出力ディレクトリ名(outdir)を指定する。

GeMoMa GeMoMaPipeline t=genome.fa r=NO o=true \
i=<reference_1_id> a=<reference_1_annotation> g=<reference_1_genome> \
GeMoMa.Score=ReAlign AnnotationFinalizer.r=NO \
threads=20 outdir=outdir

#RNA seqのbamも指定(HISAT2などでmappingして得たbam)
GeMoMa GeMoMaPipeline t=genome.fa r=MAPPED o=true \
i=<reference_1_id> a=<reference_1_annotation> g=<reference_1_genome> \
GeMoMa.Score=ReAlign AnnotationFinalizer.r=NO \
ERE.m=map.bam \
threads=20 outdir=outdir
  • t    target genome (Target genome file (FASTA), type = fasta,fa,fas,fna,fasta.gz,fa.gz,fas.gz,fna.gz)
  • r    RNA-seq evidence (data for RNA-seq evidence, range={NO, MAPPED, EXTRACTED}, default = NO)
  • o    output individual predictions (If *true*, returns the predictions for each reference species, default = false)
  • s    species (data for reference species, range={own, pre-extracted}, default = own)
  • i    ID (ID to distinguish the different reference species, OPTIONAL)
  • a    annotation (Reference annotation file (GFF or GTF), which contains gene models annotated in the reference genome, type = gff,gff3,gtf,gff.gz,gff3.gz,gtf.gz)    FILE
  • g    genome (Reference genome file (FASTA), type = fasta,fa,fas,fna,fasta.gz,fa.gz,fas.gz,fna.gz)
  • outdir    The output directory, defaults to the current working directory (.)    STRING
  • ERE.m    Parameters for selection "MAPPED". mapped reads file (BAM/SAM files containing the mapped reads, type = bam,sam) 
  • threads    The number of threads used for the tool, defaults to 1
  • GeMoMa.Score    Score (A flag which allows to do nothing, re-score or re-align the search results, range={Trust, ReScore, ReAlign}, default = ReAlign)
  • AnnotationFinalizer.r    rename (allows to generate generic gene and transcripts names (cf. parameter "name attribute"), range={COMPOSED, SIMPLE, NO}, default = COMPOSED)
  • GeMoMa.m    maximum intron length (The maximum length of an intron, default = 15000)
  • GeMoMa.sil    static intron length (A flag which allows to switch between static intron length, which can be specified by the user and is identical for all genes, and dynamic intron length, which is based on the gene-specific maximum intron length in the reference organism plus the user given maximum intron length, default = true)

Parameters for selection "DENOISE":

  • DenoiseIntrons.m     maximum intron length (The maximum length of an intron, default = 15000)    INT
  • DenoiseIntrons.me     minimum expression (The threshold for removing introns, valid range = [0.0, 1.0], default = 0.01)    DOUBLE
  • DenoiseIntrons.c     context (The context upstream a donor and donwstream an acceptor site that is used to determine the expression of the region, valid range = [0, 100], default = 10)

複数のリファレンスがある場合は、sとパラメータタグ i, a, g を対応する値で繰り返す。o=trueにすると個々の予測が別ファイルとして出力され、個々のステップ簡単かつ迅速に再実行できる。最大イントロン長を指定したい場合は、GeMoMa.m と GeMoMa.sil パラメータを検討する。RNA-seqデータがある場合、DenoiseIntrons のパラメータも確認する。

出力例

 

 

特定のコマンドのみ実行することもできます。EREコマンドの例を示す。

EREコマンド

マップされたRNA-seqリードからイントロンカバレッジを抽出する。イントロンはDenoiseIntronsコマンドでノイズ除去できる可能性がある。EREの結果は予測を改善するために使用することができ、GAFコマンドでより良い遺伝子モデルを選択するのに役立つ可能性がある。また、EREの結果は、AnnotationFinalizerコマンドによるUTRの予測に使用することができる(2018年の論文のFig.1参照)。

(STARやHISAT2などで得た)sam/bamファイルを指定する。s以外のパラメータは複数回指定できる。

GeMoMa ERE s=FR_UNSTRANDED m=input.bam outdir=OUTPUT
  • m    mapped reads file (BAM/SAM files containing the mapped reads, type = bam,sam)
  • s    Stranded (Defines whether the reads are stranded. In case of FR_FIRST_STRAND, the first read of a read pair or the only read in case of single-end data is assumed to be located on forward strand of the cDNA, i.e., reverse to the mRNA orientation. If you are using Illumina TruSeq you should use FR_FIRST_STRAND., range={FR_UNSTRANDED, FR_FIRST_STRAND, FR_SECOND_STRAND}, default = FR_UNSTRANDED)
  • mil    minimum intron length (introns shorter than the minimum length are discarded and considered as contiguous, valid range = [0, 1000], default = 0)
  • mmq    minimum mapping quality (reads with a mapping quality that is lower than this value will be ignored, valid range = [0, 255], default = 40)
  • v    ValidationStringency (Defines how strict to be when reading a SAM or BAM, beyond bare minimum validation., range={STRICT, LENIENT, SILENT}, default = LENIENT)    STRING
  • u    use secondary alignments (allows to filter flags in the SAM or BAM, default = true)
  • repositioning    repositioning (due to limitations in BAM/SAM format huge chromosomes need to be split before mapping. This parameter allows to undo the split mapping to real chromosomes and coordinates. The repositioning file has 3 columns: split_chr_name, original_chr_name, offset_in_original_chr, type = tabular, OPTIONAL)

出力例

それぞれのコマンドについてはDOcumentationで説明されています。確認してください。

引用

GeMoMa: Homology-Based Gene Prediction Utilizing Intron Position Conservation and RNA-seq Data

Jens Keilwagen, Frank Hartung, Jan Grau

Methods Mol Biol. 2019;1962:161-177

 

Combining RNA-seq data and homology-based gene prediction for plants, animals and fungi

Jens Keilwagen, Frank Hartung, Michael Paulini, Sven O. Twardziok & Jan Grau 
BMC Bioinformatics volume 19, Article number: 189 (2018)

 

Using intron position conservation for homology-based gene prediction

Jens Keilwagen, Michael Wenk, Jessica L Erickson, Martin H Schattat, Jan Grau, Frank Hartung

Nucleic Acids Res. 2016 May 19;44(9):e89

 

参考

GeMoMa で近縁種のアノテーションを移植する

https://qiita.com/drk0311/items/32aa45b85e0ad3545067

 

ProtHint

 

著者らは、真核生物ゲノムにおける遺伝子予測のための高速かつ高精度なアルゴリズムの作成に向けて、いくつかのステップを踏んできた。まず、ab initio gene findingを効率的に行うための自動化手法、GeneMark-ESを導入し、教師無しでパラメータを繰り返し学習させた。次に、GeneMark-ETでは、教師なし学習とShort RNA Readsのマッピングによって明らかになったイントロン位置の情報を統合する方法を提案した。GeneMark-EPは、タンパク質データベースという、シーケンシングプロジェクト開始前に入手可能な別の外部情報源を利用するツールである。GeneMark-EPは、ProtHintという新しいパイプラインにより、ゲノムに大量のタンパク質マッピングを行い、遺伝子候補のスプライスサイトや翻訳開始・停止サイトに関するヒントを抽出する。GeneMark-EPは、このヒントを用いて、モデルパラメータの推定を改善するとともに、予測される遺伝子の座標が最も信頼できるヒントと異なる場合(-EP+モード)、その座標を調整することができる。GeneMark-EP および -EP+ のテストでは、GeneMark-ES と比較して遺伝子予測精度の向上が見られ、GeneMark-EP+ は GeneMark-ET よりも高い精度を示した。遺伝子予測精度が最も顕著に向上したのは、大規模な真核生物ゲノムにおいてであることが確認された。

 

Githubより

ProtHintは、予測された遺伝子を参照タンパク質配列のデータベースにマッピングし、スプライシングアラインメントすることによって、目的のゲノム中のヒント(イントロン、開始コドン、停止コドンの形)を予測し、スコアリングするパイプラインである。

 

注;ProtHintはGeneMark-ESで使用されていて、GeneMark-ESの配布物自体にも含まれている。ProtHint pipelineの流れはGeneMark-EP+の論文のFigure.2で説明されている。

インストール

依存

The following non-Core Perl modules are required:

  • MCE::Mutex
  • threads
  • YAML
  • Math::Utils

Github

#perlライブラリの導入
cpanm MCE::Mutex threads YAML Thread::Queue Math::Utils

#GeneMark-ES
cd GeneMark_dir
#Install the key
zcat gm_key_64.gz > ~/.gm_key
#test run 
bash check_install.bash
#パスを通す
export PATH=$PATH:/your/path/to/GeneMark-ES/

#本体
git clone https://github.com/gatech-genemark/ProtHint.git
cd ProtHint/

> ./prothint.py -h

usage: prothint.py [-h] [--workdir WORKDIR] [--geneSeeds GENESEEDS] [--geneMarkGtf GENEMARKGTF] [--fungus] [--nonCanonicalSpliceSites] [--diamondPairs DIAMONDPAIRS] [--maxProteinsPerSeed MAXPROTEINSPERSEED] [--evalue EVALUE]

                   [--minExonScore MINEXONSCORE] [--longGene LONGGENE] [--longProtein LONGPROTEIN] [--cleanup] [--pbs] [--threads THREADS] [--version] [--prevGeneSeeds PREVGENESEEDS] [--prevSpalnGff PREVSPALNGFF]

                   genome.fasta proteins.fasta

 

ProtHint 2.6.0: Pipeline for generating genome wide footprints of homologous proteins. The set of high confidence hints is generated using default thresholds in print_high_confidence.py script. If you wish to use different filtering criteria, re-run

print_high_confidence.py script with custom thresholds.

 

positional arguments:

  genome.fasta          Input genomic sequence in FASTA format.

  proteins.fasta        Protein database in FASTA format.

 

optional arguments:

  -h, --help            show this help message and exit

  --workdir WORKDIR     Folder for results and temporary files. If not specified, current directory is used

  --geneSeeds GENESEEDS

                        Gene prediction seeds in gtf format. If this file is provided, GeneMark-ES run is skipped.

  --geneMarkGtf GENEMARKGTF

                        Same as --geneSeeds, for backwards compatibility

  --fungus              Run GeneMark-ES in fungus mode.

  --nonCanonicalSpliceSites

                        Predict introns with non-canonical splices sites. By default, ProtHint only predicts introns with GT-AG, GC-AG, and AT-AC donor-acceptor sites.

  --diamondPairs DIAMONDPAIRS

                        File with "seed gene-protein" hits generated by DIAMOND. If this file is provided, DIAMOND search for protein hits is skipped. The seed genes in this file must correspond to seed genes passed by "--geneSeeds" option. All pairs in

                        the file are used -- option "--maxProteinsPerSeed" is ignored.

  --maxProteinsPerSeed MAXPROTEINSPERSEED

                        Maximum number of protein hits per seed gene. Increasing this number leads to increased runtime and may improve the sensitivity of hints. Decreasing has an opposite effect. Default is set to 25.

  --evalue EVALUE       Maximum e-value for DIAMOND alignments hits. Default = 0.001

  --minExonScore MINEXONSCORE

                        Discard all hints inside/neighboring exons with score lower than minExonScore. Default = 25

  --longGene LONGGENE   Threshold for what is considered a long gene in Spaln alignment. Default = 30000

  --longProtein LONGPROTEIN

                        Threshold for what is considered a long protein in Spaln alignment. Default = 15000

  --cleanup             Delete temporary files and intermediate results. Cleanup is turned off by default as it is useful to keep these files for troubleshooting and the intermediate results might be useful on their own.

  --pbs                 Run GeneMark-ES and Spaln on pbs.

  --threads THREADS     Number of threads used by ES, DIAMOND, and Spaln. By default, all available threads are used.

  --version             show program's version number and exit

  --prevGeneSeeds PREVGENESEEDS

                        File with gene seeds which were used in the previous iteration. Next iteration of ProtHint is only executed for --geneSeeds which differ from --prevGeneSeeds. --prevSpalnGff is required when this option is used since results from

                        the previous iteration are reused for seeds which do not differ (Gene ids of such hints are updated to match the new seed genes).

  --prevSpalnGff PREVSPALNGFF

                        Scored hints from previous iteration.

 

gmes_petap.pl

$ gmes_petap.pl 

# -------------------

Usage:  /home/kazu/Documents/ProtHint/gmes_linux_64/gmes_petap.pl  [options]  --sequence [filename]

 

GeneMark-ES Suite version 4.71_lic

Suite includes GeneMark.hmm, GeneMark-ES, GeneMark-ET and GeneMark-EP algorithms.

 

Input sequence/s should be in FASTA format.

 

Select one of the gene prediction algorithms:

 

  To run GeneMark-ES self-training algorithm

    --ES

 

  To run GeneMark-ET with hints from transcriptome splice alignments

    --ET           [filename]; file with intron coordinates from RNA-Seq read splice alignment in GFF format

    --et_score     [number]; default 10; minimum score of intron in initiation of the ET algorithm

 

  To run GeneMark-EP with hints from protein splice alignments

    --EP           

    --dbep         [filename]; file with protein database in FASTA format

    --ep_score     [number,number]; default 4,0.25; minimum score of intron in initiation of the EP algorithm

    or

    --EP           [filename]; file with intron coordinates from protein splice alignment in GFF format

 

  To run GeneMark.hmm predictions using previously derived model

    --predict_with [filename]; file with species specific gene prediction parameters

 

  To run ES, ET or EP with branch point model. This option is most useful for fungal genomes

    --fungus

 

  To run hmm, ES, ET or EP in PLUS mode (prediction with hints)

    --evidence     [filename]; file with hints in GFF format

 

  To run algorithms with alternative genetic codes

    --gcode      [number]; default 1; supported 1 and 6/26

 

Output formatting options:

  --format       [label]; default GTF; output gene prediction in GTF of GFF3 format

  --work_dir     [folder name]; default current working directory .;

 

Masking option

  --soft_mask    [number] or [auto]; default auto; to indicate that lowercase letters stand for repeats;

                 algorithm hard masks only lowercase repeats longer than specified length

                 In 'auto' mode hard masking threshold is selected by algorithm based on the size of the input genome

  --mask_penalty [number] or [auto]; default 0.03;

 

Run options

  --cores        [number]; default 1; to run program with multiple threads

  --pbs          to run on cluster with PBS support

  --v            verbose

 

Optional sequence pre-processing parameters

  --max_contig   [number]; default 5000000; will split input genomic sequence into contigs shorter then max_contig

  --min_contig   [number]; default 50000 (10000 fungi); 

                 will ignore contigs shorter than min_contig in training 

  --max_gap      [number]; default 5000; will split sequence at gaps longer than max_gap

                 Letters 'n' and 'N' are interpreted as standing within gaps 

  --max_mask     [number]; default 5000; will split sequence at repeats longer then max_mask

                 Letters 'x' and 'X' are interpreted as results of hard masking of repeats

 

Optional parameters

  --max_intron            [number]; default 10000 (3000 fungi); maximum length of intron

  --max_intergenic        [number]; default 50000 (10000 fungi); maximum length of intergenic regions

  --min_contig_in_predict [number]; default 500; minimum allowed length of contig in prediction step

  --min_gene_in_predict   [number]; default 300 (120 fungi); minimum allowed gene length in prediction step

  --gc_donor              [value];  default 0.001; transition probability to GC donor in the range 0..1; 

                          'off' switches GC donor model OFF

  --gc3          [number]; GC3 cutoff in training for grasses

 

Developer options

  --training     to run only training step of algorithm; applicable to ES, ET or EP

  --prediction   to run only prediction step of algorithms using species parameters from previously executed training; applicable to ES, ET or EP

  --usr_cfg      [filename]; use custom configuration from this file

  --ini_mod      [filename]; use this file with parameters for algorithm initiation

  --key_bin

  --debug

# -------------------

 

 

データベース

リファレンスタンパク質配列のソースとして、OrthoDBタンパク質データベースの中で対象のゲノムに関連する分類を使用することが推奨されている。

#節足動物門の場合(ほかのリンクはGithub参照
wget https://v100.orthodb.org/download/odb10_arthropoda_fasta.tar.gz
tar xvf odb10_arthropoda_fasta.tar.gz
cat arthropoda/Rawdata/* > arthropoda_proteins.fasta

#緑色植物亜界
wget https://v100.orthodb.org/download/odb10_plants_fasta.tar.gz
tar xvf odb10_plants_fasta.tar.gz
cat plants/Rawdata/* > plants_proteins.fasta

 

 

テストラン

ゲノムのfastaファイル、proteomeファイル(上で作ったデータベースのfastaや信頼できる近縁種のproteome)を指定する。--geneSeedsをつけてGeneMark-ESからのgtfファイルも指定するとGeneMark-ESによるab initio遺伝子予測はスキップされる。このフラグを付けてないなら、初めにGeneMark-ESが実行される(パスが通っている必要がある)。結果はworkdirに保存する。

git clone https://github.com/gatech-genemark/ProtHint.git
cd ProtHint/example/
../bin/prothint.py input/genome.fasta input/proteins.fasta --geneSeeds input/genemark.gtf --workdir workdir


#plantのproteins.fastaをDBに使う。訓練も行う。
../bin/prothint.py input/genome.fasta plants_proteins.fasta --workdir workdir
  • --geneSeeds    Gene prediction seeds in gtf format. If this file is provided, GeneMark-ES run is skipped.

 

ProtHintは主に2つの出力を生成する(マニュアルより)。

  • prothint.gff   報告されたすべてのヒント(イントロン、スタート、ストップ)を含む Gff ファイル。
  • evidence.gff    Prothint.gff の高信頼度サブセットで、例えば GeneMark-EP Plus モードに適している。このセットは、ProtHint/bin/print_high_confidence.pyスクリプトのデフォルトの閾値を使用して生成されている。異なるフィルタリング基準を使用したい場合は、カスタムしきい値でprint_high_confidence.pyスクリプトを再実行する。
  • prothint_augustus.gff    BRAKERとAUGUSTUSですぐに使える出力も生成する。

出力例

workdir/

 

続いてGeneMark-ESかEP+を実行する(Github)。

gmes_petap.pl --EP workdir/prothint.gff --evidence workdir/evidence.gff --seq genome.masked.fasta --soft_mask 1000 --verbose --cores 20 --work_dir workdir
  • --EP    file with intron coordinates from protein splice alignment in GFF format
  • --ES    To run GeneMark-ES self-training algorithm
  • --evidence    file with hints in GFF format
  • --soft_mask [number] or [auto]   default auto; to indicate that lowercase letters stand for repeats;  algorithm hard masks only lowercase repeats longer than specified length  In 'auto' mode hard masking threshold is selected by algorithm based on the size of the input genome
  • --max_contig   default 5000000; will split input genomic sequence into contigs shorter then max_contig
  • --min_contig   default 50000 (10000 fungi); will ignore contigs shorter than min_contig in training 
  • --max_gap   default 5000; will split sequence at gaps longer than max_gap Letters 'n' and 'N' are interpreted as standing within gaps 
  • --max_mask    default 5000; will split sequence at repeats longer then max_mask Letters 'x' and 'X' are interpreted as results of hard masking of repeats
  • --cores    default 1; to run program with multiple threads
  • --format    default GTF; output gene prediction in GTF of GFF3 format
  • --work_dir    default current working directory .;

出力例

 

引用

https://github.com/gatech-genemark/ProtHint

 

GeneMark-EP+: eukaryotic gene prediction with self-training in the space of genes and proteins
Tomáš Brůna, Alexandre Lomsadze, Mark Borodovsky

NAR Genom Bioinform. 2020 Jun;2(2):lqaa026. doi: 10.1093/nargab/lqaa026. Epub 2020 May 13.

 

今年も大変お世話になりました。来年もよろしくお願いします。

 

AUGUSTUSの訓練と遺伝子予測のためのウェブサービス WebAUGUSTUS

2023/01/01 誤字修正

 

タンパク質をコードする遺伝子の予測は、新たにシークエンシングされアセンブルされたゲノムのアノテーションにおいて重要なステップである。AUGUSTUSは真核生物の遺伝子予測のための最も正確なツールの一つである。ここでは、AUGUSTUSの学習とAUGUSTUSによる遺伝子予測のためのウェブインタフェース、WebAUGUSTUSを紹介する。WebAUGUSTUS はユーザーのニーズに応じて、トレーニング用遺伝子構造を自動的に生成する。このステップでは、ゲノムファイルの他に、ESTタグのファイルまたはタンパク質配列のファイルのいずれかが必要である。また、外部で作成したトレーニング用遺伝子構造ファイルとゲノムファイルを提出することもできる。Webサービスでは、AUGUSTUSのパラメータを最適化し、そのパラメータで遺伝子を予測する。WebAUGUSTUS は http://bioinf.uni-greifswald.de/webaugustus で利用できる。

Tutorial

http://bioinf.uni-greifswald.de/webaugustus/predictiontutorial#param_id

 

webサービス

Bioinformatics Web Server - University of Greifswald

一番下のAUGUSTUS training submissionかAUGUSTUS prediction submissionをクリックする。旧AUGUSTUSウェブサーバー(Pubmed1, 2)へのリンクも残っている。こちらも遺伝子予測サービスを提供しているがトレーニングには対応していない。

 

 

training submission(*1)

新規ゲノムデータ用のAUGUSTUSパラメータをトレーニングするためのデータをサブミットする。パラメータが既にトレーニングされ、生物種概要の表で公開されているかどうかを確認する(link)。個人ヒトゲノムアセンブリは許可されていないことに注意。

メールアドレス、種名と(multi-)fasta形式のゲノムファイル(短くてユニークな、ホワイトスペースを含まないfastaヘッダー、DNA配列は80 文字の後に unix 形式で改行される形式)を指定する。最大100MB、250,000 scaffoldとなっている(gzip圧縮可)。cDNAファイル(*2)が提供された場合、まずこのcDNAファイルからヒントファイルが生成される。また、遺伝子構造の証拠を含むヒントファイルを gff 形式で提出することが可能。ヒントファイルは、遺伝子構造予測をサポートする外部証拠として使用される。ヒントはアライメントプログラムや情報資源(例:EST、RNA-seqデータ、ペプチド、タンパク質、...)から、適当と思われるものを自分で生成することができる。

cDNAファイルかヒントファイル、もしくは近縁種proteinどれか1つは提供する必要がある。最適な結果を得るためには、異なるファイルの組み合わせで何度もAUGUSTUSのトレーニングを開始するのではなく、1回のトレーニング実行で可能な限り多くの情報を指定する必要がある。ESTデータを所有している場合は、タンパク質配列ではなくESTの使用が推奨される。ESTを使うことで、UTRトレーニングセットを生成できる可能性が高いため。(マニュアルより)。

コメント;訓練にソフトマスクされた配列は使用しない。エラーになる。

 

 

AUGUSTUS prediction submission

リスト(link)にない生物は、(上記)トレーニング後に実行する。

ウェブサーバーで AUGUSTUS をトレーニングした場合、ジョブ確認メールに記載されているトレーニング実行のプロジェクト識別子(trainまたはpredに続く8桁の数字)を指定するだけでよい。試した時は新規生物種だったため、トレーニングを行い、その結果のparameters.tar.gzを指定した。また、その後のEVMのvalidaterスクリプトで重複フィーチャーエラーが起きたので、Allowed genee structuredの項目ではexactly one gene を選択した。

 

  • UTR予測について

UTR予測は、使用する生物種のUTRパラメータファイルが存在する場合にのみ可能。ある種のUTRパラメータファイルが存在する場合でも、それが種特異的であること、つまり、ターゲットとする種に対して実際に最適化されていることを確認する必要がある。一般的な(テンプレート)パラメータでUTRを予測するのは時間の無駄である。

UTRパラメータファイルが存在しないにもかかわらず、フォームでUTR予測を有効にした場合、WebサーバーアプリケーションはUTRを予測するという選択を、単にどのUTRも予測しないことで覆す。

デフォルトでは、AUGUSTUS はオルタナティブトランスクリプトを予測しない。

fewを選択すると、AUGUSTUSパラメータが設定され、比較的少数の代替転写産物が予測されるようになる。

  •  ストランド特異的予測

ストランドチェックボックスにチェックをつけると、AUGUSTUSは両鎖の遺伝子を予測する。チェックによって、順鎖(+)または逆鎖(-)の遺伝子のみを予測するようにもできる。

(マニュアルより)

 

コメント

空いている時はすぐにジョブが実行されますが、混雑している時もあります。例えば10日ほどジョブ待ちだった時もありますが、すぐにcomputingのステータスに入って数時間で結果が得られた時もありました(gzip圧縮で50MBのゲノム)。

引用

WebAUGUSTUS—a web service for training AUGUSTUS and predicting genes in eukaryotes 
Katharina J. Hoff, Mario Stanke
Nucleic Acids Research, Volume 41, Issue W1, 1 July 2013, Pages W123–W128

 

*1

AUGUSTUSはAb initio 遺伝子予測を行う。つまり、外部エビデンス(ESTやタンパク質アライメントなど)ではなく、コドン頻度やイントロン-エクソン長分布のような生物特有のゲノム情報をもとに数学的モデルを用いて遺伝子を同定し、そのイントロン-エクソン構造を決定する。 外部エビデンスを要求するEvidence-driven gene predictionな実装が低発現の遺伝子の予測を行えないことが多いのに対し、Ab initio 遺伝子予測は、発現していない遺伝子の予測も行えるのがメリットである。しかし新しいゲノムに適用するには、近縁な種のproteomeやcDNAで事前訓練されたパラメータファイルが必要になるのはデメリットになる。この訓練には非常に近縁な(closely related)生物種のものを使わないといけないが、近縁種であっても、エキソンーイントロン構造やコドン使用率が異なる事がある。よって対象の種自身の外部エビデンスを使って訓練するのが望ましい(このWebAUGUSTUSなら自身のゲノムからのcDNA(EST)を提供して訓練する)。AUGUSTUSに外部エビデンスは必須ではないが、このような背景から、外部エビデンスはあったほうが有利になる可能性はある。

Ab initio 遺伝子予測の実装は、Evidence-driven gene predictionとは異なり、基本的に最も可能性の高いコーディング配列を1つだけ見つけ、非翻訳領域(UTR)やオルタナティブスプライシングは報告しない(AUGUSTUSはUTR予測も任意で行えるが完全ではない)。その点を考慮し、Ab initio 遺伝子予測結果はEvidence-driven gene predictionな方法の結果と組み合わせ、精度の高い遺伝子モデルを構築していくワークフローの1ステップとして使う。

もう1つ忘れてはならないのは、基本的にゲノムサイズが大きいほど1つ1つの遺伝子サイズが大きくなるという傾向があることである(引用)(例;ヒトには多くのイントロンで分断された非常に長い遺伝子がある)。この傾向から、ゲノムサイズが大きいほど事前により連続性の高いゲノムアセンブリを構築しておかないと、より断片化された遺伝子が予測されるリスクが上がる傾向がある。ゲノムサイズが大きいほど様々な種類のリピート配列が増えアセンブルは難しくなる傾向があるので(N50などの基準で連続性が低くなる)、ゲノムサイズが大きい種ほど、遺伝子予測も難しくなるという悪循環が発生しやすい。一方で、ゲノムサイズが小さく遺伝子密度が高い種では、遺伝子予測の方法次第で隣接した遺伝子間が融合してしまった遺伝子モデルが出来やすい別のリスクもある。

品質の低い遺伝子モデルの氾濫は公共データベースにとって大きな問題である。遺伝子予測の難しさを研究者が危惧してかどうかは分からないが、3大データベースには、ゲノムアセンブリのみが登録され、アノテーション結果は登録されていない真核生物種のゲノムデータがたくさん見つかる(主観では半分以上)。生物学の研究ではゲノムを決める動機はどんな遺伝子があるか調べることにあり、それを放棄しているのは奇妙に聞こえるかもしれないが、そのようなゲノムプロジェクトも、1)自前のHPでのみアノテーションを公開しているか、2)水面下で時間をかけて現在も遺伝子予測を実行中であるか(キュレーションすると長い時間がかかる)、3)遺伝子予測結果はオフラインでラボメンのみが使用している、などの可能性がある。

 

*2

気軽に使用できる例として、近縁のモデル生物のcDNAファイルを上げることができるが、あらゆる形態のエビデンスの中で、対象それ自身の種のRNA-seq データが
遺伝子アノテーションの精度を向上させる可能性が最も高い。short readのRNA seqをヒントに使うなら、RNA seqのリードをTrinityなどでアセンブルし、TransDecorderでタンパク質をコードしている配列を探索し(raw assemblyにはnon-codingやstart-endを含まない断片化されたcDNA配列がある)、それを提供することが度々ある。

 

 

シークエンシングリードからk-merカウントを使用してテロメアリピートを予測する telomere-kmer-search

 

レポジトリより

telomere-kmer-searchは、テロメアが与えられたゲノムライブラリの中で最も豊富なダイレクトタンデムリピート配列であると仮定してテロメアリピートを予測する。出力は、頻度順に並べたダイレクトタンデムリピートユニット配列のリストである。
k-merは、まずk-mer配列となる前にjellyfish countでカウントされ、そのカウント値がFasta形式にダンプされる。その後、最小カバレッジ以上のものは保持され、最も豊富なタンデムダイレクトリピートの探索に使用される。このタスクのためのbashスクリプトの例がdo_jellyfish.shに示されている。これらのスクリプトjellyfish v2.2.10でテストされている。依存関係はCondaの環境ファイルenv.ymlに記載されている。


インストール

Github

git clone https://github.com/Swart-lab/telomere-kmer-search.git
cd telomere-kmer-search/
#依存ツール(ymlのchannnel名をnullからここではtelomere-kmer-searchに修正した)
mamba env create --file env.yml
conda activate telomere-kmer-search

> ./find_repeats_from_kmers.py -h

usage: Identify direct tandem repeat elements from k-mer counts [-h] [--counts COUNTS] [--output OUTPUT] [--maxzeroes MAXZEROES] [-k K]

 

options:

  -h, --help            show this help message and exit

  --counts COUNTS, -c COUNTS

                        Dict of counts per k-mer, in JSON format

  --output OUTPUT, -o OUTPUT

                        Prefix for output files

  --maxzeroes MAXZEROES

                        Maximum number zero-count expected k-mers to allow in reporting repeats

  -k K                  K-mer length used for counting, should be prime

 

 

 

実行方法

jellyfish でカウント後、シングルトンをフィルタリングして高頻度k-merを含むJSONファイルをfind_repeats_from_kmersに渡す。

jellyfish count -m 19 -s 10000000 -t 20 -C input.fastq -o mer_counts.jf
jellyfish dump count.jf > dumps.fa
./discard_singletons.py --input dumps.fa --cutoff 100 --histo hist.json --output cov.json
./find_repeats_from_kmers.py --maxzeroes 1 -k 19 \
--counts cov.json --output library.k19.repeats

出力例

レポジトリより

  • k-merのカウントに使用するk-merの長さは、予想される最大テロメア繰り返し単位長より長く、プライムである必要がある。
  • discard_singletons.pyに指定する最小カバレッジは、ゲノムの平均カバレッジよりやや高めで、繰り返し要素のみを保持するようにする必要がある。
  • 高頻度k-merのカウントを含むJSONファイルは、jellyfish countに使われるk-merの長さとともに、find_repeats_from_kmers.pyスクリプトに供給される。

 

引用

https://github.com/Swart-lab/telomere-kmer-search

 

関連


 

gff3出力をサポートを追加したexonerateのフォーク exonerate-gff3

2023/01/05 追記

2023 01/13 パラメータの解釈間違いを修正

 

Exonerateはペアワイズ配列比較のためのツール。DNAとcDNA(EST)、DNAとタンパク質間のアライメントを行うことができる。アライメントモデルに基づき、ギャップありアライメント、ギャップなしアライメント行うことができる。類似したソフトウエアには、Wise2Genewise)やUCSCBLAT(BLAST-like alignment tool)、最近開発されたMiniprotなどがある。これらは、開発された時期が大きく異なり、全ゲノムにスケールするか、精度(例;イントロンーエキソン境界の予測精度)などに違いがある。遅いツールを使ってゲノム全体に全タンパク質をアラインする場合、近似マッピングを求め、それから正確なアラインメントを行うなどの工夫が必要になる。

ここでは、gff3出力に対応したexonerateのフォークを紹介します。これを使うことで、スクリプトを書かなくてもexonerateから直接GFF3形式のアノテーションを得ることができます。

 

EMBL-EBIのexonerateのmanual

https://www.ebi.ac.uk/about/vertebrate-genomics/software/exonerate-manual

 

インストール

リリースからexonerate-2.3.0-x86_64.tar.gzをダウンロードしてbin/にパスを通した(ubuntu18)。

Github

git clone https://github.com/hotdogee/exonerate-gff3.git
cd exonerate-gff3/
./configure
make -j8
make install

exonerate

exonerate from exonerate version 2.3.0

Using glib version 2.26.1

Built on Jul  1 2014

Branch: unnamed branch

 

exonerate: A generic sequence comparison tool

Guy St.C. Slater. guy@ebi.ac.uk. 2000-2008.

 

 

Examples of use:

 

1. Ungapped alignment of any DNA or protein sequences:

    exonerate queries.fa targets.fa

2. Gapped alignment of Mouse proteins to Fugu proteins:

    exonerate --model affine:local mouse.fa fugu.fa

3. Find top 10 matches of each EST to a genome:

    exonerate --model est2genome --bestn 10 est.fa genome.fa

4. Find proteins with at least a 50% match to a genome:

    exonerate --model protein2genome --percent 50 p.fa g.fa

5. Perform a full Smith-Waterman-Gotoh alignment:

    exonerate --model affine:local --exhaustive yes a.fa b.fa

6. Many more combinations are possible.  To find out more:

    exonerate --help

    man exonerate

 

 

General Options:

---------------

-h --shorthelp [FALSE] <TRUE>

   --help [FALSE]

-v --version [FALSE]

 

Sequence Input Options:

----------------------

-q --query [mandatory]  <*** empty list ***>

-t --target [mandatory]  <*** empty list ***>

-Q --querytype [unknown]

-T --targettype [unknown]

   --querychunkid [0]

   --targetchunkid [0]

   --querychunktotal [0]

   --targetchunktotal [0]

-V --verbose [1]

 

Analysis Options:

----------------

-E --exhaustive [FALSE]

-B --bigseq [FALSE]

   --forcescan [none]

   --saturatethreshold [0]

   --customserver [NULL]

 

Fasta Database Options:

----------------------

   --fastasuffix [.fa]

 

Gapped Alignment Options:

------------------------

-m --model [ungapped]

-s --score [100]

   --percent [0.0]

   --showalignment [TRUE]

   --showsugar [FALSE]

   --showcigar [FALSE]

   --showvulgar [TRUE]

   --showquerygff [FALSE]

   --showtargetgff [FALSE]

   --gff3 [FALSE]

   --ryo [NULL]

-n --bestn [0]

-S --subopt [TRUE]

-g --gappedextension [TRUE]

   --refine [none]

   --refineboundary [32]

 

Viterbi algorithm options:

-------------------------

-D --dpmemory [32]

 

Code generation options:

-----------------------

-C --compiled [TRUE]

 

Heuristic Options:

-----------------

   --terminalrangeint [12]

   --terminalrangeext [12]

   --joinrangeint [12]

   --joinrangeext [12]

   --spanrangeint [12]

   --spanrangeext [12]

 

Seeded Dynamic Programming options:

----------------------------------

-x --extensionthreshold [50]

   --singlepass [TRUE]

 

BSDP algorithm options:

----------------------

   --joinfilter [0]

 

Sequence Options:

----------------

-A --annotation [none]

 

Symbol Comparison Options:

-------------------------

   --softmaskquery [FALSE]

   --softmasktarget [FALSE]

-d --dnasubmat [nucleic]

-p --proteinsubmat [blosum62]

 

Alignment Seeding Options:

-------------------------

-M --fsmmemory [64]

   --forcefsm [none]

   --wordjump [1]

 

Affine Model Options:

--------------------

-o --gapopen [-12]

-e --gapextend [-4]

   --codongapopen [-18]

   --codongapextend [-8]

 

NER Model Options:

-----------------

   --minner [10]

   --maxner [50000]

   --neropen [-20]

 

Intron Modelling Options:

------------------------

   --minintron [30]

   --maxintron [200000]

-i --intronpenalty [-30]

 

Frameshift Options:

------------------

-f --frameshift [-28]

 

Alphabet Options:

----------------

   --useaatla [TRUE]

 

Translation Options:

-------------------

   --geneticcode [1]

 

HSP creation options:

--------------------

   --hspfilter [0]

   --useworddropoff [TRUE]

   --seedrepeat [1]

   --dnawordlen [12]

   --proteinwordlen [6]

   --codonwordlen [12]

   --dnahspdropoff [30]

   --proteinhspdropoff [20]

   --codonhspdropoff [40]

   --dnahspthreshold [75]

   --proteinhspthreshold [30]

   --codonhspthreshold [50]

   --dnawordlimit [0]

   --proteinwordlimit [4]

   --codonwordlimit [4]

   --geneseed [0]

   --geneseedrepeat [3]

 

Alignment options:

-----------------

   --alignmentwidth [80]

   --forwardcoordinates [TRUE]

 

SAR Options:

-----------

   --quality [0]

 

Splice Site Prediction Options:

------------------------------

   --splice3 [primate]

   --splice5 [primate]

   --forcegtag [FALSE]

 

 

 

実行方法

ゲノムとタンパク質のfastaファイルを指定する。GFF3形式で書き出すには--gff3をつける。レポジトリでは近縁種のタンパク質を使った遺伝子モデル構築用とみられる以下の例が記載されている。

exonerate -q protein.fa -t genome.fa --model protein2genome --querytype protein --targettype dna --showvulgar no --softmaskquery yes --softmasktarget yes --minintron 20 --maxintron 3000 --showalignment no --showtargetgff yes --showcigar no --geneseed 250 --score 250 --verbose 0 --gff3 yes > out.gff3
  • --bestn 1    best hitのみ(デフォルト0
  • --percent 80  配列同一性カットオフ80% 実際のスコアが最大スコアの80%以上であるタンパク質を保持
  • --score 250   総合スコアのしきい値 250

DNA-proteinのアラインメントに利用できるモデルはいくつか存在する。ここではgenewiseに似たprotein2genomeモデルを使用している。

 

追記

alignmemtのhitファイルを作成し、EVMで使用する。その時はexonerateではGFF3形式で書き出さず、EVMのスクリプトでEVM互換GFF3形式に変換するとトラブルがない。

exonerate -q proteome.faa -t genome.fna --model protein2genome --querytype protein --targettype dna --showvulgar no --softmaskquery yes --softmasktarget yes --minintron 20 --maxintron 10000 --showalignment no --showtargetgff yes --showcigar no --geneseed 250 --score 250 --verbose 0 --gff3 no > exonerate_output
#変換
git clone https://github.com/EVidenceModeler/EVidenceModeler.git
perl EVidenceModeler/EvmUtils/misc/Exonerate_to_evm_gff3.pl exonerate_output > exonerate_output.gff3

 

引用

Automated generation of heuristics for biological sequence comparison
Guy St C Slater & Ewan Birney 
BMC Bioinformatics volume 6, Article number: 31 (2005)

 

関連


参考

Exonerate Protein2Genome Options And Output

https://www.biostars.org/p/2768/

 

こちらでは、Exonerateが遅いため、まずBLASTを実行してクエリESTやタンパク質にヒットするゲノムアセンブリの領域を探し(近似マッピング)、次にexonerateを実行してクエリEST/タンパク質をBLASTヒットの領域にギャップを考慮したアライメントするためのスクリプトが紹介されている。

http://avrilomics.blogspot.com/2013/02/using-exonerate-to-align-estscdnas-to.html