macでインフォマティクス

macでインフォマティクス

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

様々なインフォマティクスツールを簡単に実行できるサイバー環境 CyVerse

 

 Cyberinfrastructureは、直訳するとサイバー空間のインフラとなる。計算科学の分野では大規模な計算化学の課題に対する解決策を提供するもの、というような意味で使われている(wiki)。CyVerseはこのCyberinfrastructureを提供する、様々なインフォマティクスツールを実行可能なクラウド環境である。元々、米国植物科学コミュニティに貢献するために、iPlant Collaborativeという名前で2008年にNational Science Foundationによって作成されたのが始まりと記載されている。 2015年、iPlantは、ライフサイエンスへのより普遍的なサポートを行うため、CyVerseに改名された。 スーパーコンピューティング機能へのアクセスを民主化することで、科学者が将来のソリューションを見つけるための重要なリソースを提供している。

 

CyVerseについて

http://www.cyverse.org/about

CyVerse wiki

https://pods.iplantcollaborative.org/wiki/dashboard.action

Webinar(50分)

 

Discovery Environmentにアクセスする。

f:id:kazumaxneo:20180416212647j:plain

LAUNCHをクリック。

f:id:kazumaxneo:20180416212652j:plain

log inして立ち上げる。初回はユーザー登録(Register now)をクリックしてIDを手に入れる。

f:id:kazumaxneo:20180416212656j:plain

 

ブラウザの中にデスクトップが出てきたような感じの環境が出現する。Data(macのFInderに相当)をクリックすると、ウィンドウが出現する。初回はUploadからローカルのファイルをアップロードする。

f:id:kazumaxneo:20180416213450j:plain

ログアウントしても再度ログインすれば、ウィンドウの配置も含め全て復元される。 

 

 

Simple uploadでも2GBのファイルをuploadできる。例えばgz圧縮したfastqデータをuploadしてみる。

f:id:kazumaxneo:20180416213451j:plain

 

fastq.gzを指定してuploadをクリック。

f:id:kazumaxneo:20180416213738j:plain

 

左上のボタンAppsをクリックし、Apps選択ウィンドウを表示する。ウィンドウはクリックして移動したり拡大/縮小できる。知らなければ、もはやクラウドかローカルか分からない。gzと検索。

f:id:kazumaxneo:20180416214440j:plain

4つAppsが見つかった。並列化に対応したコマンドラインのgz解凍ツールのunpigzをクリック(圧縮ならpigz)。

 

unpigz解析ウィンドウが出てくるので、Inputをクリック=> アップロードしたgz圧縮ファイルを選択。必要ならOutputのパスも指定する。

f:id:kazumaxneo:20180416215130j:plain

 

Runをクリック。

f:id:kazumaxneo:20180416215428j:plain

 

進捗状況はAnalysisボタンから確認する。

f:id:kazumaxneo:20180416215630j:plain

終わるとcompletedになる。

 

指定した出力ディレクトリに解凍したfastqが出力されている。

f:id:kazumaxneo:20180416215733j:plain

 

ダブルクリックすると中身も表示できる。

f:id:kazumaxneo:20180416221421j:plain

 

 

de novo assemblyを行ってみる。Appsウィンドウでassemblyと検索。大量のツールがヒットした(まだ下にかなりスクロールできる)。

f:id:kazumaxneo:20180416220024j:plain

 

spadesと検索。複数バージョンが登録されている。ありがたい。

f:id:kazumaxneo:20180416220210j:plain

 

spades3.8を選択。inputにペアエンドのfastq.gzを指定した。

f:id:kazumaxneo:20180416220325j:plain

 

パラメータを指定する。コマンドラインと同じようなパラメータが指定できる。?をクリックすると簡単な説明も表示される。coverage cutoffを10とした。

f:id:kazumaxneo:20180416220424j:plain

 

アセンブリなどの計算負荷の高いジョブはHigh performance computing(HPC)として扱われる。詳細はこちらに記載されている(リンク)が、HPCに属するジョブはアメリカの複数の学術機関のクラスタを使い分けてランされているらしい。HPCの並列処理を実行したりして1人でキューをたくさん発生させると(または混雑時など)、ラン開始まで数日かかる可能性がある。要はスパコンの利用時と同じである。

f:id:kazumaxneo:20180416221107j:plainHPCジョブは開始まで時間がかかるため、エラーと勘違いしてラン前に何度も同じジョブを発生させてはいけない(当然ますます遅くなる)。

 

終わると該当フォルダにアセンブリされた配列が出力される。出力されるファイルはコマンドラインでspadesを実行したときと変わらない。 

 

 

大量のツールが登録されている。

mappingで検索。

f:id:kazumaxneo:20180416221712j:plain

indelで検索。

f:id:kazumaxneo:20180416221836j:plain

SNPで検索。

f:id:kazumaxneo:20180416221930j:plain

metagenomeで検索。

f:id:kazumaxneo:20180416221937j:plain

long readで検索。

f:id:kazumaxneo:20180416222115j:plain

genomeで検索。

f:id:kazumaxneo:20180416222155j:plain

RNAで検索。

f:id:kazumaxneo:20180416222237j:plain

bedで検索。

f:id:kazumaxneo:20180416222333j:plain

taxonで検索。

f:id:kazumaxneo:20180416222418j:plain

statisticsで検索。

f:id:kazumaxneo:20180416222516j:plain

bacteriaで検索。

f:id:kazumaxneo:20180416222628j:plain

filterで検索。

f:id:kazumaxneo:20180416222709j:plain

gwasで検索。

f:id:kazumaxneo:20180416222755j:plain

chipで検索。

f:id:kazumaxneo:20180416222822j:plain

errorで検索。

f:id:kazumaxneo:20180416223024j:plain

trimmingで検索。

f:id:kazumaxneo:20180416223112j:plain

 

分野によって数に大きな偏りはあるが、このようにNGSの主流の解析ツールなら非常に充実している。TopicからもAppsは探索できる。

f:id:kazumaxneo:20180416223532j:plain

 

joveで、CyVerseを使いtranscriptome解析を行う論文が発表されている。ビデオ学術誌なので動画があります。5分で概要をつかめます。

https://www.jove.com/video/55009/leveraging-cyverse-resources-for-de-novo-comparative-transcriptomics

 Appsを組み合わせて自動で進めるワークフローも構築できるようです。詳細はWebinarやwikiでチェックしてください。

https://pods.iplantcollaborative.org/wiki/dashboard.action

 

 

 

追記

GB以上の実験データをいきなり使うと時間がかかるので、最初はコントロールとなるようなファイルの小さいテストデータで試してみてください。

 

引用

Bringing numerous methods for expression and promoter analysis to a public cloud computing service.

Polanski K, Gao B, Mason SA, Brown P, Ott S, Denby KJ, Wild DL

Bioinformatics. 2018 Mar 1;34(5):884-886.

 

Leveraging CyVerse Resources for De Novo Comparative Transcriptomics of Underserved (Non-model) Organisms

Blake L. Joyce, Asher K. Haug-Baltzel, Jonathan P. Hulvey, Fiona McCarthy, Upendra Kumar Devisetty, Eric Lyons

J. Vis. Exp. (123), e55009, doi:10.3791/55009 (2017).

 

What is cyberinfrastructure?

 

indelコールの冗長性を調べる Vindel

 

 ゲノムDNAの変異は、一塩基多型(SNP)、挿入および欠失(indels)、逆位、大規模な複製/欠失、および転座などの構造変異を含む。最近の大規模なヒトゲノムシークエンシングプロジェクト[論文より ref.1]で示されているように、これらのタイプの変異の中で、indelsはヒト集団でSNPに次いで多い変異である。しかし、新たに配列決定されたヒトゲノムが利用可能になると、新規indelsの数はSNPよりもはるかに速いペースで増加する。例えば、2011年の調査によると、79の多様なヒトゲノムで同定された約2百万のindelsの63%以上がdbSNPのものと比較して新規である[論文より ref.2]。最近のインド女性のゲノムの配列決定および分析から、彼女のindesの約84%が独特のものであることを明らかにしている(3%のSNPが独特であるのに対し、シーケンシングされたゲノムデータベースのいずれにも記載されていない(論文執筆時点)。従って、SNPと比較して、indel変異を分類する研究はまだ初期段階であり、完全な目録を得るためには激しい努力が必要である。 Indelはまた、ショートリードマッピングアルゴリズムに大きな技術的困難および挑戦を提示している。第1世代のショートリードマッパーから改良され、さまざまなマッピングプログラムとindel検出プログラムが開発され、indel検出[ref.4-12]を可能にした。しかし、シード領域でindelが発生した場合(不一致のみが許可されている場合)、マッパーとindel検出プログラムはまだリードのマッピングに失敗する可能性があるが、マッピングのパフォーマンス全体にどのように影響するかは不明である。

 indel関連の研究が急速に発展していることから、同定されたindelsの品質評価は、下流の関連研究にとってますます重要になっている。 SNPのコールと比較して、indelコールは、PCR、シーケンシング、マッピング、およびバリアントコールで発生するエラーが発生しやすくなる。これらの誤りは、indel同定[ref.13]において高い偽陽性率をもたらす可能性がある。他方では、非常に厳格な基準と非常に複雑なindelコールのプロセスのために、実際のindelsも削除する可能性がる[ref.14]。最近、dbSNPのindelsに取り組んでいる間に、著者らは複数のindelsが同じ変異をもたらすが、異なる変異として扱われることに気づいた。論文の図1は、1つの挿入と1つの削除のindel例を示している(論文PDF)。どちらのタイプでも、2つのindelはバージョンの異なるdbSNPバージョンで異なるIDで表示されるが、異なる変異を引き起こす2配列をアラインメントするとリファレンスゲノムに同じ変更を引き起こすことを明確に示している。生物学的には、例に見られるように代替位置を持つindelが存在するかもしれないが、著者らの知る限りでは、実験的に真の生物学的信号を表すものを正確に知る方法はない。後で示すように、現行のdbSNPには、このような冗長なインジケータの無視できない数が存在する。注釈が異なるだけであれば、リファレンスゲノムへの変更結果の配列ではなく、indelsを冗長と呼ぶ。この冗長な情報は実際の生物学的シグナルを明らかにせず、下流の分析を誤解させる可能性がある。観察された冗長性は、アラインメントプログラムによって生成される等しく最適な配列アラインメント、すなわち、バリアント配列とリファレンスの配列が一緒にアライメントされる場合、アラインメントプログラムは複数の最適解を計算上区別できないため、異なるアラインメントプログラムが異なるindelとして報告してくる。これらのケースでは、dbSNPでキュレーションされたindelsの一般的な問題が冗長である可能性がある。しかし、冗長性の程度を十分に調査する作業は行われていない。さらに確認するために、隣接するSNPの距離分布をdbSNPの隣接するindelの分布と比較する。すべての染色体において、隣接するSNPの距離と比較して、隣接するindelsの距離はより高い割合を示すことが見出される。例示的な例として、論文の図2は、隣接するSNP間の距離のヒストグラムと、染色体22の隣接するindel間の距離を示す。両方のタイプの距離は、単調に減少する分布となっている。しかし、SNPとは対照的に、隣接indel距離= 1の数は、他のすべての距離から顕著である。これは、SNPと比較して、互いに隣接しているか、染色体上に非常に近いindelがさらに多く存在することを示している。この結果から、近隣のindelsの冗長性を調べるさらなる動機づけが筆者らに促された。

 この論文では、indelの冗長性をチェックする方法と戦略を開発する。 テストケースとしてdbSNP indelsを使用して、ヒトのindel冗長性の範囲を調べ、スタンドアロンのindel冗長性検証パイプラインであるVindelと対応するWebツールを開発する。 パイプラインの正確性をチェックするために統計分析が適用された。 indelsは病気や癌に関連しており、様々な目的で遺伝子マーカーとして使用されているため、重複したindelsをカタログ化し、計算上の値ではなく実際の生物学的信号を表す重複しない情報を持つアノテーションを開発することが不可欠である。 著者らのVindelシステムは、この目的に必要なツールを提供する。

f:id:kazumaxneo:20180416074838j:plain

公式ページより転載。冗長なindelの検出。

 

公式ページ

http://bioinformatics.cs.vt.edu/zhanglab/software/vindel/indelRedundant.php

HPにて、筆者らの計測でNCBI's dbSNPの590万のindelのうち60万近くが冗長だったと報告している。

 

ラン

Vindel weサーバー

f:id:kazumaxneo:20180416075822j:plain

テスト時は動作しなかった。サーバーが休止しているのかもしれない。

 

引用

Vindel: a simple pipeline for checking indel redundancy

BMC Bioinformatics. 2014; 15(1): 359.

Zhiyi Li, Xiaowei Wu, Bin He, and Liqing Zhang

 

生物学的に同等な可能性があるindelのフィルタリングを行う UPS-indel

 

 Indelは、DNA配列中の塩基の挿入または欠失を意味する。2番目に主要な変異であるindelsはゲノムおよびタンパク質の進化において重要な役割を果たす。シーケンシングエラー、リードのあいまいなアライメント、異なるツールによる同じバリアントの一貫性のない表現方法などの人為的要因のため、同じ変異が異なる場所で発生することがある(論文より ref.1-3)。例えば、位置100285630から100285650までの参照配列AGGAAAGAAAGAAAGAAAGAG、および位置100285632および100285650のこの領域にそれぞれ注釈されたdbSNP、rs147659011(GAAA / +)およびrs60376183(AAGA / +)に保存された2つのindelsを考える。これらのindel突然変異は、実際には異なる位置で生じてもよいが、それらは同じ変化した配列AGGAAAGAAAGAAAGAAAGAAAGAGを生じるので、生物学的に同等である。論文の supplementary table1は、dbSNPにおける重複挿入および欠失の別の例を示す。 dbSNP、DGV(Database of Genomic Variant)、Ensemblなどの多くのデータベースは大規模な研究の結果であるため、これらのデータベースには同様のケースが存在し、重大なデータ冗長性の問題を引き起こすことがある。実際、dbSNPに保存されたヒトindelsの約10%(ref.4)、Ensemblの18%(ref.1)は重複している。重要なデータベースにおけるindelの冗長性を解決することは、その後の遺伝学研究にとって重要である。それにもかかわらず、この問題はそれにふさわしい注意を払われていない。

 等価性を考慮し冗長性を解決する多数のアプローチが開発されている。 2つの異なるVCFファイル内の同じ位置、参照、および代替の対立遺伝子を共有する場合、「厳密な一致」アプローチは2つのindelで一致する。しかしながら、ref.3に示されているように、このアプローチは、同一ではない同等のインデルを見つけることができない。 「距離ベースアプローチ」は、同じ長さを持ち、±5 bp または±25 bpのような特定の距離内に存在する場合、2つのindelsを等価として扱う。しかし、この手法では、隣り合うindelsが等価でなく、distance cutoffよりも離れた等価なindelを見逃してしまうと、誤検出が発生する。明らかに、最適な距離カットオフの選択は、2つのタイプのエラーのトレードオフである。距離カットオフが小さいほど偽陽性率は減少するが、偽陰性率は増加する。(一部略) 

 本論文では、indelの普遍的な測位システムであるUPS-indelを提案する。これはすべてのindelがすべての同等のindelが発生する可能性のある位置の範囲で表されることを意味する。この表現はVCFファイルに追加され、等価なindel表現を含むUVCFファイルを出力する。このVCFファイルにVCFファイルを追加する利点は、(1)元のVCFファイル構造が変更されていないため、UVCFファイルは多くの下流プログラムと互換性があり、(2)UPS-indel表記により、異なるVCFファイル、(3)コーディング領域と非コーディング領域の両方に重複する同等のindelについて、indelコール出力に範囲列を持つと、下流のindel注釈システムは単一の位置ではなくその範囲を考慮する。要約すると、入力VCFファイルの単純な変更を行い、indel等価情報を含むUVCFファイルを出力する。

 

indelsを比較して生物学的に同等であるかどうか分析するために使用する。 冗長なindelsを見つたり、異なるツールが出力したindelを比較するためにも利用できるツールである。

 

インストール

本体 Github

https://github.com/shabbir005/ups-indel

git clone https://github.com/shabbir005/ups-indel.git 
cd ups-indel/
make

./ups_indel

$ ./ups_indel

USAGE: ups_indel REFERENCE_FILE.fa VCF_FILE.vcf OUTPUT_FILE_NAME -hd=true.

./ups_generate_redundant_indel_list 

$ ./ups_generate_redundant_indel_list 

USAGE: ups_generate_redundant_indel_list UPS_INDEL_VCF_FILE.UVCF OUTPUT_FILE_NAME.

 > ./ups_compare_uvcf_files

$ ./ups_compare_uvcf_files

USAGE: ups_compare_uvcf_files FIRST_UVCF_FILE.UVCF SECOND_UVCF_FILE.UVCF OUTPUT_FILE_NAME

java GenerateFilteredUVCFFileAfterRemovingRedundantIndel

$ java GenerateFilteredUVCFFileAfterRemovingRedundantIndel 

USAGE: java GenerateFilteredUVCFFileAfterRemovingRedundantIndel UVCF_FILE REDUNDANT_INDEL_LIST

 

ラン

、レポジトリ内のテストデータを分析する。

ups_indel example/ref.fa example/in.vcf out -hd=true

終わるとout.uvcfができる。 UPS-COORDINATEのカラムがVCFに追加される。

f:id:kazumaxneo:20180416002458j:plain

 

、重複するvcfのリストを出力する。1で出力したuvcfファイルを使う。

ups_generate_redundant_indel_list out.uvcf redundant_indel_list.txt

> cat redundant_indel_list.txt 

$ cat example/redundant_indel_list.txt 

[rs145072688, rs555500075]

[rs368127301, rs557514207]

[rs201888535, rs61158452]

uesaka-no-Air-2:ups-indel kazumaxneo$

 

、重複するコール数とユニークなコール数を出力する。1で出力したuvcfファイルを使う。

ups_compare_uvcf_files example/sample1.uvcf example/sample2.uvcf compare_list.txt 

 > cat compare_list.txt

$ cat compare_list.txt 

Number of Common Indels in example/sample1.uvcf and example/sample2.uvcf : 106

Number of Indels in example/sample1.uvcf but not in example/sample2.uvcf : 30

Number of Indels in example/sample2.uvcf but not in example/sample1.uvcf : 20

 

、重複をフィルタリングするには、javaのコマンドを走らせる(1、2を先におこなっておく)。

java GenerateFilteredUVCFFileAfterRemovingRedundantIndel out.uvcf redundant_indel_list.txt

出力はinputのファイル名に"filtered"をつけたout_filtered.uvcfが出力される。

 

引用

UPS-indel: a Universal Positioning System for Indels.

Hasan MS, Wu X, Watson LT, Zhang L

Sci Rep. 2017 Oct 26;7(1):14106.

 

VCFのコンセンサスコーラー CGES

 

 Whole-exome sequencing(WES)は、疾患に寄与する稀少変異を同定するための手頃なアプローチとなった。過去5年間で、PubMedのキーワード「exome sequencing」で索引付けされた論文の数は200倍に増加し、ヒトの遺伝学の明確な傾向を表している。生物学的メカニズムを明らかにするためのWESの有用性は、問題の表現型の遺伝的構造、配列決定技術の質、および配列の変異を同定および遺伝子型を決定するために用いられる分析方法にかなり依存する。近年、Atlas-SNP2(Challis et al、2012)、GATK(DePristo et al、2011; McKenna、2010a)、SeqEM(Martin et al、2010)を含むraw WESデータを解析するためのいくつかの方法が開発されている)、FreeBayes(Garrison and Marth、2012)、SAMtools(Li et al、2009a)、Dindel Albers et al (2011)、SOAPsnp(Li et al、2009b)、Varscan2(Koboldt et al、2012)などがある。

 これらの方法は、全エキソームおよび全ゲノム配列決定の両方を含む次世代配列決定(NGS)データの分析における実質的な努力および専門知識を表す。ここでは、CGESと呼ばれるコンセンサスコーリングアプローチを紹介する。このアルゴリズムは、シカゴ大学医学部とComputation Instituteの共同作業として開発され、4つの既存のアルゴリズム(GATKv2.8、Atlas-SNP2、SAMtoolsおよびFreeBayes)の多様なバリアントコール戦略を利用している。オープンソースの、自由に利用可能なユーザーフレンドリーな分析プラットフォームである。

 すべてのバリアントコールプログラムは、シーケンスデータの一部のコアプロパティ(リードのデプスやalleleの数など)と比較してパフォーマンスを最適化しようとするが、他のディメンションとは異なる場合がある。 GATKv2.8、SAMtools、Altas-SNP2、およびFreeBayesにコンセンサスベースのパイプラインを構築することを選択した。これらの4つのアルゴリズムは補完的なアプローチを使用するためである。たとえば、すべてのアルゴリズムには遺伝子型とindel尤度モデルが含まれているが、モデルそのものは使用される情報やそのような情報に与えられる重み付けが異なっている[O'Rawe et al 2013)とYu and Sun 2013を参照)。さらに、FreeBayesでは、サンプルにコピー数のバリエーションがある場合は明示的なパラメータ化が可能で、GATKv2.8は特定のデータセットで訓練できる洗練されたフィルタリングオプションを提供する。これらのプログラムは、タイプIおよびタイプIIの両方の誤差を低減するために、以前の変異観察、連鎖不均衡構造および構造変化を含む、可能な限り多くの情報を統合しようとするアルゴリズムのスイートを構成する。

 簡単に言えば、CGESは最初に各アルゴリズムは別々に実行し、その結果得られる4つの遺伝子型コールの集合を結合して、サイズは大きくなるが平均品質は低い3つの出力セットを作成する:コンセンサスコール(4つのアルゴリズムすべてによるコール)、部分的コンセンサスコール(例えば、3つ以上のアルゴリズムによって作られたもの)と、すべてのコール(1つ以上のアルゴリズムによって行われた呼び出し)との和集合である。 CGESはまた、各アルゴリズムからのクオリティスコアを調和させて、QCレポートおよびpublish品質のプロット(CGES-QCツール)を提供する。

 論文の図1に示すように、CGESとCGES-QCは多段階のパイプラインを形成し、複数のコーラーを実行する必要がある。Galaxyプラットフォーム(Goecks et al、2010)を使用して、パイプラインの各ブランチおよびCGES自体を形成している。 Galaxyは、NGSパイプラインを実行するために必要な計算上の専門知識を削減し、Amazon Web ServicesGoogleMicrosoftなどのパブリッククラウドで実行できるため、Galaxyを使用することでNGSデータ分析を民主化する利点がある。自閉症の被験者とその家族に収集されたexomeデータにCGESのコンセンサスコールアプローチを適用した。この調査の結果を使用して、CGESアプローチの能力を実証し、プロジェクトレベル、バリアントレベル、サンプルレベルおよびファミリーベースのメトリックを、すべてのアルゴリズムに提供する。

 

インストール

依存

  • pysam
  • pyvcf
pip install pysam pyvcf

本体 Github

https://github.com/vtrubets/galaxy.consensus

git clone https://github.com/vtrubets/galaxy.consensus.git
cd galaxy.consensus/

 > python consensus_tool/consensus_genotyper.py --help

$ python consensus_tool/consensus_genotyper.py --help

usage: consensus_genotyper.py [-h] [--site-threshold SITETHRESH]

                              [--genotype-threshold GENOTHRESH]

                              [--ignore-missing]

                              VCFS [VCFS ...]

 

Find sites and genotypes that aggree among an arbitrary number of VCF files.

 

positional arguments:

  VCFS                  List of VCF files for input.

 

optional arguments:

  -h, --help            show this help message and exit

  --site-threshold SITETHRESH, -s SITETHRESH

                        Number of inputs which must agree for a site to be

                        included in the output.

  --genotype-threshold GENOTHRESH, -g GENOTHRESH

                        Number of inputs which must agree for a genotype to be

                        marked as non-missing.

  --ignore-missing, -m  Flag specifying how to handle missing genotypes in the

                        vote. If present, missing genotypes are excluded from

                        the genotype concordance vote unless all genotypes are

                        missing.

 

 

ラン

テストラン用のデータとして、Githubレポジトリのdataディレクトリに、4つのVCFファイルが準備されている。Atlas-SNP2、GATK UnifiedGenotyper、FreeBayes、SAMtools(mpileup)のバリアントコール結果で、それぞれ132のサンプルを解析したvcfファイルとindexファイル (vcf.gz.tbi) である。

以下のようにそれぞれのvcf.gzを指定するか、vcfのディレクトリを指定してランする(data/*vcf.gz)。

python consensus_tool/consensus_genotyper.py --genotype-threshold 4 --site-threshold 4 data/atlas.vcf.gz data/freebayes.vcf.gz data/gatk.vcf.gz data/mpileup.vcf.gz output.vcf
  •  --site-threshold   Number of inputs which must agree for a site to be included in the output.
  • --genotype-threshold   Number of inputs which must agree for a genotype to be marked as non-missing.

使用するvcfはvcftoolsのvcf-sortなどでポジションソートされている必要がある。

 

 

: 現在の最新のレポジトリはエラーがあるかもしれない。pythonでよく出るエラーメッセージが出てくる。注意して使ってください。ためしてませんが、galaxyで使ったほうが良いかもしれません。

 

引用

Consensus Genotyper for Exome Sequencing (CGES): improving the quality of exome variant genotypes

Vassily Trubetskoy, Alex Rodriguez, Uptal Dave, Nicholas Campbell, Emily L. Crawford, Edwin H. Cook, James S. Sutcliffe, Ian Foster, Ravi Madduri, Nancy J. Cox, and Lea K. Davis

Bioinformatics. 2015 Jan 15; 31(2): 187–193.

 

 

samtools 使い方 mpileup ( calling SNPs ) & annotation

http://bioinfo-dojo.net/2016/07/13/samtools_calling_snps/

 

VCFのフィルタリングを行うGUIツール FMFilter

 

 遺伝病研究における次世代技術の使用が普及している。 exomeおよび全ゲノムシーケンシングが利用可能になると、データの解析と解釈が必要になる。遺伝病の研究に使えるVarSifter [論文より ref.1]、GEMINI [ref.2]、GeneTalk [ref.3]、CanvasDB [ref.4]、ExomeSuite [ref.5]、wKGGSeq [ref.6]などの既存のツールがあり、突然変異を引き起こす私たちの病気を見つけることができる。これらのツールのほとんどは、劣性および優性の症例に対処できるが、複合ヘテロ接合およびde novo変異を分析することができる公的に利用可能なプログラムがない。

 GeneTalkは、特に複合ヘテロ接合向けに設計されているが、フルバージョンは商用である。パブリックバージョンは複合ヘテロ接合解析をサポートしていない。 GEMINIおよびwKGGSeqは、変異が段階的になっている場合にのみ複合ヘテロ接合体モデルを支持し、そうでなければ二重ヒット戦略を用いる。さらに、wKGGSeqは、KGGSeqによって病原性であると予測される変異のみを考慮する。 ExomeSuiteはダブルヒット戦略も使用する。第3.2.1節では、ダブルヒット戦略が複合ヘテロ接合解析を正確に実行しないことを示し、誤りを増加させる例を提供する。 VarSifterを用いた複合ヘテロ接合症例の取り扱いも容易ではない。著者らは、既存のツールの複合ヘテロ接合解析の結果を比較した(一部略)。

ここでは、与えられた基準に従って候補バリアントを効果的にフィルタリングする、継承モデルベースのツールであるFMFilterを提供する。著者らの知る限り、FMFilterは、複合ヘテロ接合モデルとDE novoモデルを適切に扱うことができる最初の公的に利用可能なツールである(論文執筆時点)。これは、FMFilterと、2つの公開された複合ヘテロ接合症例[ref.7,8]における強力なフィルタリングオプションを使用して承認されている。また、プログラムには使いやすいグラフィカルユーザーインターフェイスがあり、プログラミングに関する知識は必要なくなっている。このプログラムは、効率的なメモリ処理技術を考慮して設計されており、通常のコンピュータを使用して非常に大きなバリアントファイルを処理するために効果的に使用できる。

 FMFilterの利点は、既存のツールと比較して必要とされる計算リソースの量の削減である。これにより、研究者はパーソナルコンピュータでツールを実行できる。主な制約はI / O境界だが、プログラムのメモリ使用量はごくわずかである。したがって、バリアントファイルサイズはプログラムの使用を制限しない。この機能により、FMFilterは、既存のツールのほとんどが開くことができないか、多くの読み込みと実行時間を必要とする場所で、大きなVCFファイルを簡単に処理できる。論文の表1では、HDDおよびSSDのストレージを考慮したプログラムの性能を示す。 FMFilterの使用法を解釈するために、異なる病気モデルから採取したサンプルvcfを検討した。サンプルには、パブリックデータ(pigo.vcf [ref.3]および1000ゲノムのvcf)およびin houseのケースの両方が含まれる。巨大ファイル例として、1000ゲノムプロジェクトからの15ギガバイトのVCFファイルと28ゲノムからなる100ギガバイトの社内VCFファイルを検討した。

 

 

 

公式ページ

http://fmfilter.sourceforge.net

マニュアル

http://fmfilter.sourceforge.net/manual.html

 

インストール

windows版とlinux版が用意されている。

SourceForge

https://sourceforge.net/projects/fmfilter/files/Fast%20Model%20Based%20Variant%20Filtering%20Tool-1.0-Linux-x86_64-Install/download

 

ラン

ここではwindows版を説明する。

 

上記からダウンロードしたexeファイルを実行する。

f:id:kazumaxneo:20180415141350j:plain

 支持に従いインストールする。

 

起動する。

f:id:kazumaxneo:20180415141437j:plain

 

BrowseからVCFを読み込む。UK10K プロジェクト(nature)のexampleデータなどが利用できればよかったが(UK10K Study Samples)、ダウンロードするには当然log inが必要だが、dbSNPのように登録申請が手間そうだったので、1000genomesのVCFを読み込んだ。1000genomesは家族、系統などの情報も消されているが、ここでは仮に両親と子供(患者)、unaffectedが以下のIDとする。複数選択は続けてクリックする。

f:id:kazumaxneo:20180415152001j:plain

発症していない両親(ヘテロキャリア)、発症した子供から、ホモの劣勢形質の原因遺伝子の変異を検出したいので、Recessiveを選択(両親にない新規変異ならばDe novo)。

ファミリーのトリオやカルテットの解析ではなくて影響ありと影響なしの2群比較なら、affectedとunaffectedだけ指定する。複合ヘテロ接合体(compound heterozygous)の解析ならCompound Heterozygousを選ぶ(See manual)。

 

ほかのフィルター条件追加できる。カバレッジ、Genotype qualityを指定できる。Othersは0だと、該当患者のみの出力になり、それ以外の数値だと最大で指定した数値までの患者以外のhitも許容するようになるらしい。populationsのデータ解析で役に立つとマニュアルに書かれている。

f:id:kazumaxneo:20180415163357j:plain

Annotationの項目でフィルタイングするには、ToolメニューからAnnotation Finderを起動し、既知バリアントにアノテーションをつけたVCFを作って、それを読み込む必要がある。(Tools => Annotation Finder)。

f:id:kazumaxneo:20180415144958j:plain

 

 右側のウィンドウは、GTフィールドに記載があれば表示される。たとえばGATKではこのフィールドがある(VCF)。

f:id:kazumaxneo:20180415165627j:plain

 マニュアルより転載。

 

ほかに、大規模プロジェクトのバリアントデータの頻度を利用してフィルタリングするMAFのウィンドウも用意されている(上の写真の右下)。MAF(リンク)のフィルタリングを行うには、Tools => Annotation Finderのアノテーション段階で、大規模プロジェクトのバリアント情報をVCFに取り込んでおく必要がある。

 

準備ができたら、右下のstartでランする。おわると指定したパスにフィルタリングされったVCFが出力される。

 

 論文で書かれている1000genome のVCFをchr1の10%だけ読み込んだが(およそ1.8GB)、エラーを起こした。GT yieldがないためか?(未検討)

 

引用

FMFilter: A fast model based variant filtering tool.

Akgün M, Faruk Gerdan Ö, Görmez Z, Demirci H.

J Biomed Inform. 2016 Apr;60:319-27.

 

ターゲットに特異的なコア配列のプライマーを設計する RUCS

 

 ポリメラーゼ連鎖反応(PCR)は、分子生物学における最も重要な科学的進歩の1つである。これは、DNAの特定の配列をコピーするための安価な技術である。 PCRは、医療、法医学、および研究のアプリケーションに不可欠なツールになっている。 PCRは、感染性病原体の検出および同定、ならびに病原性および耐性遺伝子のタイピングおよび特徴付けに使用される。

 PCR反応において予想されるアンプリコンのみを産生するプライマーを設計して見つけることは、面倒で、しばしば反復的なプロセスであり得る。プライマーが結合しようとするゲノムは巨大であり、完全一致でなくてもプライマーはDNAに結合する可能性がある。この非特異的プライミングは、PCRが検出目的のために使用される場合、陽性の結果につながる可能性がある。

 優れたプライマー候補の同定と偽陽性の減少を助けることのできるいくつかのバイオインフォマティクスツールがすでに存在している。これらのツールの中で、最も顕著なものはssGeneFinder、Primer3およびPrimerBLASTである(Ho et al、2012; Untergasser et al、2012; Ye et al、2012)。 ssGeneFinderは、BLASTを利用して特定のデータセットに固有のDNA配列を検索する(Camacho et al、2009)。 Primer3は、所定のDNA配列に対して適切なPCRプライマーを予測するためのツールを提供する。 PrimerBLASTは、特定のデータベースと一致するプライマーの除外も可能にすることにより、Primer3の有用性を拡大する(紹介)。これらのツールを組み合わせPCRプライマー候補を見つけることができるが、プロセスは最適ではなく、適切なプライマー候補を得るためにいくつかの手作業を必要とする。さらに、ssGeneFinderツールは非常に特異的であり、したがって、一塩基多型(SNP)または小さな挿入/欠失のような小さなユニークな配列を欠いてしまう。

 著者らが探した限り、PCR反応が非特異的プライミングの問題を引き起こすかどうかを予測し、ドラフトを直接作成する既存のツールは見つからなかった。これは、適切なペアが見いだされる前に、実験室で異なるプライマー対を用いた反復プロセスを用いることがしばしば必要であることを意味する。最も近いツールはFastPCRだったが、このツールをドラフトゲノムで機能させることはできなかった(Kalendar et al、2017)。この論文では、PCRプライマー設計の負担を軽減するための2つの新規な方法を紹介する。第1の方法は、陽性及び陰性ゲノムの所与のデータセットについて固有のコア配列を同定するための迅速かつ高感度な方法である。第2の方法は、Primer3のプライマー対予測と、PCRプライマープライマーセットについてのアンプリコンを、正および負のセットの基準に対して新規なin silico PCR産物の検証方法と組み合わせる。固有のコア配列を同定するこの方法をssGeneFinderと比較し、6.5-20倍高感度であることを見出した。さらに、mCR-1コリスチン耐性遺伝子を含むゲノムを標的とするプライマー対を設計するために、RUCSを使用した。予測された対のうちの3つを、PCRおよびゲル電気泳動を用いた実験的検証のために選択した。全ての3つの対は、mcr-1を含有する試料の標的長さを有するアンプリコンを首尾よく増幅し、そして陰性試料について増幅産物は産生されなかった。この論文で提示された新しい方法は、標的配列を同定するのに必要な時間を短縮し、迅速に仮想PCRを検証して、あいまいに結合するプライマーで無駄になる時間を排除することができる。

 

 webサーバ

https://cge.cbs.dtu.dk/services/RUCS/

Instructions

https://cge.cbs.dtu.dk/services/RUCS/instructions.php

 

ラン

instructionにしたがって解析してみる。公式ページに用意されているデータをダウンロードする。

 

 公式データは実践的な内容になっている。はじめにシーケンスデータをアセンブリするところから始める。 アセンブリしたターゲットゲノム(contig)に対してプライマーを設計していく。

アセンブリは、公式マニュアルではspadesの解析サーバーを利用している(spadesサーバー)。ここではローカル環境でアセンブルしてcontig配列を得る。

spades.py -k auto -1 ERR1399396_1.fastq.gz -2 ERR1399396_2.fastq.gz -o assembly

終わったら、scaffolds.fastaをERR1399396.fastaにリネームする。

positivesのフォルダとnegativesのフォルダを作成し、以下のようにファイルを配置する。

f:id:kazumaxneo:20180412175957j:plain

positivesのフォルダには、ターゲットに設定した全てのゲノム配列を入れる。negativesのフォルダにはネガコンのゲノムを入れる。このファイル配置にしたがって、positiveで特異的増幅が期待できるプライマーペアが探索される。

 

webサーバーに移動する。

https://cge.cbs.dtu.dk/services/RUCS-1.0/

 

テストでは、増幅サイズを100-300に変更する。

f:id:kazumaxneo:20180412180344j:plain

リファレンスはpositiveのゲノム(ERR1399396.fasta)に指定する。 

f:id:kazumaxneo:20180412180405j:plain 

 最後にファイルを指定する。

f:id:kazumaxneo:20180414091546j:plain

ボタンはなくて、uploadが終了すると自動で解析画面にジャンプする。サーバーが混雑している場合、メールアドレスを記載するよう促される。テストした時は、結果のメールが来るまで2日かかった。

また、不具合かどうか不明だが、Safariではうまくページジャンプできなかった。chromeに切り替えると自動ページジャンプした。

 

結果

みつかったペアが表示される。

f:id:kazumaxneo:20180414092703j:plain

3−5列は以下のような定義に基づいている。

、Sensitivity - fraction of positive genomes producing the amplicon

、Specificity - fraction of negative genomes NOT producing the amplicon

、Noise - number of produced amplicons significantly different to the target size*

 

+をクリックすると配列が表示される。

f:id:kazumaxneo:20180414092657j:plain

このほか、コア配列、検討されたプライマーなどのファイルがダウンロードできる。

 

引用

RUCS: rapid identification of PCR primers for unique core sequences.

Thomsen MCF, Hasman H, Westh H, Kaya H, Lund O.

Bioinformatics. 2017 Dec 15;33(24):3917-3921.

 

ノーマライズしてVCF間の比較時のバイアスを減らす BAN

 

 Variant Call Format(VCF)は、遺伝的変異および遺伝子型に関する情報を格納するためのタブ区切りのテキスト形式である(論文より Petr et al、2011)。 VCFファイル中の変異のレコードは、リファレンスDNA配列を試料DNAのシーケンスに変換する情報を記憶する。最も簡単な形式では、各バリアントレコードには、リファレンス(POS)に変更を加えなければならない位置と、 リファレンス(REF)に存在するサブシーケンスの2つのサブシーケンスが含まれる。 2つは、サンプルの配列(ALT)を作成するために他の部分配列を置き換える代替の配列である。実際、VCFファイルは、リファレンス配列とサンプル配列との間の関係のテキスト記述であり、一方、参照配列への試料配列のアライメントは、その関係の図式表現である。シーケンスはいくつかの方法で整列できるため、複数のVCFファイルで同じサンプルシーケンスをさまざまな方法で表現できる。そのようなVCFファイルは、異なる種類のバリアントを含んでいても、暗黙的に同等と見なされる。(一部略 論文の図1(リンク)を使った説明)。

 2つのVCFファイルを比較する場合、バリアントの標準化された表現が欠如していることが問題となる。 VCF比較の1つの重要なアプリケーションは、バリアントコールワークフローの評価と査定である(Liu et al、2013; O'Rawe et al、2013; Pirooznia et al、2014)。これらのワークフローは、Mapping、Alignment Filtering、Variant Calling、Variant Filtering(Van der Auwera et al、2013)などのいくつかのプロセスで構成されている。アルゴリズムまたはプロセスで使用されるパラメータに適用される変更は、最終出力VCFファイルに影響する。最終的なVCFファイルは前のすべてのプロセスの精度を反映しているため、最終的なVCFファイルと実際のセットとを比較して、どれほど類似しているかを確認する必要がある。この類似性は、ワークフローの正確性を評価する指標と考えることができる。ただし、VCFファイル内の同じサンプルシーケンスの表現が異なるため、類似性が不明確であり、評価が難しい場合がある。

 異なるVCF比較ツールを使用して2つのVCFファイルを比較すると、比較器は同じ入力ファイルに対して同様の結果を出力することが期待されるが、各VCFファイルの何百万ものバリエーションを持つ実際のデータに暗黙的な等価性(図1のものなど)が存在するため、全一致する結果はめったに得られない。コンパレータ(comparators: 比較水準器)は、2つのバリアントセット(または2つの個別のバリアント)が同等かどうかを判断するために異なる基準を検討するため、VCFファイル間のさまざまな同値数を識別する。まったく同じセットのバリアントが表されている場合、すべての比較プログラムがそれらを識別して同様に報告する。したがって、セクション3で説明したメトリックを使用してさまざまなVCF比較ツールの出力間の不一致を測定することによって、入力VCFファイルがどれだけ正常に正規化されるかを判断できる。VCFファイルが正常に正規化されれば、暗黙の等価性がより少ない場合、低い不一致が観察される。特定のVCF正規化方法で正規化された数組のVCFファイルを比較すると、比較器間の平均不一致は、その特定のVCF正規化方法の有効性を示す。この論文では、Best Alignment Normalization(BAN)と呼ばれる新しいVCF正規化手順が導入され、これまで文献で提案されている正規化手順と比較して不一致は著しく低くなっていた。

 

 

 Variant Call Format(VCF)は、変異に関するデータを保存するために広く使用されている。バリアントコールのワークフローは、ショートリードシーケンスからの潜在的な変異を検出し、それらをVCF形式で報告する。バリアントの呼び出し元の精度を評価するには、標準のバリエーションセットを含む参照VCFファイルと出力を正しく比較することが重要になる。しかしながら、VCFファイルの比較は、個々のゲノム変異体をいくつかの異なる方法で表すことができ、また異なるソフトウェアによって独自の方法で報告されるとは限らないため、複雑な作業である。
 この論文では、より正確なVCFファイルの比較をもたらすBest Alignment Normalization(BAN)と呼ばれるVCF正規化方法を紹介している。 BANは、VCFファイルのすべてのバリエーションをリファレンスゲノムに適用してサンプルゲノムを作成し、このサンプルゲノムをリファレンスゲノムにアライメントして変異を再コールする。BANの目的は、VCF比較時に正確な結果を得ることであるため、異なるVCF比較器の出力間の不一致が少なくなるように、より良い正規化方法を定義する。

 

公式ページ

https://sites.google.com/site/banadf16/

 

Start Guide Video & How to run BAN

 

インストール

ダウンロードするには、右上のダウンロードアイコンをクリックする。

https://drive.google.com/file/d/0B-5UdOXcwlPcYWgweHhDNHJQRTA/view

マニュアルPDFもダウンロードされる。 

#解凍
tar -xvf BANv1.tar

ban.shを開き、各コマンドのパスを修正する。Bioperlとmummer2vcf_m.plはダウンロードに含まれている。以下のようにした。

f:id:kazumaxneo:20180413224946j:plain

 >  bam.sh

ban.sh

 

 

======================================

Best Alignment Normalization (BAN)

======================================

Command:               

Input fasta file name: 

Input VCF file name:   

Normalization with available phasing (Input VCF must be phased)

======================================

 

 

USAGE:

Program prep/hap/dip ref.fasta input.vcf R/P

 

 

ラン 

はじめにリファレンスのindexを作成する。.dictとfasta.faiができる。picardとsamtoolsで作ってもよい。

ban.sh prep refrence.fasta

 

haploidのVCFをノーマライズする。

ban.sh hap reference.fasta input.vcf

解析が終わると、VCFにしたがって修正されたリファレンスゲノム、正規化されたVCF、正規化されたVCFにしたがって修正されたリファレンスゲノム、そして2つの修正されたリファレンスゲノムをmummerツールの1つnucmerで比較したファイルが出力される(manual PDFにも少し記載されている。) 。

オリジナルのVCFと正規化したVCFに全く違いがなければ、最後の比較ファイルの出力はゼロバイトになる。

 

 diploidのVCFをノーマライズする。

ban.sh dip reference.fasta input.vcf

 

引用

Improved VCF normalization for accurate VCF comparison.

Bayat A, Gaëta B, Ignjatovic A, Parameswaran S

Bioinformatics. 2017 Apr 1;33(7):964-970.