macでインフォマティクス

macでインフォマティクス

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

高カバレッジな細菌ゲノムのdenovoゲノムアセンブリツール HGA

assembly bacteria k-mer

 

　デノボゲノムアセンブリにはgreedy strategy、string overlap graph、そしてde Bruijn graphの3つの主なアプローチがある。greedy strategyは、シードリードを選択し、最大のオーバーラップが可能になるまで貪欲に拡張していくことによって機能する。このアプローチは、SSAKE [論文より ref.1]、SHARCGS [ref.2]、VCAKE [ref.3]などの初期のアセンブラで採用された。残念ながらこのアプローチは、リピートおよびシーケンシングエラーによって引き起こされるあいまいさを考慮せず、多数の誤ったアセンブリエラーを招く。string overlap graphのアプローチは、ノードとエッジとしてのリードと、すべてのノードの対を接続するエッジを有するグラフを、重複が最小の状態で構築することに基づいている。オーバーラップグラフを作成するには、計算集約的な全ペアごとの比較ステップが必要となる。グラフを構築した後、リードのレイアウトを計算し、コンセンサス配列をmultiple アライメントを用いて決定する。このアプローチは、Newbler [ref.4]やSGA [ref.5]のようなアセンブラで実装され、CABOG [ref.6]はSangerと454のシーケンシングのような長いリードに対してより効率的な方法になっている。

　de Bruijn graphモデル[ref.7]に基づく第3のアプローチは、ABySS [ref.8]、ALLPATHS-LG [ref.9]、Euler-USR [ref.10]、MaSuRCA [ref.10] (OLC+de Bruijn graph) 、SoapDenovo2 [ref.12]、SPAdes [ref.13]およびVelvet [14]など、NGSデータを対象とするアセンブラで最も一般的に使用されている[ref.11]。 de Bruijn graphの作成は、リードを長さk（k-mersと呼ばれる）の全部分文字列を集めることから始まる。それから k-mer配列aのk-1の長さがk-mer配列bの長さk-1の長さと一致しk+1の配列を得られるならば、aとbを重複させることによって拡張する。 de Bruijnグラフは線形時間で構築できるが、格納するには非常に大量のメモリが必要となり、通常はoverlap graphよりもはるかに大きくなる。 de Bruijn graphを作成した後、各アセンブラでは、ゲノムの反復やヘテロガスな部位によって引き起こされるcyclesやbulgesなどのグラフ構造をヒューリステックなアプローチで単純化する。最後にアセンブラは、最終的なコンティグを形成するde Bruijn graphの単純なパスを選択する。ゲノムアセンブリのアルゴリズムの詳細については、[ref.15-17]に記載されている。

　GAGE-Bの研究 [ref.20]では、バクテリアゲノムについて、高いカバレッジ（100-300×）のショートインサートライブラリをシーケンスして、MaSuRCAおよびSPAdes を単独で実行することで高いアセンブリ品質が得られることを示している。本論文のHGAは、このような高いカバレッジを利用するために、独立したサブセットのリードをアセンブルし、できたサブセットのアセンブリを結合し、最後にオリジナルのリードを使い結合されたコンティグを再アセンブリする。より高いカバレージと短いkmerよりも、低いカバレッジと短いkmerを使用する方がtips、bubbles、bulges、cycles、false branchは少なくなる。結果として得られるコンティグは、コンティグがより短くなるが、エラーがより少なくなる。この方法を実行するため、HGAは最初に読み込み全体を複数のパーティションに分割し、各パーティション内のカバレッジを低くしている。ただし低いカバレッジを使用することによって得られるコンティグは、カバレッジが高いよりも短くなり、より多くのギャップを有する。これを解決するために、各パーティションのアセンブリから生成された全てのコンティグをマージまたはアセンブルする。全てのパーティションから生成されたコンティグを組み立てることで、共通のコンティグ（より正確なコンティグである可能性が高い）を選択し、任意の重複または偽（部分または完全）コンティグをフィルタリングしている。100bpのIllumina HiSeqと250bpのIllumina MiSeqのデータからなる7つのGAGE-Bバクテリアデータセットを用いて、この方法論を経験的に評価している。結果は、HGAがN50および修正N50メトリックに基づいてアセンブリの品質を大幅に向上させることを示している。

  

インストール

環境

依存

  • Velvet（make 'MAXKMERLENGTH=111' 'LONGSEQUENCES=1'をつけてビルドしていること、不明ならランのシェルスクリプトを真似する）

本体 Github

https://github.com/aalokaily/Hierarchical-Genome-Assembly-HGA

git clone https://github.com/aalokaily/Hierarchical-Genome-Assembly-HGA.git
cd Hierarchical-Genome-Assembly-HGA/
python HGA.py

 

$ python HGA.py

 

Hierarchical Genome Assembly version 1.0.0

--Prerequisite-- 

Velvet should be installed with option 'LONGSEQUENCES=1' in the make command, to allow Velvet to accept contigs (long sequences) as input; as well and 'MAXKMERLENGTH=111' as the default max kmer length is 31

 

--Parameters--: 

 -velvet  path      path to the directory where velvet bineries (velveth, velvetg) are located 

 -spades  path      path to the directory where spades.py file is located

 -PA      string    select "SPAdes" or "velvet" to choose which assembler will be used to assemble the partitions

 -P12     file      file of the interlaced fastq file of paired reads that will be used in the partiton step, usually raw or clean (corrected) reads

 -R12     file      file of the interlaced fastq file of paired reads that will be used in the re-assembly step, usually raw or clean (corrected) reads

 -ins     int       insert size of the fragments

 -std     int       standard deviation of the inset size

 -P       int       number of partitions

 -Pkmer   int(odd)  kmer value for assembling the partitions

 -Rkmer   int(odd)  kmer value for re-assembly step 

 -t       int       number of threads to be used for re-assembly step using SPAdes assembler, default 1.

 -out     Path      output path

 -h                 print option menue

 

 

ラン

テストランを行うスクリプトSample_test_cholera.shとSample_test_axonopodis.shが準備されている。これを動かしてみる。

chmod u+x Sample_test_cholera.sh
./Sample_test_cholera.sh

GAGE-B（リンク）からV_cholerae_HiSeq.tar.gzがダウンロードされ、解凍された後interleaveのfastqに分離される。spadesとvelvetがインストールされ、以下のコマンドが実行される。

curr=`pwd`

python HGA.py -velvet $curr/velvet_1.2.10 -spades $curr/SPAdes-3.0.0/bin -PA velvet -P12 $curr/data.fastq -R12 $curr/data.fastq -ins 335 -std 35 -P 4 -Pkmer 31 -Rkmer 51 -t 20 -out ./HGA
  • -P12   file of the interlaced fastq file of paired reads that will be used in the partiton step, usually raw or clean (corrected) reads
  • -R12    file of the interlaced fastq file of paired reads that will be used in the re-assembly step, usually raw or clean (corrected) reads 
  • -ins    insert size of the fragments
  • -std    standard deviation of the inset size
  • -P    number of partitions
  • -Pkmer    (odd) kmer value for assembling the partitions
  • -Rkmer    (odd) kmer value for re-assembly step
  • -t    number of threads to be used for re-assembly step using SPAdes assembler, default 1.
  • -out    output path

 終わると出力ディレクトリにいくつかのファイルができる。part〜はサブセットのアセンブルディレクトリになる。上では-P 4なので、４つに分割してアセンブルしている。

$ ls -alth HGA/

total 848M

-rw-r--r--  1 uesaka user 174K Mar 15 12:02 HGA.log

drwxr-xr-x  7 uesaka user 4.0K Mar 15 12:02 HGA_combined

drwxr-xr-x  7 uesaka user 4.0K Mar 15 11:33 HGA_merged

drwxr-xr-x  2 uesaka user 4.0K Mar 15 11:06 combined_contigs

drwxr-xr-x  2 uesaka user 4.0K Mar 15 11:05 merged_contigs

drwxr-xr-x  2 uesaka user 4.0K Mar 15 11:05 part_4_assembly

drwxr-xr-x  2 uesaka user 4.0K Mar 15 11:04 part_3_assembly

drwxr-xr-x  2 uesaka user 4.0K Mar 15 11:03 part_2_assembly

drwxr-xr-x  2 uesaka user 4.0K Mar 15 11:02 part_1_assembly

-rw-r--r--  1 uesaka user 214M Mar 15 11:01 data._part_4.fastq

-rw-r--r--  1 uesaka user 214M Mar 15 11:00 data._part_3.fastq

-rw-r--r--  1 uesaka user 212M Mar 15 11:00 data._part_2.fastq

-rw-r--r--  1 uesaka user 209M Mar 15 10:59 data._part_1.fastq

 

 

論文公開時はspade3.0を使っているが、現在ではspades3.11まで更新されている。現在のspadesと比較してみる。データは先ほどのコレラのtrimmed.fastqを使う。パラメータは上のテストと同じにした。spadesは以下のパラメータ設定でアセンブルを行った。

spades.py -k auto -t 30 --careful -1 reads_1.trimmed.fastq -2 reads_2.trimmed.fastq -o output

Questの分析結果。

f:id:kazumaxneo:20180315142207j:plain

このデータセットではspades(青線)の方がややcontiguityは上になった。

 追記

spadesのk-merを51に変えても変化がなかった。

追記

HGAのk-merサイズを長くするとcontiguityが悪くなった（*1）。

 

引用

HGA: de novo genome assembly method for bacterial genomes using high coverage short sequencing reads.

Al-Okaily AA

BMC Genomics. 2016 Mar 5;17:193.

 

*1

使えないと言っているわけではないことに注意してください。どんなツールもアルゴリズムの特性にあった使い方ができていなければパフォーマンスは発揮できません。