macでインフォマティクス

macでインフォマティクス

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

バクテリアのsub-populationsのレアバリアント検出ツール VarCap

bacteria small indel structural variations (SV) SNV sub-populations 2017

 8/20 sambambaコマンドミス修正、varscan2バグに対応する迂回コマンド追加

 

　1つの原核生物種のheterogeneityな集団の遺伝子型決定（genotyping）は、一般的な選択圧下での集団（populations）組成および動態に関する微生物学的問題に対処するためにますます重要な方法になってきている。 このアプローチは、例えば、実験室進化（experimental evolution : EE）実験（ Barrick＆Lenski、2013 ）および宿主 - 病原菌システムの研究（ Gardy et al、2011 ; Bos et al、2011 ; McElroy、Thomas＆Luciani、 2014 ）。 次世代シークエンシング（NGS）技術の最近の進展により、短期間の高いカバレッジでのシーケンスが可能になったが、リード長が短いという制限がある。

　ショートリードシーケンシングからゲノムを構築する古典的アプローチは、最も豊富な genotypeをゲノムコンティグおよびscaffoldsに再構築することが好ましい。 ハプロタイプ頻度情報を検索するためには、アセンブリまたはリファレンスゲノム上にリードをマッピングする必要がある。 次に、得られたアラインメントに対してバリアントコールを実行する。 予測されたバリアントは、バリアント部位の不十分な連結のために全ハプロタイプ再構成が可能でない場合、ハプロタイプまたは対立遺伝子にフェージングすることができる。 しかし、変異予測は、InDelsや置換などのシークエンシングエラーによって誤検出につながる可能性がある。 リードは、その短さのためにマッピング中に誤ってアライメントされる可能性があり、したがって、偽陽性バリアントコール（ Li、2014 ）につながる可能性がある。 シークエンシングエラーは、クオリティフィルタリングとエラー訂正によって部分的に削減できる（ Yang、Chockalingam＆Aluru、2013 ）。 結果として、イルミナの置換エラー率は1％以下に低下する可能性がある。InDelホモポリマーのエラーはストレッチの長さに比例して対数的に累積することが示されており、それによって確実に識別することができる。

　進化する集団では、アレルのheterogeneityが予想される（ 図1 ）。 これまでの原核生物の遺伝子型解析研究のほとんどは、クローナルなバクテリアカルチャーのリシーケンシングによって行われてきた（Maharjan et al、2013 ; Blount et al、2012）。主に、 Pool-seqと名付けられた非クローナルな集団のディープシーケンシングの手法が、メタゲノムコミュニティのプロファイリング（ Qin et al、2010 ）と、わずかであるが、アレル頻度の特徴付けのため延長されて使用されている（ Eyre et al、 2013 ; Khan et al、2011 ; Köser et al、2012 ; Pulido-Tamayo et al、2015 ）。 しかし、heterogeneityな集団における非クローナルなバリアントのgenotypingは依然として困難である（ DePristo et al、2011 ; Nielsen et al、2011 ; Kofler＆Schlötterer、2014 ; Pulido-Tamayo et al、2015 ）。

　異なるハプロタイプまたは対立遺伝子頻度の最も完全な画像を得るためには、Pool-seqを高いカバレッジで行い、SNP、InDelsおよび構造変異（SV）の全てのタイプのバリアントを検出することが基本になる。 全てのタイプのバリアントを検出するには、異なるアプローチを利用する複数のバリアントコールソフトウェアツールを統合することが有効である。

　SNPを特定するためによく使用されるツールは、SAMtools / bcftoolsとGATKである（ Li et al、2009 ; McKenna et al、2010 ）が、これらのツールは二倍体生物のバリアントを検出するという前提で開発されており、haploidの原核生物の検出能力を制限している。 よって著者らは、低カバレッジや高カバレッジのデータのSNP頻度を予測できるより一般的なツールである、VarScan2（ Koboldt et al、2012 ）、LoFreq-Star（ Wilmら、2012 ）、Breseq（ Barrick et al、2014 ）（紹介）、FreeBayes（ Garrison＆Marth、2012 ）を考えた。 これまでに報告された評価から、ここでは、Breseqよりも感度の点で優れているLofreq-Starを使用した（ Wilm et al、2012 ）（Lofreq紹介）。また、既知のプールサイズの真核生物データのPool-seq実験で広く使用されるFreeBayesを評価したが、 これらのツールは、すべてアライメントファイルまたはmpileupファイルで動作し、リードとマッピングのクオリティスコアとストランドバイアスのフィルタを使用してSNPを確実に検出することがわかり、プールサイズが未知のバクテリア集団の分析に使えることができると考えた。 さらに、SAMtools / bcftoolsとVarScan2とFreeBayesを使用して、小さなInDelを識別することもできる。 Pindel（ Ye et al、2009 ）はpattern growth アルゴリズムを使用して、1 bpから10 kbまでの小さなInDelと大きいInDelを検出できる。 BreakdancerとDelly（ Chen et al、2009a ; Rausch et al、2012 ）は、インサートサイズの偏差を利用しているため、大規模なInDelと転座、duplication、逆位などの300bpを超える構造変化（SV） バリエーションを見つけることができる。比較するための代案手法として、マップされたリードに頼らず デノボアセンブルされたコンティグを使用するCortex_var（ Iqbal et al、2012 ）がある。（一部略）。これらのアプローチは全て、二倍体ゲノムから低頻度のgenotypeを有するmultiploid populationsに至る様々な度合いの接合性（zygosity）について設計されている。

　実験室進化実験における原核生物集団のgenotypingは、通常、高いカバレッジを有する多くのNGSデータセットに基づいたものである。 したがって完全自動化されたマッピングとバリアントコールのソフトウェアが求められている。それは、感度が高く且つ正確であり、少ないsubpopulationsを検出し、またすべてのタイプのバリアントを検出できるものでなくてはならない。 著者らの知る限り、このようなソフトウェアのワークフローはこれまでに公開されていない。 この研究では、合成データを使ってバリアントコールを評価し、感度と精度を比較した。 これにより、原核生物集団の正確かつ迅速なgenotypingのためのワークフローVarCapを開発し、検証することができた。 最後に、アメーバのバクテリア共生（bacterial symbiont of amoebae）の長期に渡る実験室進化実験にVarCapを適用した。

 

VarCap Github


Methods

シミュレーションデータ

  • バリアント検出対象のデータセットは、NCBI Refseqデータベース（ Pruitt et al、2012 ）からダウンロードしたリファレンスゲノムに異なるタイプのバリアントを無作為に挿入することによって作成された（ table S1 ）。バクテリアゲノムのSNP / InDel比がしばしば15〜20と報告されたので、小さなInDelについてはSNP / InDel比=10を使用し、すべてのInDelについては20を使用した（ Moran、McLaughlin＆Sorek、2009 ; Chen et al、2009b ）。 またラージサイズのInDelも導入した。これは、水平遺伝子伝達のプロセスを模倣するためである。 構造変化がバクテリアゲノムの進化にとって重要であると報告されているので、検出ソフトウェアにはチャレンジングだが、転座、重複および逆位もわずかに導入した。
  • 135の様々なタイプのバリアントを含むデータセットと、特定のタイプのバリアントだけを50-100含むデータセットを作成した。 混合タイプのデータセットの135の変異型は、100のSNP、10のsmall InDel、10のlarge InDelおよび5の転座、5個のduplication（1回の二重複製を含む）、および5の逆位からなる（これをsim_135VARとする。 table S1 ）。
  • 100のSNPは、単一のSNPおよび変異のホットスポットとして配置した。ホットスポット内のSNPの最大数は4であり、ランダムに4〜60塩基の距離内に配置した。 large InDelのサイズは、5〜2,000ヌクレオチドの間で無作為に選び、転座、duplicationおよび逆位のサイズは300〜2,000ヌクレオチドの間で無作為に選んだ。
  • ALFSim（pubmed）はゲノム進化シミュレーターであり、より遠くに進化したsub-populationsの進化をシミュレートするために使用された（ Dalquen et al、2012 ）。 ゲノムアノテーションからコーディング及び遺伝子間ヌクレオチド配列を抽出し、ALFSimの入力とした。 ALFsimの出力から、全ヌクレオチドに0.8％違いがあるシミュレーションのsubspeciesを選択し、その結果、21,000のSNP、100のInDelおよび3のduplicationが得られた。 得られたfastaファイルは、シーケンシングリードのシミュレーション、heterogeneityな集団の構築、およびバリアント予測に用いられた。
  • 100ヌクレオチド（nt）のシーケンシングリードのシミュレーションにSimSeq（ https://github.com/jstjohn/SimSeq ）を使用した。 リードは、250ntのインサートサイズおよび10,20および30％のインサートサイズSDでシミュレートした。 変異したリファレンスからシミュレートしたリードと元の変異なしリファレンスからシミュレートしたリードを混合することによって、minor allele frequencies（MAF）（wiki）をシミュレートした。 MAF 40, 20, 10, 4％で作成した。
  • シミュレートされたリードのクオリティは、FastQCで分析し、 Prinseq-liteおよびTrimmomaticによってフィルタリングおよびトリミングが行われた。サイズ10のスライディングウインドウで3 '側からクオリティスコアが20以下の塩基部分を除去し、40 n以下になった場合リードを除去した。 それから標準設定でbwa-mem（bwa-0.7.5a、 Li、2013 ; Li＆Durbin、2009 ）を使用してマッピングし、bamファイルにした。 SAMtoolsとPicard Toolsを使用して、SAMからbamファイルへ変換し、Cortex_varの入力はこのbamから変換したfastqファイルを使った。

リアルデータ

  • Amoeba symbiont Protochlamydia amoebophila（サンガーシーケンスでゲノム決定済み（Horn et al。、2004））の実験室で約9年培養した２つの独立したカルチャーをイルミナGAIIプラットフォーム（100bp PE、250bpインサート、3,000xカバレッジ） SRA：SRR5123091）でシーケンシングした。 分析のために、得られたデータをそれぞれ250倍のカバレッジにランダム抽出して分割し、異なる存在量のsub-populationsのバリアントを検出するために利用した。

バリアントコール

  • True positiveおよびfalse positiveを評価するために、40％、20％、10％、5％および2％の存在量でSNP、InDelおよびSVを含有する人工のnon-clonal なpopulationssをシミュレートした。SNV検出にSAMtools / bcftools（0.1.18、 Li et al、2009 ）、GATK-lite（Genome AnalysisTKLite-2.2-8, McKenna et al., 2010 ）、VarScan2（2.3.6、 Koboldt et al、2012 ） 、LoFreq（0.6.1、 Wilm et al、2012 ）およびLoFreq2（lofreq-star 2.0.0 beta 1、 https://github.com/CSB5/lofreq ）を使った。 Small InDel検出には、VarScan2、Pindel（024t、 Ye et al、2009 ）を使用した。 Large InDelと構造変化（SV）については、1〜1,000ntの間でよく機能するPindelと、breakdancer（breakdancer-1.1_2011_02_21、 Chen et al、2009a ）とdelly（0.0.11、 Rausch et。その他、2012 ）を使った。 さらに、アセンブルされたコンティグをと比較することによって変化を検出できるcortex_var（CORTEX_release_v1.0.5.14、 Iqbal et al、2012 ）を使用した。 複合ワークフローの感度と精度は、recall = TP /（TP + FN）、precision = TP /（TP + FP）として計算した。 TP、FPおよびFNはバリアントごとに測定し、例えばSNPおよび大きな欠失は同じウエイトを与える。

バリアントの最小存在量の設定

  • 一部のバリアントコーラーの検出には、4〜8のリードにバリアントが存在していることがコールの条件になっている。ここではminimum absolute abundance （MAA）として8リードを設定する。 カバレッジはゲノムに沿ってわずかに変化するため、MAAに加え、Minimum Relaxy Abundance（MRA）も定義した。MRAは、全カバレッジと比較した変異リードのパーセンテージを意味する。 8リードのMAAは、400×カバレッジでは2％のMRAに対応する。

 

結果１　バリアントコールベンチマークによる適したツールの絞り込み

　現時点で、Procaryotesゲノムの全タイプのバリアントを同時に検出できるソフトウェアツールは存在しない。 そのため、著者らはSNP、InDelおよび構造変異（SV）の検出ツールを別々に評価した。最初に、最も基本的な変異モデル（Methodsのsim_135VAR）を利用して、非常に高感度で確実にレアバリアントを検出できるツールを評価した。その結果を元に、sub-populationsのバリアント検出に適したツールを選抜した（論文　図2）。

（各ツールのSNPs、indel、SV検出パフォーマンスの詳細は省略）

 

結果２ false positiveへの対策

　VarCapは、SNP検出にLoFreq-StarとVarscan2、small indel検出にVarscan2とPindel、を使用する。理由は、方法が異なるが感度が同じほど高かったためである。より大きなSVについては、Pindelの使うpattern growth、スプリットリード、ペアエンドリード情報のコンビネーションアプローチが高い感度をもたらすことを観察した。この方法は、規定された限界（1〜10kb）内で良好に機能する。より長いバリアントを検出するには、ペアドエンド情報のみを使用するBreakdancerが有効である（検出解像度を犠牲にして）。Cortex_varは感度が劣っていたが、de-novoアプローチを使用することにより、検出されたバリアントに関するより多くの情報を明らかにしていた。この情報は、バリアントのタイプ、位置、長さまたはシーケンスを正しく識別するために使用できる。以上の観察から、large InDelsにはPindel、Breakdancer、Cortex_varを使い、SVにはBreakdancer、Delly、Pindel、Cortex_varを使用する。

  • False positives（FP）は、通常、1％の割合以下で存在するシークエンシングエラーによって発生するため、FPはこの割合以下の発生が期待される。シーケンシングエラーの影響を調べるために、MRA 2％と1％に焦点を当てた。1％のMRAでは、ほぼ全ての変異が検出できるようになったが、SNPsのFPは大きく増加し、1Mbあたり80のFPが観察された（表1）。一方、他のタイプのバリアントのFPは1FP / Mbで、ほぼ変化なしだった。これらの結果は、偽陽性SNPsがシークエンシングエラーによって引き起こされることを明確に示している。

  • MRA 2％と1％では、小さなInDelのFPのコールがリピート内からのものであることを見出した。リピートの類似性のためにリードがマッピングされ、そこでFPコールを引き起こす。

    FP数がMAAとカバレッジにどのように影響されるかを評価するため、シミュレーションのsim_135VARデータセットで、80〜1600×カバレッジと4〜20MAAにそれぞれ変化させてFP数を調べた（図5）。図では、MRAのカットオフ（全カバレッジと比較した変異リードのパーセンテージ）に含まれるセルが色で表現されている。<MRA1なら1600xのFPは全て除かれる。160xでは ≤MRA2でもFPが除ききれない。このことから、著者らは、FPを除くには、MRAカットオフに加えてMAAカットオフを使用して、低いカバレッジでもFPコールを回避する必要があることわかった、と書いている。

  • リードのミスマッピングがFPの１つの原因として報告されている（Li、2014）。 したがって、リファレンスゲノムが不完全だと、類似の領域にマッピングされたリードがそこでFPをコールする可能性がある。人工的に短くしたリファレンスゲノムに変異のないのないリードをマッピングし、MRA2％でバリアントコールすると、1Mbあたり∼180 FP SNPs/75 FP検出された。これらは 存在量が異なり（20％、8％、3％）、ホットスポットにグループ分けされた（図６A）。正しいリファレンスにマッピングすると偽陽性変異は観察されなかった（図６B）。 すなわち、誤ってマッピングされたリードがFPバリアントコールを引き起こすという著者らの仮定を強く支持している。
  • FP検出を防ぐためにバリアントサイトのカバレッジ情報を使用するが、カバレッジ情報だけでは、低頻度FPの解決には完全ではない（図６Cベン図の"COV"）。すべてのFPは小さなクラスターまたはホットスポットとして存在していたので、インサートサイズの2倍の長さのスライディングウインドウサイズで4つ以上のSNPを含む領域を探し、FPを引き起こす可能性のあるブレイクポイント（BP）を識別した。これらのフィルタを適用することで、FPコールを識別し除外することができた。
  • 逆位を詳細に見ると、それらはほとんど逆位であると特定されてなかったことに気づいた。しかし、逆位の開始位置および終了位置はブレイクポイントとしてマークされていた（表S2）。ブレイクポイントになっていたのは、ペアリードの片方のみがマッピングされ、順方向または逆方向のリードのみが蓄積されるために発生する。これらは、より大きな構造変化のより一般的な指標である。したがって、これらのコールはバリアントの部分的な解決を表してる。
  • 最終的にFPを特定、除外するために、マッピング後に次のフィルタリングを適用する：シーケンシングエラーによって引き起こされるFP のSNPコールを回避するために、MRAを2％で適用する。リピート領域にマッピングされたリードによるFPを回避するために、ほぼ同一の領域を予めマスクする。真のバリアントのタイプの不完全な検出によって引き起こされるFPを解決するために、より小さなバリアントより大きなバリアントを優先する。したがって、小さなバリアントはより大きなバリアントにアサインする。
  • 精度と堅牢性を得るために、バリアントあたりの予測のintersectionを使う必要があった。したがって、少なくとも2つの異なるツールでバリアントコールをサポートするものを選抜した。このステップにより、低MRAカットオフ（1％）での精度が上がったが、検出率はわずかに低下した（表1）。この発見は、最近の論文でも支持されている（pubmed）。（inversionsとlarge insertionsはこの条件を使っていない）

 

結果３　異なるゲノムでの検出率

  • ゲノムのG + C含量やサイズなどの異なる特性により、バリアントコールの感度と精度が影響を受ける可能性があるので、6つの異なるゲノムでワークフローが評価した。これらの生物は、26〜72％のGC、0.68〜8.66Mbのゲノムサイズを有する5つのバクテリアおよび1つのアーキアから構成されている。MRA2％またはMAA 8リードで解析すると、これまでの結果と一致し、シミュレーションされたバリアント（> 90％）のほとんどが検出された。MRAとMAAの違いとして、2％のMRAフィルタリングではGCまたはゲノムサイズに依存しないが、MAA８リードフィルタリングで高いGC含量と大きなゲノムで変異検出感度が低下した（図7）。このため、ゲノムの特性に影響されないMRAを変異検出の最小カットオフとして使用するというのは妥当であることが再び裏付けられた。これはカバレッジの低い領域に固定MAAが必要であることを排除するものではない。

 

結果４　異なる進化をしたpopulationsのバリアント検出パフォーマンス

　ポジティブ選択の下でより遠くに進化したpopulationsが出現すると、そのpopulationsはより多くの変異を蓄積する。ALFSim（Dalquenら、2012）を用い、P. amoebophilaゲノムにより遠い進化的変化（SNP、InDelsおよびduplications）を加えシミュレートした。 進化したゲノムは、およそ21,000のSNP、100のInDelおよび3つの遺伝子重複を含み、オリジナルリファレンスとの類似性が約99％になった。

  • VarCapによる変異の感受性を、4％の低量の sub-populationsで評価した。 MRAは3％、2％、1％、MAAは8回リードだった（400×で2％のMRAに相当）。最低必要量のカットオフに応じて、全てのSNPの90％〜99％、InDelの74％〜94％、3つのduplicationsのうち2つ検出できた。MRAを3％から2％に下げると、SNPの真の陽性検出率は98％に上昇し、偽陽性率は0.3％以下のままであった。しかし、MRAをさらに1％に下げると、FPが400 FP / Mbに増加したが（2％の数十倍）、TPレートは99％に増加した（図8A）。 MRA2％では、FPなしに全てのsmall InDel（サイズ= 1）を含むすべてのInDelの90％以上を検出できた（図8A）。重複に関しては、MRAを問わず2つのMRAを見つけられたが、最短のMRAは欠落していた（図8C SV（DUP））。これらの発見は、進化したゲノムがかなり異なる場合であっても高い精度を達成できることを意味している。しかし、新しい気づきは、最近のduplicationsが反復領域に類似しているため、誤ってアライメントされたリードによるFP検出につながる可能性があるということである。そのため、FPの可能性が領域として重複領域のタグ付けをワークフローに組み込んだ。

 

結果５　lab-evolutionのリアルデータへの適用

  • 進化するpopulations内の変異頻度を予測するために、構築した変異コールワークフローを使い、P. amoebophilaの長期培養株を調べた。 20％から2％までの異なるMRAカットオフを使用し、2％までの頻度で変異を明らかにした（図9A、外リング）。 2％のMRAで、34のSNP、20のInDelおよび17のSVからなる合計71のバリアントを検出した。 SNPおよびSNPは、SNPEffを用いてアノテーションを付けた。変異の約83％がコード領域内に位置していた（表S3）。
  • VarCapのコールのバリデーションのために、4％、11％、28％の豊富さで予測された3つのバリエーションを追試した。3つのバリアントコールを含む領域それぞれPCRで増やしてベクターにクローニングし、各々16クローンをSangerシーケンシングした（表2、ファイルS1-S2）。 3つのバリエーションをすべてを検出でき、VarCapソフトウェアの予測が正しいことを確認できた。

Discussionより

In order to overcome false positives by misplaced reads, we removed variants that at least fulfill two of the four following rules: 

  • (I)  Either variants lie within regions with a coverage of 20% above the average
  • (II)  If there is a break position detected at or within read length of the variant site
  • (III)  If they lie within a repeat region
  • (IV)  If more than five variants lie within the length of one insert size. 

The efficiency for FP removal for each rule may differ among experiments as they dependent on organism, experiment setup, sequencing and reference quality. Therefore, we strongly suggest to use all rules in combination for a most flexible removal of FP predictions due to misplaced reads.

 

テスト 

Galaxyにアクセスしてログインし、このリンクのマニュアル:

http://galaxy.csb.univie.ac.at:8080/static/welcome.html

に従いVarCapを検索したが見つからなかった。VarCapのフィルタリングツールのみ検索から見つかるが、説明と違う気がする。GalaxyのMainサーバー以外にあるのか不明。

ここでは、別個にツールをランして VarCapのフィルタリングツールのジョブを走らせてみる。

f:id:kazumaxneo:20180819150423p:plain

 

１、delly

各SVタイプを順番にコールする。マージコマンドもあるが、ここでは実行しない。  

delly call -t DEL -o delly_del.bcf -g ref.fa pairend.bam 
delly call -t DUP -o delly_dup.bcf -g ref.fa pairend.bam 
delly call -t INV -o delly_inv.bcf -g ref.fa pairend.bam 
delly call -t BND -o delly_bnd.bcf -g ref.fa pairend.bam 
delly call -t INS -o delly_ins.bcf -g ref.fa pairend.bam

bcftools view delly_del.bcf > delly_del.vcf 
bcftools view delly_dup.bcf > delly_dup.vcf 
bcftools view delly_inv.bcf > delly_inv.vcf 
bcftools view delly_bnd.bcf > delly_bnd.vcf 
bcftools view delly_ins.bcf > delly_ins.vcf

 

２、LoFreq（紹介）

lofreq call -f ref.fa -o lofreq.vcf pairend.bam

 

３、 varscan2（wiki）

samtoolsのpileupの互換で高速な sambambaを使う（紹介）。

sambamba mpileup pairend.bam -t 12 --samtools -f ref.fa > pileup
#またはsamtoolsでパイルアップするなら
samtools mpileup -f ref.fa pairend.bam > pileup 

#varscanのpileup2snpにエラーがありvcfヘッダが出力されない。代わりにmpileup2vcfを使う。
varscan mpileup2vcf mpileup --output-vcf 1 --variants 1 \
--min-coverage 8 --min-reads2 2 --min-avg-qual 15 --min-var-freq 0.01 \
--min-freq-for-hom 0.75 --p-value 99e-02 \
> varscan.vcf


#フィルタリング1 SNV
vcftools --vcf varscan.vcf --remove-indels --recode --recode-INFO-all --out varscan_SNV.vcf
#フィルタリング2 indel
vcftools --vcf varscan.vcf --keep-only-indels --recode --recode-INFO-all --out varscan_indel.vcf


#エラーがないなら、varscanのSNVコールは
VarScan pileup2snp mpileup --output-vcf 1 \
--min-coverage 8 --min-reads 2 --min-avg-qual 15 \
--min-var-freq 0.01 --min-freq-for-hom 0.75 \
--p-value 99e-02 \
> varscan_snv.vcf

#エラーがないなら、varscanのsmall-indelコールは
VarScan pileup2indel mpileup --output-vcf 1 \
--min-coverage 8 --min-reads 2 --min-avg-qual 15 \
--min-var-freq 0.01 --min-freq-for-hom 0.75 \
--p-value 99e-02 \
> varscan_indel.vcf
  • --min-coverage Minimum read depth at a position to make a call [8]
  • --min-coverage Minimum read depth at a position to make a call [8]
  • --min-reads2 Minimum supporting reads at a position to call variants [2]
  • --min-avg-qual Minimum base quality at a position to count a read [15]
  • --min-var-freq Minimum variant allele frequency threshold [0.01]
  • --min-freq-for-hom Minimum frequency to call homozygote [0.75]
  • --p-value Default p-value threshold for calling variants [99e-02]
  • --strand-filter Ignore variants with >90% support on one strand [1]
  • --output-vcf If set to 1, outputs in VCF format
  • --variants Report only variant (SNP/indel) positions (mpileup2cns only) [0] 

"--min-reads2 2"にすると非常に高感度になる。

 

４、Pindelでsmall indelとSVコール（リンク）。dokcer imageを使う。

docker pull opengenomics/pindel

#現在のディレクトリとdataを共有dirにしてラン
docker run -m "50g" --rm -itv $PWD/:/data opengenomics/pindel

#configファイル作成 
printf "pairend.bam_400_sample1" > config
sed -e 's/_/ /g' config > pindel_config && rm config #<space>を入れたかったので少し遠回りした

pindel -f ref.fa -T 16 -i pindel_config -o pindel
cat pindel_* > pindel_merge

#BP call（BP = unassigned breakpoints）
pindel2vcf -R ref.fa -r ref.fa -p pindel_BP -v pinde_BP.vcf -e 4 -d 20180819
#LI call（LI = large insertion）
pindel2vcf -R ref.fa -r ref.fa -p pindel_LI -v pinde_LI.vcf -e 4 -d 20180819 pindel2vcf -R ref.fa -r ref.fa -p
#merge call
pindel2vcf -R ref.fa -r ref.fa -p pindel_merge -v pindel_merge.vcf -e 4 -d 20180819


#参考 mergeからSVタイプを絞り込むならvcffilterを使う
vcffilter -f " SVTYPE = INS | SVTYPE = DEL " pindel_merge.vcf > pindel.indel.vcf

 

５、breakdancerでSVコール（リンク）。dokcer imageを使う。

docker pull molecular/breakdancer
docker run --rm -itv $PWD:/data/ molecular/breakdancer

cd /data/
bam2cfg.pl pairend.bam > breakdancer.cfg
breakdancer-max breakdancer.cfg > breakdancer_output

 

６、cortex_varでバリアントコール（GIthubリンク）。(commnad)

#インストールしてないなら導入(centos7に入れた)
git clone https://github.com/iqbal-lab/cortex.git
cd cortex/
bash install.h
make cortex_var  #binにcortex_varのバイナリ（k=31）ができる。

./cortex_var_31_c1 --format FASTQ --pe_list pair1.fq.gz,pair2.fq.gz \
--kmer_size 31 --mem_width 100 --mem_height 15 --max_read_len 100 \
--output_supernodes contigs.fa

cortexでエラーが起きる。osを変えてやり直してみる。

 

 

引用

Variant profiling of evolving prokaryotic populations
Markus Zojer, Lisa N. Schuster, Frederik Schulz, Alexander Pfundner, Matthias Horn, Thomas Rattei

PeerJ. 2017; 5: e2997.Published online 2017 Feb 16.

 

 

参考

varscan

https://www.biostars.org/p/236408/