macでインフォマティクスmacでインフォマティクススマートフォン用の表示で見るmacでインフォマティクスNGS関連のインフォマティクス情報についてまとめています。PindelとPlatypusのフローを修正。freebayes、lumpyの誤りを修正。lumpyの流れを見やすく修正。構造変異検出ツールのパフォーマンスを比較したペーパーが出ている。data実際に導入して、パフォーマンスを比較してみる。以下のツールを紹介する。brewsci/scienceほとんどのスクリプトはインプットとしてソート済みのbamファイルを使う。indelとSNV検出編を確認してください。indel検出前にbamのリアライメントとBQSR補正を行い、できる限り正確なアライメントを行ったbamを出力する流れを書いています。R1.fqとR2.fqがfastqファイルで、ref.faがリファレンスとする。リファレンスのindexも必要とするものが多い。faidxでfasta.faiを作っておく。$fこれらの結果を1つのフォルダにまとめておき、各スクリプト実行時にそのフォルダを参照するようにする。追記--bamの作成、リアライメント、フィルタリング--bamの作成方法の例を示す。caller比較の論文で使われている方法になる(リンク)。ALN.BQSR.bamALN.BQSR.bamとその.baiファイルが得られる。callerの入力に使う。--分類--ショートリードから構造変化を予測する手法は大きく4つに分けることができる。<--)からindelを予測する。<--)部位にインバーションがあるだろうと考える方法論である。read(DPR)情報に基づき構造変化を予測する。アライメントのエラーもありうるので、エラーを考量して異常なペアが1箇所に2組以上集積している部位をコールしている。この手法には大きな問題点がある。indelは検出できない。インサートサイズの分布によって感度が変化する。構造変化後の配列に関する情報も得られないこともデメリットとなる。assemblyを組み合わせて検出する手法がいくつか報告されている。2017年現在Breakdancer-maxがダウンロードできるが、cpanでたくさんのperlモジュールをインストールする必要があって面倒だった。breakdancer")。解析手順は、まずbamファイルをサンプリングしてインサートサイズなどを分析したファイル(.cfgファイル)を作り、それからコールを行う。viewサンプリング自体は最初の数百行だけ行うらしく、すぐに終わる。してもランしたければ、オプション"-v"をdefaultの1から2-10まであげると一応ランはできる。breakdancer_output人のchr20で試すと解析は10分で終わり、500箇所ほどコールされた。(1-20bp)を予測するツール。リードを分割してアライメントし、そのアライメントのパターンから構造変化を予測している。MEMやRNA-seqで使うtophatなども使っている方式となる。mappingの代表的な論文なのでもう少し手順を説明する。1、ペアリードの片側だけマップされるペアを探索する。この片側のリードのアライメントが肝なので、ユニークでかつミスマッチがないパーフェクトマッチのアライメントが行われる。2、相方リードがマップされないのは、この相方リードが構造変化部位を横断するようにマップされているからと考える。この相方リードのアライメント部位を探すため、アライメントできた方のリード(アンカーと呼んでいる)からインサートサイズ平均の2倍の範囲内でマッピングされなかったリードの方がマッチする部位を探す。alignmnet)。3、5'側と3'側のアライメント部位にずれがあれば、即ちそのずれのサイズがdeletionのサイズになる(下図a)。5'側と3'側のアライメント部位にズレはないが、真ん中のピースが浮いた状態でマッピングされていれば、その真ん中のピースがinsert配列そのもののはずである。4、2つ以上のペアで同じ情報が検出された部位をコールする。片方のリードを正確にアライメントし、その相方の位置も限定することで偽陽性を抑え且つ解析をスピードアップしている。原理上リードの長さによって検出可能なindelの長さの上限が決まる。初期の30-40bpのリードでは20bpのinsertionが上限ということらしい。以下の図は解析の流れ。インストールの仕方が書いてある。macOSX環境ではインストールできなかったのでcentOSにインストールした。opengenomics/pindel---------------------デモランは以下の通りに行う。outputランに必要なものはリファンレス配列(.fasta)とそのindexファイル(.fasta.fai)。そのリファレンスにmapして作ったbamファイルとそのindexファイル(.bam.bai)。このあたりは他のツールと変わらない。pindelは情報を記載したコンフィグファイルを作って、それを元にランする。サンプルファイルたと、コンフィグファイルに相当するのはsimulated_config.txtファイルとなる。ファイル内には以下のことが書かれている(コピペして使うなら、スペース2つなど入れてるので、タブ1つに置換して下さい)。SAMPLE3データが1つなのでこれを以下のように変更した。SAMPLE1この600はインサートサイズ。output終わると構造変化の種類ごとにアライメントファイルが出力される。結果を見るにはpindel2vcfで変換する。pindel2vcfはビルドしてでできたpindel/に入っている。versionに準拠しているらしい。-eは必須パラメータではないが、デフォルトの1では高感度すぎて明確な変異なしのシミュレートデータ(50カバレッジ)でも1000近くの偽陽性コールが検出される。5になっていたのでここではそれに習った。目的やリードのカバレッジによってチューニングが必要と思われる。略語は以下のように説明されている。breakpoints最終的には以下のようなフローで全タイプのSVを検出できる。pindel.indel.vcfvcffilterがないなら、vcftlibsのレポジトリを--recursiveをつけてクローンする。#たくさんツールがある。ここではディレクトリ全体にパスを通す。approachesダウンロードサイト。インストールは、解凍したディレクトリの中で以下のコマンドを打つだけでよい。./install入力はbamファイルとなる。approaches.解析の流れはこのURLに丁寧に書いてある。macOSX環境でのインストールも可能。homebrewでインストールできる。erictdawson/lumpy-svソフトクリップのマッピングファイルから構造変化を予測する手法。memなどのマッパーを使う。デプスのないデータなどでもsensitivityが高いらしい。複数サンプルで差分をとって精度を上げるやり方も搭載している。テストランはtestディレクトリにあるtest.shを実行する。./test.sh問題なければNA12878、honeybee、riceのindel予測結果がvcf形式で出力される。解析時は下のような流れでコールする。sample.vcf解析できるまで手間取った。最後の肝心のlumpyexpressがモジュールを読み込めずエラーになる。最終的に強引にrootで進めることでなんとかランできた。samtoolsのソートでエラーになる人は、>(リダイレクト)とリダイレクトの後ろ側の.bamを決してランする。これはsamtoolsのバージョンによっては>と.bamが入らないバージョンがあるため。上の図に載っているように複数サンプルの差分や同時解析もできる。1、normalサンプルについて上のコマンドを実行し、SVに関係するリードのbamを得る。2、tumorサンプルについて上のコマンドを実行し、SVに関係するリードのbamを得る。multi_sample.vcflumpyはspeedseqパイプラインでも使われています。aln))、他のツールの結果もまとめてくれます。speedseqも検討してみてください(リンク)。Hosting.下のような流れでコールする。1500ただし結果はひどくしょぼい。進め方かパラメータを間違ってる可能性がある。data.deletionなどを予測する。コピーナンバーが変わるリピートのdeletionや、depthが少ないデータにも強いらしい。sampleなど)両方を使った差分解析を特徴とするが、単体のデータでも解析は可能。insertionはターゲットではないので注意する。こちらからダウンロードできる。centOSサーバーにインストールした。基本的にはgitからダウンロードしてmakeするだけ。