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