macでインフォマティクス

macでインフォマティクス

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

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

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

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

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

 

現在、多くの微生物種について全ゲノム配列が利用可能になっているが、既存の全ゲノムアラインメント手法では、複数の配列を同時に比較することができないため、その能力に限界がある。ここでは、何千もの微生物株を迅速かつ同時解析するためのコアゲノムアライメントおよび可視化ツールの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

  

インストール

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
conda create -n harvest -y
conda activate harvest
#parsnp (link)
conda install -c bioconda -y parsnp
#harvest tools (format conversion)(link)
conda 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 -g ref.gbk -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)

 

 

 テストデータのラン。

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

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

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

 

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(またはマウスのホイール)を2本指で上下することで、スムーズな拡大縮小が可能。

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

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

 

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.