macでインフォマティクス

macでインフォマティクス

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

高速でスケーラブルなパンゲノミクスツールボックス PIRATE

2019 GigaScience docker pan-genome tRNA rRNA 高速なツール

2022/10/17 追記

 

　天然の細菌集団内の遺伝子分布をカタログ化することは、進化の過程や適応の遺伝的基盤を理解するために不可欠である。全ゲノム配列決定技術の進歩により、公開データベースに登録される細菌ゲノムの量は膨大なものとなっている。このような大規模なゲノムデータセットから、遺伝子やその他の特徴を抽出するためのソフトウェアが求められている。本発表では、パンゲノミクスツールボックスPIRATE (Pangenome Iterative Refinement and Threshold Evaluation)を紹介する。これは、バクテリアのパンゲノムから、幅広い配列類似度の閾値を越えてオルソログ遺伝子ファミリーを同定・分類するツールボックスである。PIRATEは、最近のスケーラブルなソフトウェア開発により、数千の分離株を迅速に調査することができる。PIRATEは、アミノ酸やヌクレオチドの同一性の閾値を幅広く設定し、遺伝子（またはその他の注釈付き特徴）をクラスタリングし、その情報を使ってパラロガス遺伝子ファミリーや推定分裂・融合イベントを迅速に同定する。さらに、PIRATEは有向グラフを用いてパンゲノムを順序付けし、対立遺伝子変異の指標を提供し、各遺伝子ファミリーの配列発散を推定することができる。

　PIRATEはサンプル数と計算資源に比例してスケールアップし、大規模なゲノムデータの解析が可能であり、他の一般的なツールとも比較できることを実証した。PIRATEはバクテリアのパンゲノムを解析するための堅牢なフレームワークを提供し、クローン性の高い種から汎発性の高い種まで幅広く対応できる。

 

Wiki

https://github.com/SionBayliss/PIRATE/wiki/Usage

 

インストール

python3.7の環境を作ってテストした。

Github

#conda (link)
mamba create -n pirate python=3.7 -y
conda activate pirate
mamba install -c conda-forge -c bioconda pirate -y

#導入に苦労したのでdockerイメージを上げておきます
docker pull kazumax/pirate
docker run -itv $PWD:/data -w /data --rm kazumax/pirate
> conda activate pirate

> PIRATE -c #check

# PIRATE -c

 

Running PIRATE on test files:

 

 - PIRATE completed with no errors

 

 - tests completed

 

 - temporary files saved to /tmp/5YHN04MH6E

> PIRATE -h

# PIRATE -h

Usage:

            PIRATE -i /path/to/directory/containing/gffs/

 

     PIRATE input/output:

     -i|--input     input directory containing gffs [mandatory]

     -o|--output    output directory in which to create PIRATE folder

                    [default: input_dir/PIRATE]

 

     Global:

     -s|--steps     % identity thresholds to use for pangenome construction

                    [default: 50,60,70,80,90,95,98]

     -f|--features  choose features to use for pangenome construction.

                    Multiple may be entered, seperated by a comma [default: CDS]

     -n|--nucl      CDS are not translated to AA sequence [default: off]

     --pan-opt      additional arguments to pass to pangenome_contruction

     --pan-off      don't run pangenome tool [assumes PIRATE has been previously

                    run and resulting files are present in output folder]

 

     Paralog classification:

     --para-off     switch off paralog identification [default: on]

     --para-args    options to pass to paralog splitting algorithm

                    [default: none]

     --classify-off do not classify paralogs, assumes this has been

                    run previously [default: on]

 

     Output:

     -a|--align     align all genes and produce core/pangenome alignments

                    [default: off]

     -r|--rplots    plot summaries using R [requires dependencies]

 

     Usage:

     -t|--threads   number of threads/cores used by PIRATE [default: 2]

     -q|--quiet     switch off verbose

     -z             retain intermediate files [0 = none, 1 = retain pangenome

                    files (default - re-run using --pan-off), 2 = all]

     -c|--check     check installation and run on example files

     -h|--help      usage information

 

 

実行方法

PIRATE をアミノ酸 %ID の閾値の範囲 ("-s 50,70,90,95") で実行し、パラログを分類し、出力テーブルをディレクトリに生成する。単一の閾値なら-pを使う("-p 98")。個々の遺伝子配列をMAFFTでアラインメントし、コアゲノムアラインメントまで行うなら"-a"を付ける。ランするにはGFF3ファイル（full）を指定する。ファイルの末尾にヌクレオチド配列を含むGFF3アノテーションファイルを指定する必要がある（Prokkaが生成するフォーマットのため、prokkaでアノテーションを付けることが推奨されている）。パラログ分類をオフにするには"--para-off"を付ける。より高速なdiamodを使うには-k "--diamond"を付ける（タンパク質のみ）。

PIRATE -i gff_files/ -s "50,70,90,95" -o output_dir -a -r -t 20
  • -i    input directory containing gffs [mandatory]
  • -o   output directory in which to create PIRATE folder [default: input_dir/PIRATE]
  • -s    % identity thresholds to use for pangenome construction [default: 50,60,70,80,90,95,98]
  • -p    single % identity threshold to use for pangenome construction [default: 98]
  • -a    align all genes and produce core/pangenome alignments [default: off]
  • -r     plot summaries using R [requires dependencies]
  • --para-off   switch off paralog identification [default: on]
  • -t    number of threads/cores used by PIRATE [default: 2]

 

出力例

 

 

PIRATE.pangenome_summary.txt - パンゲノムに含まれる遺伝子の数と頻度の簡単なサマリー。


PIRATE.gene_families.ordered.tsv  - すべての遺伝子ファミリーを表形式でまとめたもの。各遺伝子ファミリーにつき1エントリーで、パラログ分割の段階で分離されたファミリーはアンダースコアと数字で示される（例：g0001_1、g0001_2）

 

PIRATE.unique_alleles.tsv - 各遺伝子ファミリーのユニークアレルを表形式でまとめたファイル。ユニークアレルとは、%identityの閾値が高い遺伝子座の新規MCLサブクラスターと定義される。

PIRATE.gene_families.ordered.tsv とPIRATE.unique_alleles.tsv の各列の詳細はレポジトリで説明されています。

 

binary_presence_absence.fasta/nwk

遺伝子ファミリーの有無のバイナリ表現データを使いfasttreeで生成したツリー。そして、その生成に使用したfastaファイル。

 

pangenome.gfa - 遺伝子ファミリー間の全てのユニークな接続を表現したGFAネットワークファイル（GFFファイルから抽出）。Bandageで読み込み、可視化することができる。

 

modified_gffs/ 

PIRATE用に標準化されたGFF3ファイル

 

feature_sequences/

各遺伝子ファミリーの全アミノ酸・ヌクレオチド配列（MAFFTでアラインメント）が格納されている。

 

representative_sequences.faa、representative_sequences.ffn - 各遺伝子ファミリーの代表的なマルチファスタ配列を塩基配列（ffn）およびアミノ酸配列（faa）で表したもの。遺伝子ファミリーごとの最長配列を代表配列として選択（ゲノムはアルファベット順）。分離株と遺伝子ファミリーの情報がFastaヘッダに含まれる。

 

core_alignment/pangenome_alignment.fasta - MAFFTを用いて作成したコアおよびフルパンゲノムの遺伝子ごとのヌクレオチドアラインメント。遺伝子はPIRATE.gene_families.ordered.tsvファイルにより順序付けされている。また、パンゲノムが翻訳CDSから作成された場合は、遺伝子のコドン構造を保持するために、アミノ酸配列から逆翻訳してアラインメントが作成されている。

 

注 - 遺伝子ファミリーの遺伝子量／コピー数が1以上のゲノムでは、アラインメント中のseqeuenceを?sに置き換えています。

 

PIRATE_plots.pdf (”-r”を付けてランした時)

 

  • （論文より）遺伝子は、正選択と純化選択の影響により、異なる速度で変異を獲得する。そのため、特定の同一性の閾値を超えたら同じファミリーに属さなくなるという定義は困難である。PIRATEは、バクテリアのパンゲノムから、幅広い配列類似度閾値でオルソログ遺伝子ファミリーを同定・分類する（この機能により、単一の閾値では損失してしまうコア遺伝子を逃しいくくなる。精度を上げるのは細かく数値を指定する。上限は記憶が正しければ98％。下限はどこまで下げてもよいのかは難しいところ。同種なら95%など高い値に止めたりするが、別種同属の距離まで離れると難しい）。
  • パラロガス遺伝子ファミリーを推定遺伝子分裂/融合イベントまたは遺伝子重複のいずれかに分類する。
  • アミノ酸同一性ではなく、ヌクレオチド同一性を用いてCDSの特徴でパンゲノムを作成するには”-n”を付ける。
  • tRNA および rRNAを使って PIRATE を実行するには"-f"を付ける。例えば-f "rRNA,tRNA"。
  • 遺伝子アライメントを再作成し、ゲノム、アレル、目的の遺伝子に対してフィルタリングを行うことができるsubset_alignments.plが用意されている。
  • 遺伝子ファミリー／アレルに対する代表配列を特定するselect_representativeコマンドが用意されている。
  • 各アレル/遺伝子ファミリーの存在/非存在のバイナリテーブルに変換するPIRATE_to_Rtab.plが用意されている。入力はPIRATE.*.tsvファイル。
  • PIRATEファイル形式をroaryのgene_presence_absence.csvファイル形式に変換するPIRATE_to_roary.plが用意されている。入力はPIRATE.*.tsvファイル。これにより、PIRATEの出力を、roaryの下流解析ツール(phandango/scoaryなど)に対応させることができる。また、PanXユーザーインターフェース、SCOARY、Microreact、Phandangoなど、下流の解析に使用するソフトウェアへの入力として使用することができる。

気づいた事

  • 各遺伝子について調べたい時は出力ディレクトリにあるサブディレクトリ；feature_sequences/にアクセスする。グループ化された遺伝子ごとに配列がmulti-fasta形式で塩基配列とアミノ酸配列でまとめられており、allele_nameでアクセスできる。例えばg014314_000019.aa.fastaには18個のタンパク質配列が含まれる。
  • 属レベルのような広範な分類でコア遺伝子を探索する場合、閾値を下げて探索しないと配列の違いが激しいオルソログが見つけられなくなる。そのため閾値を十分に下げて探索せざるを得なくなるが、その場合、homologousなタンパク質（オルソログとパラログ）やanalogousなタンパク質（収斂進化によって独自に進化した同じ機能を持つタンパク質）と、偶然配列が似ているためにヒットした機能が異なるタンパク質間、の間の区別はますます難しくなる。その場合、まず出力テーブルのmin_length(bp)やmax_length(bp)、average_length(bp)列に注目する。通常、同じ機能を持つ進化的起源が同一のタンパク質の長さは大きく変わらない。minやmaxの長さの数値が平均の列の数値と大きく違う場合、異なる機能のタンパク質が混ざり込んでいる可能性がある。またfeature_sequences/にあるMSAファイルを取り出し、混在が疑われるタンパク質を加えて系統推定することも考えられる。

 

追記

PIRATEはlocus IDを標準的な形式（ゲノム名[アンダースコア]locus番号）に変更して扱う。元の名前に戻すにはsubsample_outputs.pl スクリプトを使う。

subsample_outputs.pl -i PIRATE.gene_families.ordered.tsv -g modified_gffs/ -o
PIRATE.gene_families_ordered_rename.tsv --field "prev_locus"

 

 

引用
PIRATE: A fast and scalable pangenomics toolbox for clustering diverged orthologues in bacteria
Sion C Bayliss, Harry A Thorpe, Nicola M Coyle, Samuel K Sheppard, Edward J Feil

Gigascience. 2019 Oct 1;8(10):giz119

 

関連