BWA, GATK and snpEff + GE (human (hs37d5) only, PE)

Introduction

このパイプラインは、Illuminaで読まれたペアエンドのデータを入力として、SNV, INDELを抽出するパイプラインです。マッピングはBWAで行い、変異の抽出はGATKによって行います。セルイノベーションプロジェクトで作成した公共データベースからSNV, INDELの頻度情報を付加することもできます。さらにsnpEffによってアミノ酸置換の生じる変異を抽出し、PolyPhen2, PROVEANによるアノテーションを付与します。




Input format	FASTQ (Illumina)
Library layout	For paired-end read only
Species	Human (hs37d5) only
Execution time	About a few days (30M paired-end reads[76bp+76bp])


Results
・SNPのリスト(VCF形式)
・サンプル間の近縁関係を示す図


GATK UnifiedGenotyperによる変異抽出

Table of contents

• BWA, GATK and snpEff + GE (human (hs37d5) only, PE) [IN: fastq (paired-end) (max 10)]
  
  
  • Introduction
  • Table of contents
  • Inputs
  • Outputs
  • How to run this pipeline
  • Result explanation
  • Comments
  • Use case
  • Related Information

Inputs

・ファイル名にはピリオド「.」、ハイフン「-」やアンダーバー「_」以外の特殊文字、日本語文字などは使用できません。
・ペアエンドの配列ファイルを自動的に認識するため、「xxx_1.fastq」と「xxx_2.fastq」などの決まった配列名のパターンを必要とします。
詳しくはこちらの定義をご覧ください。

FASTQ-formatted DNA sequence file (paired-end)

FASTQ Format Specification

Outputs

実行例1




1．SNV, INDELのリスト(VCFファイル) GATK UnifiedGenotyperによる変異抽出	2．SNV, INDELの統計情報 snpEffによる統計情報
3．個体間の類似性解析(クラスタリング) 変異情報を数値化しクラスタリング	個体間の類似性解析(主成分分析) 変異情報を数値化し主成分分析
4．Genome Explorerによる可視化 変異箇所の確認	5. フィルタリングしやすいように整形されたタブ区切りファイル SNV, INDELリスト

How to run this pipeline

MASERの一般的な操作方法のページも併せてご参照ください。(アップロードについて、パイプライン実行について）

1．入力ファイル(ここではFASTQファイル)をアップロードします。アップロード時に、「Data Label」などを入力する画面がありますが、このパイプラインを用いる場合には、「Sample Name」をきちんと設定し、「Application Type」を「WGS(たとえデータはExomeだとしても)」を選択すると良いです。
サンプル名には、アルファベット＋数字＋アンダーバー「_」で名前を付けてください。


2．本パイプラインは「Resequencing」の項目の中に配置されています。


3．上記1．で適切にSample Nameを入力していない場合は、入力するサンプル名をオプションで「Sample Name」の欄に直接設定することもできます。(データ登録時にきちんとしたサンプル名を入力していた場合には、「Sample Name」欄を「DEFAULT」のままにして頂ければ、登録時のサンプル名が使用されます。データ登録時にSample Nameが空のままで、「DEFAULT」の場合には、ファイル名がサンプル名として利用されます。)


2．領域を限定して変異を抽出したい場合には、この「-L」のオプションを設定してください。
例1：　X染色体のみにする場合、「X」
例2：　1番染色体の10,000 - 20,000の間を調べる場合、「1:10000-20000」


3．そのほかのGATKのパラメーターはここで設定してください。
また、既知とするSNV, INDELのデータベースを選択してください。デフォルトではセルイノベーションで作成したデータベースとなります。

Result explanation

1．実行結果を説明します。


2．2番のファイルを開くと、snpEffによって作成されたSNV, INDELの統計情報を見ることができます。


snpEffによる統計情報


3．アウトプットの3番のファイルは、サンプルの血縁関係などを推測するため、SNVのあり、なしを数値化し、クラスタリング、主成分分析を行った結果の図です。(ホモ置換を2、ヘテロ置換を1、置換なしを0と数値化し、群平均法によるクラスタリング。)


変異情報を数値化しクラスタリング

変異情報を数値化し主成分分析



4．アウトプットの1番のファイルは、SNV, INDELのリスト(VCFファイル)です。
　上から順に、SNV.vcf、INDEL.vcf、SNV+INDEL.vcf、SNV.important.vcf (各種擬陽性フィルターをパスし、アミノ酸orスプライシング近傍に変異ありに限定)、INDEL.important.vcf(各種擬陽性フィルターをパスし、アミノ酸orスプライシング近傍に変異ありに限定)のVCFファイルがあります。


GATK UnifiedGenotyperによる変異抽出


5．実際に遺伝様式が、単一遺伝子による優性遺伝を示すレアな遺伝病の場合のフィルタリング方法を示します。まずは、アプトプットの5番のファイル「snp_indel.filtered.snpeff.simple.tsv (アミノ酸コード領域orスプライシング近傍の変異でアミノ酸に影響を与えうる変異に限定したファイル)」をダウンロードします。
それを表計算ソフトで開きます。例としてExcel2013で解析する場合の例を示します。
Excelで開く場合に、ドラッグ＆ドロップによってエクセルで上記tsvファイルを開きますと、変異の有無のセルが意図しない変換を受けるため、必ず「ファイルを開く」から上記tsvファイルを開いてください。


6．すると「テキスト　ファイルウィザード」が立ち上がります。

上図のようになっていることを確認し「次へ」をクリックします。

区切り文字がタブだけ選択されていることを確認し、次へをクリックします。

最後に「データのプレビュー」欄を右にスクロールすると、FORMATの次の列からサンプル名が表示されているかと思います。サンプルの列を選択し、「G/標準」から「文字列」へと変更してください。

7．ファイルが開かれた後は、フィルターの機能を利用して絞込みを行っていきます。


8．以下の絞込みは一例ですが、まずクオリティの低い変異を除外します。

G列の「数値フィルター」の「指定の値以上」を選び、｢30」(99.9％正しいとされるクオリティです)を入力します。

9．次に1000ゲノムデータベース中で頻度の高い変異を除外します。

D列の「数値フィルター」の「指定の値以下」を選び、「1」(1％)を入力します。場合によっては、「0｣のみを対象とするケースも多いかと思います。
この1000ゲノムデータベース中の変異に関しては、当パイプラインと同じパラメーターで変異を抽出した1000ゲノムデータベースを用いており、通常の1000ゲノムデータベースと比べ、ツールのアーティファクトを除く際により効果があります。

10．次に罹患個体で存在する変異に絞り込みます。

N列以降で、罹患個体のフィルターを開き、「0/0 (変異なし)」、「(空白セル)　(リードのデプスが不十分で変異の有無が判定できない箇所」のチェックを外します。
ちなみに、他の項目の意味は次のとおりです。
「0/1 (ヘテロ変異)」、「1/1 (ホモ変異)」、「0/2 (ヘテロ変異ですが、同じ位置に他の変異がある場合、2つ目はこのように「2｣となります。)」

11．さらに正常個体では存在しない変異に絞り込みます。

同じくN列以降で、正常個体のフィルターを開き、「0/0 (変異なし)」以外のチェックを外します。

12．これで優性遺伝の候補変異が400個程度以下に絞り込まれたのではないかと思います。
両親と子どもをトリオで読んだ劣性遺伝のExome解析であれば、候補は10個程度になることが多いです。
これらの変異がタンパク質に影響を与える可能性が高いのかどうか、外部のWEBサイトを利用して調べることもできます(ページ下部参照)。

13．実際に絞り込んだ結果について、Genome Explorerで見てみることで、変異の最もらしさを目で見ることができます。4番の結果を開きGenome Explorerを立ち上げます。まずは「Zoom in」を押して1 base/pxまで拡大して下さい。次に「Search」タブを開き、By Locationの項目を入力すると、変異のある場所にジャンプすることができます。


変異箇所の確認

Comments

現在はパイプラインの実行結果に、PolyPhen2, PROVEANのスコアが付与されるため、下記の手順を行う必要はありませんが、参考情報として残しています。

• PROVEAN+SIFTによるアミノ酸への影響度評価の手順
1．まずは、入力となるデータを作成します。入力データは例えば、19番染色体の12,502,998塩基目がCからTに変わっていた場合、「19,12502998,C,T」となります。
そのようなデータをまとめて作るため、G列を選択し、挿入をクリックします。

たとえばこの例ですと、一番上の変異は210行目の変異であるため、G210セルに次のように入力します。
「=A210&","&B210&","&E210&","&F210」

入力が終わったら、G列の一番下まで選択し、Ctrl+Dを押して残りのセルにも入力します。

これらの範囲をコピーし、PROVEANのサイトを開きます。

2．PROVEANのサイトに移動したら、先ほどコピーした内容を貼り付けます。

下にスクロールして、送信をクリックします。


3．100件程度の変異であれば1分程度で計算が終わります。

ここで、「View result table (condensed version)」をクリックすると結果が見られます。

各変異ごとに、PROVEAN, SIFTそれぞれで変異が有害であるかどうかのスコアが表示されます。


• エラー情報
・入力ファイルのリード数が極端に少なかったり、特定の染色体上の領域に偏ったりしたデータである場合、データの補正が上手くいかずエラーとなる場合があります。

Use case

Related Information

Public


SAMtools	(Original)	(NGS Surfer's wiki)
GATK	(Original)
snpEff	(Original)	(NGS Surfer's wiki)
Picard	(Original)
BWA	(Original)	(NGS Surfer's wiki)