~indel.vcfが出力ファイル。thread数bedファイルを指定することで解析領域を限定することも可能。マニュアルには300bpまでのinsertionをコールできるとあるが、検証したところ確かに300bpまでの長さのInsertionを検出できていた。Novoalignなどの多くのアライナーをサポートしている。brewやcondaでインストール可能。varscanVarscanはpileupデータを使うため、初めにソート済みbamファイルをpileupする必要がある。output.txtオプションとして以下のようなものがある。[0]pileupしたファイルは大きいので、パイプを使い1ライナースクリプトで進めるやり方がHPに紹介されている。approaches.マニュアルリンクmacOSでも動く。brewでインストール可能。(genbank形式)、GFF3に対応。アウトプット遺伝子のアノテーションがついたgff3かgenbankファイルがあるなら、検出結果を遺伝子の部位情報付きで出してくれるので非常に便利。genbankファイルは、NCBIのftpサイトか、もしくは該当ページの情報を配列付きで表示させて全テキストコピーしてもいい。gff3ファイルはUCSCなどからダウンロードできる。pair2.fq.gz感度調整のオプションはないが、-pをつけると多型検出ができる。helpのメッセージを載せておく。MODE解析が終わるとhtmlのまとめファイルができるのが特徴。内部でRなどの描画ライブラリを動かしているぽい。上のようにアミノ酸置換まで含めて報告されている。変異部位をクリックするとアライメント結果も表示される。coveragedeletion予測時は下のようにリードのデプスを見ることもできる。indelの抜け漏れはこの中で一番少なかった。インフォマティクスのツールに不慣れなユーザーにも使いやすいと思われる。動作が複雑なためか他のソフトの安定性かわからないが、時折解析中にエラーを起こす。ゆう時は大体パラレルにBreseqを動作させていた時だったりする。複数走らせないほうよい。バクテリア解析用に設計されており、真核生物ゲノムだとかなり重たい。vcf出力はdata/output.vcfに出力される。不完全な変異は証拠が不十分になるため、最終出力に残らないことがある。不完全変異も調べるならoutput/output.gdをチェックする。追記4breseqだけで完結するように設計されている。そのため人の目で解釈しずらいVCFは出さない。他のツールと連携させる前提なら、VCFなど出さないbreseqはオススメしない。-pをつけると、クローナルなサンプルのサンプルのバリアント検出から、ミクスチャーのサンプルのバリアントコールに変更できる。研究室進化などのサンプルで、heterogeneticalなアレルのバリアント検出を行いたい時などにも使うことができる。approaches.ヒトゲノムをターゲットにしている。腫瘍サンプル特異的な変異の差分検出もできる。オプションにもdad(父)とかmam(母)とかある。人のサンプルを使うことがあれば試してみたい。一応ランの流れを書いておく。scalpel解析にはbedフォーマットのポジション情報が必要(リンク)。bcftoolsを使ってbamファイルから変換することができる(SEQanswers)。bamファイルはあらかじめソートしてindexをつけておく。Runにはデータベースファイルが必要。ref.fa上のコマンドを実行するとoutdir/フォルダに色々なファイルができる。その中にvariants.db.dirというファイルができている。これがデータベースファイル。データベースファイルができないエラーもあるらしい。ref.fa人データでテストしたところ、エラーなく最後まで走ったが、コールがゼロになる。下にリンクを追加しました。ランにはsamファイルが必要。推奨される作り方がマニュアルに書いてあるのでそれに従う。backtrackを使う。aln_pe.samテスト用のsamとfaが用意されている。下記URLからダウンロードする。https://sourceforge.net/projects/breakseek/files/toy_testset/以下はそのファイルで進行している。alnでsamを作ったら、sam_prep.pyでクロモソームごとに分割する。test.sam終わるとfasta内のクロモソームごとにディレクトリが作られ、中にsamができる。このsamファイルを使ってランする。インサートサイズのSD。正規分布ならインサートサイズの10-20%くらいらしい。dataset...が出てから、解析途中は何もlogが出てこない。初めはランに失敗したと思っていたが、強制終了せず気長に待ったところ、1時間近くかかった。3Mのバクテリアゲノムを使っただけでこれは長い...。プログラムの最適化の度合いが低いのかもしれない。alnでなくてもOK。bwaのindexアルゴリズム。"bwasw"か"is"。approaches.performance比較の論文で高評価されているツール。動作に必要なものがgitのダウンロードページに記載されている。(http://www.scipy.org/)Anaconda環境なら以下のコマンドでワンライナーインストールできる(bioconda)。scanindel----------------------------------------------------------------------------インストール当時はcondaに頼らず導入した。上の依存の大半はbrewとpipコマンドですぐ導入できたが、Inchwormはmacで動作しなかった。14.04にインストールした。brewで導入できるコマンドでないものがあれば入れておく。listでインストール済みリストのチェック。list下のコマンドのうちlistででないものがあれば導入。freebayespysam、PyVCF、Biopython、SciPy、NumPyもpythonのpipコマンドで導入できる。PyVCFはpythonのSNV検出スクリプト。freezeないものがあれば導入。assemblerはTrinitiyの中のアセンブラで知ってる人も多いと思う。TrinityのソースコードHPからダウンロードしてインストールする。上記ファイルをzipでダウンロードして解凍。中にあるInchwormディレクトリに移動。Inchwormをビルド。installbin/内に実行ファイル"inchworm"ができているはず。これにパスを通すかusr/bin内にコピー。これでようやくランの準備が整った。fileを読み込み、その内容に従って解析される。そのため依存ツールのリファレンスindexのパスを書いたconfig.txtと、サンプルのパスを書いたsample.txtを準備しておく必要がある。以下のようなファイルを用意した。ref.2bitで作り、上の内容通りgenome/に置いておく。Human?chr19.fastaで作成し、fastaと一緒にgenome/に配置しておく。Freebayesは同じfastaを指定するだけでO.K。fastqのパスを指定する。今回は解析するカレントディレクトリにそのままfastqを置いたので、上のようになる。あとは以下のコマンドを打つとScanindelが実行される。callまで自動で進む。-------参考1----------------------------------------------------------------------scanindelはリードとリファンレスの相同部位検索にblatを用いている。blatはblastに類似した相同部位検索プログラムで、indexファイルの作成がいらないなどの特徴をもつ。詳細はこちらblatはスタンドアローンとwebサーバー版(UCSC)がある。最初にbinaryファイルを作るところは同じで、以下のコマンドを打つ(上の通り)。ref.2bitそのあとサーバー版ではgfServerコマンドを使いポートを開ける。&localhostは自分で指定する名前。17774は自分で指定するポート番号。----エラーが出たとき-----ただしたくさんのスクリプトを統合しているためか、エラーが出てきたときはlogファイル(gfserver.temp.log)を開きどこで止まったか確認して進める。memのプロセスでエラーが出た。input_sorted.bamにしないとsortred.bamができないバージョンがあるかららしい。この件はgithubのScanindelのところでもコメントで表示されている。memのソートでエラーが出ていた場合、pythonスクリプトを以下のように修正。ScanIndel.pyをテキストエディタで開く。>"+output+".bwa.sorted.bam"に変更。'.temp.sorted.bam'に変更すれば動く。-------参考2----------------------------------------------------------------------ScanIndelは内部でblatを動かし、自動で上のコマンドを走らせポートを開けてリファレンスとのマッチを行う。そのため、blat検索時のポート番号の指定などのオプションを備えている。ScanIndelのランに失敗すると、gfServerが停止せずポートが開きっぱなしになる。このままだと、次回ランでエラーになる。ポート番号を--gfServer_portオプションで指定するか、ポートを閉じておく。2587これで停止できる。その他、途中で止まるとログファイル(gfserver.temp.log)も残っている。これも次回ランでエラーを起こすぽいので、エラー内容を確認したら消しておく。indelを検出しようとすると、エラーログなしに解析が止まることもたくさんあった。indelを少し加えると、ランは最後まで行くようになった。indelも正確に検出されていたので、何かしら小さなバグだと思われる。ランが途中で止まる人は、コントロールで大きなindelを加えてやり直してみてもいいかもしれない。その他、シミューレーションデータで検証時、短いindelを検出しようとすると途中で止まることがあった。indelを混ぜ込むとランは正常に終了したので、対処療法的にバグ回避はできたが原因は不明。approaches.linuxのみサポートされている。PRISM_PATH=$PWDPRISMは大きく分けて4つのステップからなる。それらのtoolはtoolkitディレクトリ中に保存されている。実際のランはrun_PRISM.shを走らせるだけでよい。/home/user/PRISM_1_1_6/sample/reference/chr100.fasamファイルとfaファイルはフルパスで指定する。ランが成功すると、カレントディレクトリにPRISM_outputができ、中にdiscordantなリードペアのクラスターがsam形式で出力される。これが構造変化の種類ごとに分けて出力されるのでかなりたくさんのsamができる(詳細はREADMEに書いてある)。肝心のindelはsplit_all.sam_ns_rmmul_cigar_sorted_svに保存されている。パスを通してからのランは下のように行う。output-Iと-Oでそれぞれ作業ディレクトリと出力ディレクトリを明示する。インサートサイズとSDを大きく間違うと途中で止まることがあるみたいので注意する。approaches.ダウンロードサイト。mantatumorとnormalを比較して、高い感受性でガン変異を検出する機能を持つ。engineのような分散コンピュータ環境にも対応しているらしいが、ゲノムサイズがそれほど大きなければローカル環境で十分高速に動く。はじめにbin/の中にあるconfigManta.pyにパスを通しておく。ランは、はじめにコンフィグレーションファイルを作り、それを元に解析する2段階方式dである。1、コンフィグディレクトリの作成。ペアエンドのbamファイルの指定。リファレンスfastaの指定。(local|sge)。ローカルなら"local"を指定。SGE環境なら128まで対応。seqの変異解析を行うことも可能になっている。oputputRNAモードでは1データのみ現在は対応。Topにリンクを追加しました。感じたことシミュレーションデータでホモの構造変化の検出テストを行うと、Fermikit、ScanIndelが最も検出率が高く、コール内容も正確だった。次点でBreseq、pindelの検出率が高かった。Pindelは特にレアな構造変化の検出に強かった。これはリードの個数を指定して感度をカバレッジに応じて変化できる点が効いていると思われる。indelの検出だと、Platypusとpindelの検出感度が優れていた。2つとも、デフォルトではサポートするリードが2つでバリアントを検出する設定のため、ディープにシーケンスしたデータだとエラーも拾いやすい。screeningではオプションをつけて感度を下げて検出する工夫も必要と思われる。補足ツールによって異なる重み付けでSVはコールされており、どうやって統合するかについても論文で議論されています。1つペーパーのリンクを貼っておきます。indelを検出しVCFで出力し、さらに変異を修正したfastaも自動出力してくれます。追記トップにもリンクを張りましたが、Parliament2は簡単に動かすことができるとても使いやすいツールです。smoove注目記事はてなブログをはじめよう!kazumaxneoさんは、はてなブログを使っています。ブログを報告するスターをつけました引用をストックしましたストック一覧を見る閉じる引用するにはまずログインしてくださいログイン閉じる引用をストックできませんでした。再度お試しください閉じる限定公開記事のため引用できません。