macでインフォマティクス

macでインフォマティクス

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

細菌ゲノムとプラスミドの系統に基づく比較ゲノムパイプライン GEnView

 

 ある細菌遺伝子のゲノム座を株や種を超えて比較することで、後天的な移動性、異なる分類群間での保存の度合い、あるいは遺伝子の水平伝播事象の示唆など、その進化に関する洞察を得ることができる。現在までに数千の細菌ゲノムが利用可能であるが、多数のゲノムの個々の遺伝子座の比較を容易にするソフトウェアは存在しない。GEnView は、多数の細菌ゲノムの遺伝子座を比較解析するための Python ベースのパイプラインで、現在 NCBI Assembly および RefSeq データベースに登録されている 840.000 以上のゲノムおよびプラスミドに分類群別に自動アクセスできるほか、 NCBI に登録されていないローカルゲノムを処理し、GEnView が作成する対話的可視化と豊富なメタデータファイルによってゲノム配列検索とその遺伝子環境の解析が可能にする。

GEnViewはPython 3で実装されている。ダウンロードと使用方法は、https://github.com/EbmeyerSt/GEnView の GLP3 に掲載されている。

 

tutorial

Home · EbmeyerSt/GEnView Wiki · GitHub

 

インストール

依存

We recommend to set up a virtual environment and install all dependencies to run GEnView there.

  • Python >=3.6
  • pandas >=1.0.3
  • sqlite3 >=2.6
  • Biopython >=1.68

The following softwares should be located in your $PATH:

  • Diamond >=0.9.24
  • Prodigal >=2.6.3
  • CD-hit >=4.7
  • mafft >=7.310
  • FastTree >=2.1.9

Github

mamba install -c bioconda genview -y

>  genview-makedb -h

usage: genview-makedb [-h] -d TARGET_DIRECTORY -db DATABASE [-p PROCESSES] [-id IDENTITY] [-scov SUBJECT_COVERAGE] [--is_db] [--uniprot_db UNIPROT_DB] [--uniprot_cutoff UNIPROT_CUTOFF] [--taxa TAXA [TAXA ...]] [--assemblies] [--plasmids]

                      [--custom] [--save_tmps] [--acc_list ACC_LIST] [--flanking_length FLANKING_LENGTH]

 

Creates sqlite3 database with genetic environment from genomes containing the provided reference gene(s).

 

optional arguments:

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

  -d TARGET_DIRECTORY, --target_directory TARGET_DIRECTORY

                        path to output directory

  -db DATABASE, --database DATABASE

                        fasta/multifasta file containing amino acid sequences of translated genes to be annotated

  -p PROCESSES, --processes PROCESSES

                        number of cores to run the script on

  -id IDENTITY, --identity IDENTITY

                        identity cutoff for hits to be saved to the database (e.g 80 for 80% cutoff)

  -scov SUBJECT_COVERAGE, --subject_coverage SUBJECT_COVERAGE

                        minimum coverage for a hit to be saved to db (e.g 80 for 80% cutoff)

  --is_db               database containing IS, integrons, ISCR sequences

  --uniprot_db UNIPROT_DB

                        path to database for annotation of surrounding sequences in ".dmnd" format. If unspecified default uniprotKB database will be downloaded to target directory

  --uniprot_cutoff UNIPROT_CUTOFF

                        % identity threshold for annotating orfs aurrounding the target sequence, default 60

  --taxa TAXA [TAXA ...]

                        taxon/taxa names to download genomes for - use "all" do download all available genomes, cannot be specified at the same time as --acc_list

  --assemblies          Search NCBI Assembly database 

  --plasmids            Search NCBI Refseq plasmid database

  --custom              Search custom genomes

  --save_tmps           keep temporary files

  --acc_list ACC_LIST   csv file containing one genome accession number per row, cannot be specied at the same time as --taxa

  --flanking_length FLANKING_LENGTH

                        Max length of flanking regions to annotate

> genview-visualize -h

usage: genview-visualize [-h] -gene GENE -db DB -id ID [-nodes NODES] [-taxa TAXA [TAXA ...]] [--force] [--compressed]

 

Extract, visualize and annotate genes and genetic environments from genview database

 

optional arguments:

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

  -gene GENE            name of gene/orf to extract and visualize

  -db DB                genview database created by genview-create-db

  -id ID                percent identity threshold for genes to extract

  -nodes NODES          should nodes be connected to genome with solid line (solid), connected by dashed line (dash) or no connection (none)

  -taxa TAXA [TAXA ...]

                        list of genera and/or species to extract

                        By default all taxa are extracted

  --force               Force new alignment and phylogeny

  --compressed          Compress number of displayed sequences, helpful with large number of identical sequences

 

 

実行方法

ここではGEnViewのチュートリアル通り進めます。

 

1、NCBIから分類を指定してゲノムをダウンロードする。

> cat PER.fa

>PER-1
MNVIIKAVVTASTLLMVSFSSFETSAQSPLLKEQIESIVIGKKATVGVAVWGPDDLEPLL
INPFEKFPMQSVFKLHLAMLVLHQVDQGKLDLNQTVIVNRAKVLQNTWAPIMKAYQGDEF
SVPVQQLLQYSVSHSDNVACDLLFELVGGPAALHDYIQSMGIKETAVVANEAQMHADDQV
QYQNWTSMKGAAEILKKFEQKTQLSETSQALLWKWMVETTTGPERLKGLLPAGTVVAHKT
GTSGIKAGKTAATNDLGIILLPDGRPLLVAVFVKDSAESSRTNEAIIAQVAQTAYQFELK
KLSALSPN

各タンパク質配列のヘッダー名はできるだけシンプルにし、'|'を含まないようにする。複数配列使う時は、似たような名前のヘッダー名にする。例えば、PER型β-lactamaseの異なる変異体をPER-1, PER-2, PER-6と命名し、後でデータベースからPERという共通名で抽出する。

genview-makedbコマンドを使ってNCBIから指定した種のゲノム(指定すればプラスミドも) を自動的にダウンロードし、指定した参照タンパク質について検索を行う。ここではPER-1ファミリーβ-ラクタマーゼをコードするタンパク質配列を参照タンパク質とする(アミノ酸配列がチュートリアルに載ってます)。参照タンパク質と出力ディレクトリ、分類のrheinheimera(Genus)を指定する。アミノ酸同一性の閾値は80%、10 CPU指定。

genview-makedb -d out_dir -db PER.fna -p 10 -id 80 -scov 80 --taxa 'rheinheimera' --assemblies
  • -d     path to output directory
  • -db   fasta/multifasta file containing amino acid sequences of translated genes to be annotated
  • -p     number of cores to run the script on
  • -id    identity cutoff for hits to be saved to the database (e.g 80 for 80% cutoff)
  • --taxa    taxon/taxa names to download genomes for - use "all" do download all available genomes, cannot be specified at the same time as --acc_list

このコマンドでは、ダウンロードしたゲノムは、DIAMOND blastxを用いて、参照遺伝子とユーザー定義のIDおよびカバレッジ閾値で検索される。ゲノムから参照遺伝子が見つかると、その遺伝子と上下流20kbまでが抽出され、Prodigalを使用してORFが特定される。これらは、Uniprotナリッジベースのカスタムバージョンとに対してユーザー定義のID閾値でDIAMONDを使用して検索される。結果は、SQLite3 を使用してデータベースに保存される。

 

出力例

f:id:kazumaxneo:20220315004804p:plain

結果はSQLite データベースファイルgenview_database.dbに保存される。このファイルには、ターゲット遺伝子がどのゲノムに存在され、予測され、アノテーションされた遺伝子なのか、全ての情報が含まれている。ジョブが終わるまで5分くらいかかった。

--taxa 'all'を指定すると全データ(4TB)がダウンロードされ、処理だけで数日かかるので、推奨されていない。

 

補足

複数の分類群や、特定のゲノムアセンブリだけダウンロードしたい時は、1行に1つのアクセッションを書いたアクセッションリストを--acc_list オプションで提供する。--acc_list は--taxaオプションとは排他的関係にあり、同時指定はできない。

GCA_000217935.2
GCA_003990335.1
GCA_001275035.1
GCA_000986865.1
GCF_003862465.1
GCA_004005375.1

genview-makedb -d out_dir -db PER.fna -p 10 -id 80 -scov 80 --assemblies  --acc_list list.txt

 

 

2、配列の系統に基づいた視覚化を行う。遺伝子名(ヘッダのPER-1)とアミノ酸同一性の閾値を指定する。

genview-visualize -db outdir/genview_database.db -gene PER-1 -id 80
  • -gene     name of gene/orf to extract and visualize
  • -db      genview database created by genview-create-db
  • -id       percent identity threshold for genes to extract
  • -taxa   list of genera and/or species to extract. By default all taxa are extracted

第2ステップでは、ユーザーが指定した遺伝子名についてデータベースから相補的なファイルが抽出される。それぞれの遺伝子とその遺伝子環境は、FASTA形式のファイルとして保存される。これらの配列は MAFFTでアライメントされ、このアライメントに基づく系統樹がFastTree v2.1.9で計算される。

 

SQLite データベースがあるパスにディレクトリができる。

f:id:kazumaxneo:20220315010850p:plain

interactive_visualization.html

f:id:kazumaxneo:20220315010835p:plain

Optionから視覚化パラメータを調整できる。

f:id:kazumaxneo:20220315010944p:plain

optionからORFを拡大したり、画面に見えている範囲の画像を出力できる。

f:id:kazumaxneo:20220315115428p:plain

 

出力ファイル

  • annotation_meta.csv - 対象遺伝子の抽出範囲に含まれるアノテーション遺伝子の情報(名前、位置、配列など)

    visualization_meta.csv - 対象遺伝子の抽出範囲に含まれる注釈付き遺伝子の情報(名前、位置、配列など)を、対話型HTML可視化出力ファイルの作成に使用する形式で格納

    yourgenename_contexts.fna - ゲノムごとに抽出された配列(標的遺伝子+遺伝的環境)を含むFASTAファイル

    yourgenename_contexts.unique.fna - 重複を除去した各ゲノム毎の抽出配列(標的遺伝子+遺伝的環境)を含むFASTAファイル

    yourgenename_contexts.unique.aln 抽出したユニークな配列のアライメント

    yourgenename_contexts.unique.tree FastTree2により作成されたツリーファイル

    interactive_output.html ツリーと抽出された配列(標的遺伝子+遺伝的環境)を対話的に可視化したもの(上の写真)。

 

引用

GEnView: a gene-centric, phylogeny-based comparative genomics pipeline for bacterial genomes and plasmids
Stefan Ebmeyer, Roelof Dirk Coertze, Fanny Berglund, Erik Kristiansson, D G Joakim Larsson

Bioinformatics. 2021 Dec 24;38(6):1727-1728