シーケンシングデータのハプロタイプを可視化し、リードを分類する 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）。このアプローチは、バリアント密度に依存しないが、ローカルまたはグローバルのハプロタイプを推論することはできない。１サンプル内で同じバクテリア種に複数の株がある場合、ゲノム分析を行うベストな戦略を決定するためには、シーケンシングデータのハプロタイプを視覚化するツールが必要である。

　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で動作確認した。

解析の流れ

１、リードをアライメントしてbamを出力する（公式例）。

２、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

３、 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としているが、拡張子なしでも認識する）。

４、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

公式より転載

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

　右の図は、バリアント組み合わせに基づいてユニークなバリアントプロファイルを"flow"で表している。先ほどの４リードサポートがあるオレンジの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