macでインフォマティクス

macでインフォマティクス

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

数百から数千のバクテリアゲノムを比較する Harvestスイート

ゲノム比較 (comparative genomics) GUIツール 結果の視覚化 (visualization) bacteria 2014 Genome Biology 分子系統樹 multiple sequence alignment (MSA)

2019 12/9 インストール更新

2020 5/31 インストール方法修正

2020 6/8 挿入文章変更、タイトル変更

2021 6/23 vcf input追記

2023/05/14オプション名の誤り修正

 

現在、多くの微生物種について全ゲノム配列が利用可能になっているが、既存の全ゲノムアラインメント手法では、複数の配列を同時に比較することができないため、その能力に限界がある。ここでは、何千もの微生物株を迅速かつ同時解析するためのコアゲノムアライメントおよび可視化ツールのHarvestスイートを紹介する。Harvestには、高速なコアゲノムマルチアライナーであるParsnpとダイナミックなビジュアルプラットフォームであるGingrが含まれている。これらを組み合わせることで、インタラクティブなコアゲノムアラインメント、バリアントコール、組換え検出、系統樹を提供する。シミュレーションと実データを用いて、既存の手法の精度を維持しつつ、我々のアプローチが他の追随を許さないスピードを発揮することを実証した。Harvestスイートはオープンソースであり、http://github.com/marbl/harvest から自由に利用できる。

 

公式ページ

http://harvest.readthedocs.io/en/latest/#

PDFマニュアル

https://media.readthedocs.org/pdf/harvest/latest/harvest.pdf

 

特徴

  • Parsnpは、近縁種のドラフトゲノムと完全長ゲノムを入力として、コアゲノムアラインメントを行い、出力としてマルチアラインメント（XMFA）、バリアント（VCF）、コアゲノム系統樹（Newick）、Gingr入力フォーマット（GGR）を返す。
  • Parsnpの主な利点は、バリアント（SNP）コールの頑健なフィルタリング、出力としての複数のアラインメント、優れたスピードにある。Parsnpは、16コアのサーバーで200-300の細菌株を30分未満でアライメントでき、数時間で~1000のアライメントが可能になっている。
  • Gingr (http://github.com/marbl/gingr) を使うと、 Parsnp の出力からマルチアラインメント、バリアント、コアゲノムアラインメントから推定される系統樹をインタラクティブに表示することができる。

 　

インストール

Github

 harvestパッケージは公式サイトからダウンロードする。

http://harvest.readthedocs.io/en/latest/

 

parsnpだけならwgetでダウンロードできる (上記のダウンロードにparsnpは入っている）。またcondaでもインストール可能。

wget https://github.com/marbl/parsnp/releases/download/v1.2/parsnp-OSX64-v1.2.tar.gz
tar -xvf parsnp-OSX64-v1.2.tar.gz

#bioconda
mamba create -n harvest -y
conda activate harvest
#parsnp (link)
mamba install -c bioconda -y parsnp
#harvest tools (format conversion)(link)
mamba install -c bioconda -y harvesttools

> parsnp

$ parsnp -h

|--Parsnp v1.2--|

For detailed documentation please see --> http://harvest.readthedocs.org/en/latest

usage: parsnp [options] [-g|-r|-q](see below) -d <genome_dir> -p <threads>

 

Parsnp quick start for three example scenarios: 

1) With reference & genbank file: 

 >parsnp -g <reference_genbank_file1,reference_genbank_file2,..> -d <genome_dir> -p <threads> 

 

2) With reference but without genbank file:

 >parsnp -r <reference_genome> -d <genome_dir> -p <threads> 

 

3) Autorecruit reference to a draft assembly:

 >parsnp -q <draft_assembly> -d <genome_db> -p <threads> 

 

[Input parameters]

<<input/output>>

 -c = <flag>: (c)urated genome directory, use all genomes in dir and ignore MUMi? (default = NO)

 -d = <path>: (d)irectory containing genomes/contigs/scaffolds

 -r = <path>: (r)eference genome (set to ! to pick random one from genome dir)

 -g = <string>: Gen(b)ank file(s) (gbk), comma separated list (default = None)

 -o = <string>: output directory? default [./P_CURRDATE_CURRTIME]

 -q = <path>: (optional) specify (assembled) query genome to use, in addition to genomes found in genome dir (default = NONE)

 

<<MUMi>>

 -U = <float>: max MUMi distance value for MUMi distribution 

 -M = <flag>: calculate MUMi and exit? overrides all other choices! (default: NO)

 -i = <float>: max MUM(i) distance (default: autocutoff based on distribution of MUMi values)

 

<<MUM search>>

 -a = <int>: min (a)NCHOR length (default = 1.1*Log(S))

 -C = <int>: maximal cluster D value? (default=100)

 -z = <path>: min LCB si(z)e? (default = 25)

 

<<LCB alignment>>

 -D = <float>: maximal diagonal difference? Either percentage (e.g. 0.2) or bp (e.g. 100bp) (default = 0.12)

 -e = <flag> greedily extend LCBs? experimental! (default = NO)

 -n = <string>: alignment program (default: libMUSCLE)

 -u = <flag>: output unaligned regions? .unaligned (default: NO)

 

<<Recombination filtration>>

 -x = <flag>: enable filtering of SNPs located in PhiPack identified regions of recombination? (default: NO)

 

<<Misc>>

 -h = <flag>: (h)elp: print this message and exit

 -p = <int>: number of threads to use? (default= 1)

 -P = <int>: max partition size? limits memory usage (default= 15000000)

 -v = <flag>: (v)erbose output? (default = NO)

 -V = <flag>: output (V)ersion and exit

> harvesttools #link

$ harvesttools 

harvesttools options:

   -i <Gingr input>

   -b <bed filter intervals>,<filter name>,"<description>"

   -B <output backbone intervals>

   -f <reference fasta>

   -F <reference fasta out>

   -g <reference genbank>

   -a <MAF alignment input>

   -m <multi-fasta alignment input>

   -M <multi-fasta alignment output (concatenated LCBs)>

   -n <Newick tree input>

   -N <Newick tree output>

   --midpoint-reroot (reroot the tree at its midpoint after loading)

   -o <Gingr output>

   -S <output for multi-fasta SNPs>

   -u 0/1 (update the branch values to reflect genome length)

   -v <VCF input>

   -V <VCF output>

   -x <xmfa alignment file>

   -X <output xmfa alignment file>

   -h (show this help)

   -q (quiet mode)

 

 

実行方法

解析の流れ

parsnpでゲノム同士のマルチプルアラインメントを行い、variants（SNVやindel）、MLによる分子系統樹を出力する。それらをGUIアプリのGingrに読み込ませて解析結果を一覧表示する。

 

parsnpによるアラインメント。

-gでリファレンスゲノムのgenbankファイルを指定する。また、-dで指定したディレクトリの中に、比較したいバクテリアのゲノム（またはアセンブルして作ったcontigを全て入れておく）。このような感じになる。

f:id:kazumaxneo:20171126204534j:plain

標準条件でのラン。

parsnp -r ref.fasta -d contig/ -p 8
  • -g　Gen(b)ank file(s) (gbk), comma separated list (default = None)
  • -d　(d)irectory containing genomes/contigs/scaffolds
  • -p　number of threads to use? (default= 1)
  • -r     (r)eference genome (set to ! to pick random one from genome dir)

 

 テストデータのラン。

 比較ゲノムデータ（harvest-master/docs/content/parsnp/mers42.tar.gzとmers49.tar.gzを解凍、またはチュートリアルよりダウンロード ）

リファレンス（mers42と49を解凍してその中にあるref/）

parsnp -r EMC_2012.fasta -d mers42/ -c
  •  -c　(c)urated genome directory, use all genomes in dir and ignore MUMi? (default = NO)

 

2021 6/23 追記

harvesttoolsを使ってユーザーが実行済みのmulti-sample VCFを選択、ggrファイル出力

harvesttools -v freebayes-joint-call.vcf -g reference.gb -o out.ggr

 

GUIアプリGingrをラウンチ。

f:id:kazumaxneo:20171126203634j:plain

open Gingr.app

 

 > FIle openでparsnp.ggrを選択。

f:id:kazumaxneo:20171126203716j:plain

 

結果が表示される。

f:id:kazumaxneo:20171126203840j:plain

 

MacBookのTrackpad（またはマウスのホイール）を２本指で上下することで、スムーズな拡大縮小が可能。

塩基が表示されるまで拡大。

f:id:kazumaxneo:20171126204134j:plain

縮小。上の方にあるcDNAのアイコンをダブルクリックすれば、cDNA全長が収まるサイズに縮尺を自動調整してくれる。

f:id:kazumaxneo:20171126204308j:plain

 左の系統樹の見たい部分をクリックすれば、その枝部分を拡大表示できる。全体に戻すにはもう一度クリックするだけ。

 

File -> export variants(VCF)を選択し、VCFファイルを出力。

f:id:kazumaxneo:20171128230047p:plain

 全ての変異が行列で出力される。そのポジションで変異があった株は1、変異がなかった株は0と表示されている。

 

他のテストデータ

https://www.cbcb.umd.edu/software/harvest

 

  • ログ出力には、コアゲノムアライメントに含まれるゲノムの割合を示すカバレッジ値が示されている。これは、Muscle aligned-regions と maximal unique matches (MUMs)を含んでいる。コアゲノムアライメントのサイズは、あるゲノムのカバレッジ値にそのゲノムの長さを乗じることで算出できる。
  • コアには必ず全ゲノムが含まれるので、参照ゲノムの選択は重要ではない。優柔不断な方や、genomeディレクトリに含まれるすべてのゲノムが同程度の品質である場合は、パラメータ '-r !' を使用してランダムにリファレンスを選択することができる。通常、完成したゲノムやクローズドゲノムが、アセンブリのアーチファクトやコンタミを含まない高品質なものであることを確認するために、参照株として使用される。
  • デフォルトでは、parsnpはリファレンスとゲノムディレクトリ内の各ゲノムとの間のMUMi距離を計算する。MUMi distance <= 0.01 のゲノムはすべて含まれ、それ以外のゲノムは破棄される。ゲノムディレクトリに存在するすべてのゲノムを含めるようにするには、コマンドラインパラメータとして '-c' を指定する。
  • parsnpの目的は、指定されたクレードのコアゲノムで見つかったすべての有益なシグナルを捕捉することにある。すべてのゲノムで共有されていない領域のSNPは報告されない。さらに、アライメントが不十分な領域で見つかったSNPも破棄される。

 

2020 6/8追記（自分の理解が間違ってました）

macでparsnpを実行すると同じフォルダにあるいくつかのfastaファイルだけを無視することがありました。同じデータをLinuxでランした時は起きなかったのでバグと考えています。注意して下さい。

引用

The Harvest suite for rapid core-genome alignment and visualization of thousands of intraspecific microbial genomes.

Treangen TJ, Ondov BD, Koren S, Phillippy AM.

Genome Biol. 2014;15(11):524.

 

2021 5/9

こちらも読んでみて下さい。