macでインフォマティクス

macでインフォマティクス

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

パンゲノム解析を行う roary

2020 3/19 4/6 スクリプト修正

2020 3/19 4/10 サンプル数が多い時のオプション追記

  

 現在、典型的な原核生物のポピュレーションシーケンシング研究は、数百または数千の分離株で構成されている。 これらのデータセットを調べることで、原核生物ゲノムの遺伝構造に関する詳細な洞察を得ることができる。 ここではRoaryを紹介する。Roaryは、大規模なパンゲノムを迅速に構築し、コア遺伝子とアクセサリー遺伝子を特定するツールである。 Roaryは、結果の精度を損なうことなく、標準デスクトップ上で数千の原核生物サンプルのパンゲノムの構築を可能にする。 単一のCPU Roaryを使用すると、13 GBのRAMを使用して4.5時間で1000の分離株からなるパンゲノムを生成でき、複数のプロセッサを使用してさらに高速化できる。RoaryはPerlで実装されており、http;//sanger-pathogens.github.io/RoaryオープンソースGPLv3ライセンスの下で自由に利用できる。

 

manual

https://sanger-pathogens.github.io/Roary/

 

 

インストール

macos10.14のanaconda3.7環境でテストした。

本体 Github

conda create -n roary -y
conda activate roary

#bioconda (link)。fasttreeなども導入される(yamlファイル)。
conda install -c bioconda -y roary

roary -h 

$ roary -h

Use of uninitialized value in require at /Users/kazu/anaconda3/envs/roary/lib/site_perl/5.26.2/darwin-thread-multi-2level/Encode.pm line 61.

 

Please cite Roary if you use any of the results it produces:

    Andrew J. Page, Carla A. Cummins, Martin Hunt, Vanessa K. Wong, Sandra Reuter, Matthew T. G. Holden, Maria Fookes, Daniel Falush, Jacqueline A. Keane, Julian Parkhill,

"Roary: Rapid large-scale prokaryote pan genome analysis", Bioinformatics, 2015 Nov 15;31(22):3691-3693

    doi: http://doi.org/10.1093/bioinformatics/btv421

Pubmed: 26198102

 

Usage:   roary [options] *.gff

 

Options: -p INT    number of threads [1]

         -o STR    clusters output filename [clustered_proteins]

         -f STR    output directory [.]

         -e        create a multiFASTA alignment of core genes using PRANK

         -n        fast core gene alignment with MAFFT, use with -e

         -i        minimum percentage identity for blastp [95]

         -cd FLOAT percentage of isolates a gene must be in to be core [99]

         -qc       generate QC report with Kraken

         -k STR    path to Kraken database for QC, use with -qc

         -a        check dependancies and print versions

         -b STR    blastp executable [blastp]

         -c STR    mcl executable [mcl]

         -d STR    mcxdeblast executable [mcxdeblast]

         -g INT    maximum number of clusters [50000]

         -m STR    makeblastdb executable [makeblastdb]

         -r        create R plots, requires R and ggplot2

         -s        dont split paralogs

         -t INT    translation table [11]

         -ap       allow paralogs in core alignment

         -z        dont delete intermediate files

         -v        verbose output to STDOUT

         -w        print version and exit

         -y        add gene inference information to spreadsheet, doesnt work with -e

         -iv STR   Change the MCL inflation value [1.5]

         -h        this help message

 

Example: Quickly generate a core gene alignment using 8 threads

         roary -e --mafft -p 8 *.gff

 

For further info see: http://sanger-pathogens.github.io/Roary/

query_pan_genome -h

# query_pan_genome -h

Usage: query_pan_genome [options] *.gff

Perform set operations on the pan genome to see the gene differences between groups of isolates.

 

Options: -g STR    groups filename [clustered_proteins]

         -a STR    action (union/intersection/complement/gene_multifasta/difference) [union]

         -c FLOAT  percentage of isolates a gene must be in to be core [99]

         -o STR    output filename [pan_genome_results]

         -n STR    comma separated list of gene names for use with gene_multifasta action

         -i STR    comma separated list of filenames, comparison set one

         -t STR    comma separated list of filenames, comparison set two

         -v        verbose output to STDOUT

         -h        this help message

 

Examples: 

Union of genes found in isolates

         query_pan_genome -a union *.gff

         

Intersection of genes found in isolates (core genes)

         query_pan_genome -a intersection *.gff

         

Complement of genes found in isolates (accessory genes)

         query_pan_genome -a complement *.gff

 

Extract the sequence of each gene listed and create multi-FASTA files

         query_pan_genome -a gene_multifasta -n gryA,mecA,abc *.gff

 

Gene differences between sets of isolates

         query_pan_genome -a difference --input_set_one 1.gff,2.gff --input_set_two 3.gff,4.gff,5.gff

 

For further info see: http://sanger-pathogens.github.io/Roary/

iterative_cdhit -h

# iterative_cdhit -h

Usage: iterative_cdhit [options]

Iteratively cluster a FASTA file of proteins with CD-hit, lower the threshold each time and extracting core genes (1 per isolate) to another file, and remove them from the input proteins file.

 

Required arguments:

         -m STR   input FASTA file of protein sequences [_combined_files]

 

Options: -p INT   number of threads [1]

         -n INT   number of isolates [1]

         -c STR   cd-hit output filename [_clustered]

         -f STR   output filename for filtered sequences [_clustered_filtered.fa]

         -l FLOAT lower bound percentage identity [98.0]

         -u FLOAT upper bound percentage identity [99.0]

         -s FLOAT step size for percentage identity [0.5]

         -v       verbose output to STDOUT

         -h       this help message

 

For further info see: http://sanger-pathogens.github.io/Roary/

docker imageも準備されている。

docker pull sangerpathogens/roary

 

ここではcondaでは導入されない結果を視覚化するスクリプトも使う。Marco Galardiniさんのコードでroaryレポジトリのcontributor/に含まれている。

 依存

  • python (2 or 3)
  • Biopython
  • numpy
  • pandas
  • matplotlib
  • seaborn

#足りないものだけ導入。ここでは2つ
conda install -c bioconda -y biopython seaborn

#roaryレポジトリをとってくる
git clone https://github.com/sanger-pathogens/Roary.git

python Roary/contrib/roary_plots/roary_plots.py

$ python Roary/contrib/roary_plots/roary_plots.py 

usage: roary_plots.py [-h] [--labels] [--format {png,tiff,pdf,svg}] [-N SKIPPED_COLUMNS] [--version] tree spreadsheet

roary_plots.py: error: the following arguments are required: tree, spreadsheet

 

もう1つ、こちらはTorsten Seemann さんが書かれたコード。同じくcontributor/に含まれている。

perl Roary/contrib/roary2svg/roary2svg.pl

$ perl Roary/contrib/roary2svg/roary2svg.pl 

Usage: Roary/contrib/roary2svg/roary2svg.pl [options] gene_presence_absence.csv > pan_genome.svg

  --help          This help.

  --verbose!      Verbose output (default '0').

  --width=i       Canvas width (default '1024').

  --height=i      Row height (default '20').

  --taxacolumn=i  Column in gpa.csv where taxa begin (default '14').

  --colour=s      Colour of core cells (default 'DimGray').

  --sepcolour=s   Colour of horizontal separators/borders (default 'LightGray').

  --acconly!      Only draw accessory (non-core) genes (default '0').

  --border!       Add outline border (default '0').

 

 

実行方法

ランするのに必要なのは、prokkaで自動アノテーションして得たGFF3ファイル。NCBIからダウンロードしたGFF3は、遺伝子アノテーションの後にゲノム配列が付いていないため、使用することはできない。GenBank(full)をダウンロードし、bp_genbank2gff3.pl スクリプトlink)でGFF3に変換して使用する。その他の注意点として、結果の解釈しやすさの関係で、それぞれのサンプルの gene IDs(locus tag)はユニークがあることが望ましい。prokkaを走らせる時に"--locustag xxx"をつけて実行する。

I, annotation

prokkaを使い、カレントのfastaファイルの自動アノテーション付けを行う  。例えば*fnaを 検索してforで回す。

 1、スクリプト作成

#ヒアドキュメントで以下のスクリプトを保存。"*.fasta"を対象にアノテーション付けする。
#locus tag部分は適宜修正

#!/bin/sh 
#prooka annotation

for file in `\find *fasta -maxdepth 1 -type f`; do
 genome=${file}
folder=${file%.fasta}

prokka $genome --outdir prokka_${folder} --addgenes --cpus 12 --quiet
#rename
mv prokka_${folder}/PROKKA*.gff ${folder}.gff
done

prokka_annotation.shとして保存した。

 

 2、prokka_annotation.shを実行

#prokka導入(インストール済みならスキップ)
conda create -n prokka -c bioconda prokka -y
conda activate prokka

#prokka実行
bash prokka_annotation.sh

 

II, roary

 1、roaryの実行。

roary -e --mafft -p 8 -f output_dir *.gff

#sanple数が数百あるならgroup_limitを増やす
roary -e --mafft -p 40 --group_limit 500000 -f output_dir *.gff

#dockerだと
docker run --rm -itv $PWD:/data/ -w $PWD sangerpathogens/roary
> roary -e --mafft -p 8 -f output_dir *.gff

 2、FastTree (manual)の実行 (下の順番で指定しないとエラーになった) 。

FastTree -gtr -nt core_gene_alignment.aln > tree.newick

 3、 視覚化1

python Roary/contrib/roary_plots/roary_plots.py \
tree.newick gene_presence_absence.csv

例;192のサルモネラゲノムの視覚化結果

f:id:kazumaxneo:20200217171822p:plain

f:id:kazumaxneo:20200218194901p:plain

 

"A piechart of the pan genome, breaking down the core, soft core, shell and cloud." githubより

f:id:kazumaxneo:20200218194857p:plain

summary_statistics.txtも確認する。

 4、 視覚化2

perl Roary/contrib/roary2svg/roary2svg.pl gene_presence_absence.csv > pan_genome.svg

12ゲノムの結果

f:id:kazumaxneo:20200218194335p:plain

 

ⅡI, query_pan_genome 

コア遺伝子(その種で共通して存在)とアクセサリ遺伝子(株間で共通して存在)、サンプルに特異的な遺伝子を、共通、交差、ユニークなどに分けてリスト出力する。使用するには前もって1のroaryを実行し、clustered_proteinsを用意しておく必要がある。また、gffファイルも必要なので、roaryを指定ディレクトリに出力しているなら、*gffかclustered_proteinsのいずれかについて、カレントへのシンボリックリンクを張っておく必要がある。

#全ゲノムで見つかる遺伝子(core gene)
query_pan_genome -a intersection *.gff -o intersect_outout

#2ゲノム以上(≥2)で共通の遺伝子
query_pan_genome -a union *.gff -o union_outout


#全ゲノムには共通していない遺伝子 (accessory genes)(1回しか見つからないという意味ではない)
query_pan_genome -a complement *.gff -o complement_outout

roaryと違って該当ファイルのみ出力される。

comp_outputの例

f:id:kazumaxneo:20200214142254p:plain

 

query_pan_genome -a complement を走らせて5つのゲノムアノテーションを比較したところ、3ゲノムにしか見つからない遺伝子lldRがあったとする。

f:id:kazumaxneo:20200219104429p:plain

 

この遺伝子配列を取り出す。

使用した全てのgffを対象に”query_pan_genome -a gene_multifasta ”を走らせれば簡単に取り出せる。

該当する遺伝子のタンパク質配列を抽出

query_pan_genome -a gene_multifasta -o test -n lldR *.gff

#複数同時も可能
query_pan_genome -a gene_multifasta-o out -n gryA,mecA,abc *.gff

 出力されるtest_lldR.faを開く。

> cat test_lldR.fa

>genome1_00013

MKSATSAQRPYQEVGAMIRDLIIKTPYNPGERLPPEREIAEMLDVTRTVVREALIMLEIK

GLVEVRRGAGIYVLDNSGSQNTDSPDANVCNDAGPFELLQARQLLESNIAEFAALQATRE

DIVKMRQALQLEERELASSAPGSSESGDMQFHLAIAEATHNSMLVELFRQSWQWRENNPM

WIQLHSHLDDSLYRKEWLGDHKQILAALIKKDARAAKLAMWQHLENVKQRLLEFSNVDDI

YFDGYLFDSWPLDKVDA*

>genome2_00013

MKSATSAQRPYQEVGAMIRDLIIKTPYNPGERLPPEREIAEMLDVTRTVVREALIMLEIK

GLVEVRRGAGIYVLDNSGSQNTDSPDANVCNDAGPFELLQARQLLESNIAEFAALQATRE

DIVKMRQALQLEERELASSAPGSSESGDMQFHLAIAEATHNSMLVELFRQSWQWRENNPM

WIQLHSHLDDSLYRKEWLGDHKQILAALIKKDARAAKLAMWQHLENVKQRLLEFSNVDDI

YFDGYLFDSWPLDKVDA*

>genome3_00013

MKSATSAQRPYQEVGAMIRDLIIKTPYNPGERLPPEREIAEMLDVTRTVVREALIMLEIK

GLVEVRRGAGIYVLDNSGSQNTDSPDANVCNDAGPFELLQARQLLESNIAEFAALQATRE

DIVKMRQALQLEERELASSAPGSSESGDMQFHLAIAEATHNSMLVELFRQSWQWRENNPM

WIQLHSHLDDSLYRKEWLGDHKQILAALIKKDARAAKLAMWQHLENVKQRLLEFSNVDDI

YFDGYLFDSWPLDKVDA*

 

 

複数ゲノム間で異なる遺伝子を探す。

query_pan_genome -a difference --input_set_one 1.gff,2.gff --input_set_two 3.gff,4.gff,5.gff

 

こちらで紹介されているツールを使うと、sqliteのデータベースを作り、パンゲノム中の指定した遺伝子配列を抽出できる。例えばコア遺伝子rpoAに着目して、全ゲノムからのrpoA遺伝子配列を抽出する。

 

roaryの出力はpanX、Gubbins、scoary、FriPan、PanVizGeneratorなど様々なツールの入力に指定できる。詳細はroaryのHP参照(link)。 

引用

Roary: rapid large-scale prokaryote pan genome analysis
Andrew J. Page, Carla A. Cummins, Martin Hunt, Vanessa K. Wong, Sandra Reuter, Matthew T.G. Holden, Maria Fookes, Daniel Falush, Jacqueline A. Keane, Julian Parkhill 
Bioinformatics, Volume 31, Issue 22, 15 November 2015, Pages 3691–3693

 

関連


 

参考