macでインフォマティクス

macでインフォマティクス

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

シーケンシングデータのハプロタイプを可視化し、リードを分類する HapFlow

2018 11/3 誤字修正 

2019 3/18 freebayes追記

 

 ハイスループットシーケンシング技術の出現により、バクテリア集団のシーケンシングのような新しい実験的アプローチが可能になった。感染は、しばしば同じ種の複数の株を含んでおり(Darch et al、2015; Taylor et al、1995)、これは治療方法決定(Cohen et al、2012)に重要な意味を持つ。mixed-strain の集団を分析するためのいくつかの方法が開発されている。ShoRAH(Zagordi et al、2011)は、全体的なハプロタイプの最小セットを再構成し、推測されたハプロタイプの頻度を推定する。これは、オーバーラップしたリードによってリンクするため十分な密度のバリアントが必要になる。ドミナントの株およびマイナーな株からの感染を同定するために 2ステップの最尤アプローチも報告されている(Eyre et al、2013)。このアプローチは、バリアント密度に依存しないが、ローカルまたはグローバルのハプロタイプを推論することはできない。1サンプル内で同じバクテリア種に複数の株がある場合、ゲノム分析を行うベストな戦略を決定するためには、シーケンシングデータのハプロタイプを視覚化するツールが必要である。

 Savant(Fiume et al、2010 pubmed)、Tablet(Milne et al、2010)、Consed(Gordon and Green、2013 pubmed)など、多くの優れたリードアライメント視覚化ツールが存在する。これらのツールは、各リードが線、またはベースの行として表される方法でリードを配列する。このレイアウトは、バリアントまたはミスアライメントされたリードを識別するのには満足できるものだが、リードに存在するハプロタイプを識別するのには理想的ではない。これは、リードは密接にまとめられているため、離れた位置のバリアントが同じリードペアにあるかどうかを判断することは困難であるためである。さらに、ハプロタイプによってグループ分けされていないので、シーケンシングデータでハプロタイプがどの程度頻繁に表されるかを同定することが困難である。

 HapFlowは、存在するハプロタイプをより容易に同定するため、リードアラインメントデータを抽出することによってこれらの問題を解決する。 HapFlowは、シークエンシングアーチファクトの同定、サンプル中に存在する最小の株数の同定、およびシーケンシングデータ単独でローカルハプロタイプまたはグローバルハプロタイプの決定が可能かどうかを判断するのに役立つ。

 

HP(プログラムの他、テストデータ、マニュアルもダウンロードできる)

http://mjsull.github.io/HapFlow/files.html

wiki

https://github.com/mjsull/HapFlow/wiki

 

インストール

HPからダウンロードする。CUI環境で動くpythonのパッケージと、GUIバージョンがある(windows版、mac版、linux版あり)。

HapFlow - downloads

ここではmac osxで動作確認した。  

  

解析の流れ

1、リードをアライメントしてbamを出力する(公式例)。

 

2、Freebayesなどのploidyを指定してバリアントコールできるツールを使い、バリアントコールを行う(freebayesはcondaで導入できる)。

例えばploidyの値を3に増やしてバリアントコールする。

freebayes -f ref.fa -p 3 -F 0.2 -C 5 input.bam >aln.vcf

#もっと厳しくするなら
freebayes -f assembly.fasta -p 3 -F 0.2 -C 50 -q 30 -m 30 --min-coverage 100 -i -u input.bam >aln.vcf
  • -F     Require at least this fraction of observations supporting an alternate allele within a single individual in the in order to evaluate the position. default: 0.2
  • -p    Sets the default ploidy for the analysis to N. default: 2
  • -C    Require at least this count of observations supporting an alternate allele within a single individual in order to evaluate the position. default: 2

 

3、 Flow Fileを作成するため、HapFlowを実行する。

./HapFlow

コマンドを叩くとHapFlowのウィンドウが出現する。

f:id:kazumaxneo:20181102004906p:plain

メニューバーのCreate flow fileをクリック

f:id:kazumaxneo:20181102212815p:plain

ウィンドウが表示されるので、bamとvcfを指定する。

f:id:kazumaxneo:20181102212451p:plain

ジョブ開始はOKをクリック。ランが終わると.flwファイルが出力される(上の画像ではoutoutをflowとしているが、拡張子なしでも認識する)。

 

4、File => Load flow file  から.flwファイルを読み込み可視化する。

f:id:kazumaxneo:20181102213129p:plain

ここで読み込んでいるのは、HPからダウンロードしたテストデータ1のtutorial_1.flw(real data from mixed-strain infection of Koala)

 

図の表現方法は、公式wikiで丁寧に説明されている。 

https://github.com/mjsull/HapFlow/wiki/Interpreting-data-using-HapFlow

f:id:kazumaxneo:20181103062406p:plain

公式より転載

 左の図は、リードがリファレンスゲノムにアライメントされていることを表している。縦の線は、そこにバリエーションがあることを意味する。同じポジションに赤と青のバリアント(アレル)があり、それが3ポジション存在している。リードによってバリアントのパターンは異なるが、左から順にのバリアントプロファイルを持つリードが一番多い。左図右端の縦線は、バリアントのプロファイルによるグループ分け結果を表していて、こののバリアントプロファイルのリードはオレンジの縦線割り当てられている。

 右の図は、バリアント組み合わせに基づいてユニークなバリアントプロファイルを"flow"で表している。先ほどの4リードサポートがあるオレンジのflowが一番太く、ーの組み合わせの紫のflowが一番薄い。このように、flowの厚みはそのユニークなバリアントプロファイルの比率を表している。flowの矢印の向きはリードの向きと対応している。

 

次の図

f:id:kazumaxneo:20181103070400p:plain

公式より転載

 ペアエンドシーケンシングの場合、ペアのリードは同じテンプレートに由来しているので、より長い範囲でバリアントプロファイルをカバーできる。flow図では、ペアエンドリードは上図のように破線で接続される。そのほか、ゲノムは一番上に青のboxで表わされ(右の図)、バリアントのリファレンス上での位置は、青boxの中のオレンジboxで表わされる。オレンジ領域を拡大したのが一段下のオレンジboxで、バリアント位置はオレンジbox中の縦線で表わされる。

この図のflowをよく見ると、dark green flowの左端矢印は長さが半分になっているが、これはバリアント位置がリードの上流末端に相当することを意味している。同じく、リード下流末端のバリアントも矢印の長さが半分になる(Light greenのflowの右側矢印)。末端バリアントのみ区別した表記があるのは、リード末端はミスアライメントに由来することが多いためと述べられている。flowはバリアントを共有している度合いによってグルーピングされ、にたバリアントプロファイルを持つflowは同系色の色で表現され、近傍にまとめられる(上で言えば、赤とオレンジのflowグループ、濃い緑と薄緑のflowグループ)。

 

 

f:id:kazumaxneo:20181103112200j:plain

 一番下には、リファレンスのポジションとリファレンスのアレル(*付きの方)が表示される。

 

f:id:kazumaxneo:20181103112817j:plain

この領域では、左端2番目の363598のアレルでメジャーなもの(青色のflows) は、次(左端3番目)の363732ではマイナーなアレルに変化している。

 

view > hide gappedで、ギャップのラインを非表示。

f:id:kazumaxneo:20181103113649j:plain

また、A、D、W、Sそれぞれのキーを押すと、flowの幅と高さを1段階ずつ変更できる。

Dを数回押して横長にした。

f:id:kazumaxneo:20181103114310j:plain

 

tools > goto baseでポジションにジャンプできる。

f:id:kazumaxneo:20181103114502j:plain

 

 

異なるtaxonomicグループが混ざっていることがわかれば、taxonomicグループごとにリードを分割することもできる。そのためには、自分でリードをアサインしていく必要がある。

左端から順にダブルクリックして、OTUの番号を割りふる。まずはOTU1をアサインする。左端から順に同じ系統のカラーをダブルクリックしていく(1つでも飛ばすと番号がつかない)。

f:id:kazumaxneo:20181103120041j:plain

全て終わったら、 OTU2に切り替え、同様に番号をアサインしていく。

f:id:kazumaxneo:20181103120246j:plain

 

アサインが終わった状態

f:id:kazumaxneo:20181103120333j:plain

公式より

 

tools > Write OTUsで、アサインしたグループごとにbamを分割出力する。

f:id:kazumaxneo:20181103120422j:plain

 

出力された。

f:id:kazumaxneo:20181103120501j:plain

 

少し大きなゲノムでも表示できるようです。詳細はwikiのperfrmanceを読んで下さい。質問も受け付けているようです。

https://github.com/mjsull/HapFlow/wiki/Performance

引用
HapFlow: visualizing haplotypes in sequencing data
Sullivan MJ, Bachmann NL, Timms P, Polkinghorne A

Bioinformatics. 2016 Feb 1;32(3):441-3