2020 7/29 論文追記
2021 10/18 パラメータ追記
細菌クローン間の遺伝子レパートリーの違いを見つけることは、病原性、抗生物質耐性、および代謝能力などの表現型の違いの遺伝的根拠を理解するために重要である。また、遺伝子の有無に関する情報を含む系統解析は、細菌集団において遺伝子が失われ、獲得されるペースとメカニズムについての情報を提供するのに役立つ[ref.1]。したがって、遺伝子の存在または非存在のゲノムワイド解析は、細菌の進化と順応をよりよく理解するために必要である。
遺伝子の有無の同定に利用可能な複数のオープンソースバイオインフォマティクスツールが存在する。各ツールには、それぞれ独自の利点と制限がある。ここでは、入力としてアセンブリされたゲノムを使用するツールに焦点を当てる。いくつかのツール、例えばPanSeq（Laing et al, 2010）は、所与の参考文献中に存在する遺伝子の存在または非存在について特異的に試験するため、クエリのゲノム配列をリファレンスゲノムへアラインメントすることに基づく方法である。このアプローチは、プロファージ遺伝子およびプラスミド遺伝子がvariable bacterial genetic contentの大部分を構成するため限られている[ref.2]。他のツール、例えばRoary [ref.3]、SaturnV [ref.4]、PanDelos [ref.5]、panX [ref.6]、BPGA [ref.7]およびEDGAR [ref.8]は、入力されたゲノムアセンブリからパンゲノムを構築し、次にそれぞれに存在する遺伝子セットを決定する。これら両方のアプローチの性能は、クエリのゲノム配列の完全性に依存し、そして全ての前述のツールは、シーケンシングおよびアセンブリ不完全性のためにわずかに遺伝子がかけているが、ほぼ完全ゲノムの分析用に設計されている。よって高度に断片化されたアセンブリ（しばしばショートリードデータに基づくde novoゲノムアセンブリの産物）、これは前述のアプローチではアセンブリ不完全性のため遺伝子がabsentであるとする多数の偽陽性の結果を招く、に対して設計された分析ツールが必要とされている。
ここでは、GenAPIという同種の細菌クローンの遺伝子量の同定と比較のためのコマンドラインツールを紹介する。遺伝子断片のみがアセンブリされている場合、または遺伝子のアセンブリが複数のコンティグに分割されている場合を説明するために、断片化されたゲノムアセンブリをうまく処理するGenAPIを開発した。シーケンシングおよびアセンブリの不完全性に対するcompensationは、trueの遺伝子欠失コール率を維持しながら、falseの遺伝子欠失コール率を最小限に抑えることを証明している。
GenAPIはUnix環境で動作するように設計されており、BLAST +、CD-HITおよびBedtoolsソフトウェアを必要とする[ref.9、10、11]。オプションの遺伝子有無マトリックスの可視化および系統樹の生成には、Rの pheatmap library [ref.12, 13]、およびRAxML [ref.14]ソフトウェアが必要である。 GenAPIの段階的なワークフローを図1Aに示す。ここで、アノテーションつきゲノム配列ファイルを入力とし、90％の同一性と80％の配列長の重複要件でCD-HIT-ESTにより遺伝子をクラスター化してパンゲノムを生成する 。 パンゲノムの遺伝子と各ゲノム配列間の最良のアラインメントが同定される。 最後に、遺伝子はユーザー提供の閾値に従って存在するまたは存在しないと定義される（デフォルト：遺伝子存在の定義は、最良のアラインメントで、98%の同一性で25%以上のカバレッジ、または90%のカバー率で50%のカバレッジがある）。 より詳細な方法の説明は補足で提供されている。
当初、このプログラムはRoaryの代替としてvery closely related な不完全細菌ゲノムのため書かれた。 このツールは、複数のリアルデータセットとシミュレーションデータセットで検証されており、真に存在しないすべての遺伝子を識別し、偽陽性コールは非常に少なくなっている。
インストール
ubuntu18.04LTSでテストした（docker使用。ホストOS; macos10.14）。macosでは正常動作しないので注意する。
依存
- BLAST >=2.6.0+
- CD-HIT >=4.6.1
- Bedtools >=2.26
Requirements for optional visualizations:
- Heatmaps:
- R >=3.2.5
- pheatmap >=1.0.10
- Phylogenetic tree (minimum 4 samples are required for the tree to be created):
- RAxML >=8.2.11
conda install -c bioconda -y blast cd-hit bedtools
#ヒートマップ画像生成やクラスタリングを行うならpheatmap、RAxMLが必要。 pheatmapはR内で
> install.packages("pheatmap", type = "source")
#RAxML (linuxならcondaで導入可)
conda install -c bioconda -y raxml
本体 Github
git clone https://github.com/MigleSur/GenAPI.git
cd GenAPI/bin/
#bin/の他のスククリプトも使うのでパスを通しておく
export PATH=$PATH:$PWD
> genapi -h
# genapi -h
Usage: GenMat [options] [--name <analysis name>]
-n, --name Chosen analysis name.
Default: date in yyyy-mm-dd format
-p, --threads Number of threads used for running the analysis.
Default: 1
-a, --coverage Alignment coverage for the shorter sequence for CD-HIT-EST.
Default: 0.8
-c, --identity Sequence identity threshold for CD-HIT-EST.
Default: 0.9
-x, --geneCov1 First minimum alignment length threshold.
Default: 0.25
-y, --geneCov2 Second minimum alignment length threshold.
Default: 0.50
-w, --geneIden1 First minimum gene identity threshold.
Default: 0.98
-z, --geneIden2 Second minimum gene identity threshold.
Default: 0.90
-l, --geneLen Minimum gene length. Shorter than the threshold genes are
excluded from the analysis.
Default: 150
-t, --tree Create a phylogenetic tree using gene presence-absence
matrix. Requires RAxML to be installed.
Default: False
-m, --matrix Create a gene presence-absence matrix visualization. Requires
Rscript and pheatmap library to be installed.
Default: False
-v, --version Print the tool version.
-h, --help Print this message.
実行方法
コマンドを実行する。カレントパスにprokkaで自動アノテーション付けして得たGFFファイル（*1）を用意しておけば自動で認識される。
#setting.shをコピーしておく
cp GenAPI/bin/settings.sh .
#実行。自動でGFFファイルを認識する
genapi
#系統樹とヒートマップも出力する
genapi -t -m -p 20
- -t Create a phylogenetic tree using gene presence-absence matrix. Requires RAxML to be installed. Default: False
- -m Create a gene presence-absence matrix visualization. Requires Rscript and pheatmap library to be installed. Default: False
- -p Number of threads used for running the analysis. Default: 1
- -a Alignment coverage for the shorter sequence for CD-HIT-EST. Default: 0.8
- -c Sequence identity threshold for CD-HIT-EST. Default: 0.9
- -x First minimum alignment length threshold. Default: 0.25
- -y Second minimum alignment length threshold. Default: 0.50
- -w First minimum gene identity threshold. Default: 0.9
- -z Second minimum gene identity threshold. Default: 0.90
- -l Minimum gene length. Shorter than the threshold genes are excluded from the analysis. Default: 150
出力
遺伝子の有無のマトリックスファイル（５サンプルからなる小さなテストゲノム使用）
アクセサリー遺伝子（variable_genes）のみのヒートマップも出力される。
より大きなデータセット
heatmap_plot_all_genes
出力の詳細はGithubで説明されています。
感想
GenAPIはランタイムが長くなる事と引き換えに、遺伝子のpresence/absenceの判定精度が高くなっています。手持ちの配列でテストすると間違いがゼロでした。Roaryなどで同じことをすると遺伝子のpresence/absenceの判断を間違うことが結構あります（特にabsenceの偽陽性が多い）。莫大なサンプル数を扱うのでないなら、類似ツールの中で、GenAPIは良い選択肢になると思います。
引用
GenAPI: a tool for gene absence-presence identification in fragmented bacterial genome sequences
Migle Gabrielaite, Rasmus L. Marvig
bioRxiv preprint first posted online Jun. 3, 2019;
2020 7/29 追記
GenAPI: a tool for gene absence-presence identification in fragmented bacterial genome sequences
Migle Gabrielaite & Rasmus L. Marvig
BMC Bioinformatics volume 21, Article number: 320 (2020)
*1
GFFにはgene featureのアノテーション情報とゲノム配列の両方の情報が含まれていなければならない。prokkaで自動アノテーション付けを行っていれば、自動で生成される。
関連
