macでインフォマティクス

macでインフォマティクス

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

small indelとSNV検出のワークフロー

SNVやsmall indel検出については精度の高いワークフローがすでに確立されている。例えば下記のニューヨーク大のHP

https://gencore.bio.nyu.edu/variant-calling-pipeline/

には、SNVとsmall indel検出ワークフローが記載されている。流れを説明すると

 

bam作成と変異のコール(raw call)

PCR duplicateの除去

リアライメント(indel周辺のアライメント修正)

SNVとindelに仕分けしてfalse callぽい部位のフィルタリング

raw callを元にしたbase quality補正 (BQSR)

変異の再コール

SNVとindelに仕分けして再度フィルタリング

SNVのアノテーション (indelは不可)

 

となる。ワークフロー前半では、一度アライメントしてからPCRバイアスを除き、indel周辺のアライメント修正、エラーの可能性が高い部位のフィルタリング、塩基のクオリティ補正 (BQSR)を行っている。変異の検出前にアライメント精度をできる限り高めることで、偽陽性になりそうな部位をあらかじめ排除しているわけである。同様の手順はindel検出時のアライメント方法による感度差を比較した論文

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4619817/pdf/BMRI2015-456479.pdf

中でも取られており(論文のFigure 1参照)、妥当な方法と言えそうである。

解説付きで上記のワークフローを書いておく。

  

 準備編はこちら(必要なツールなど)

 

--検出の流れ--

ここではリファンレンスをref.fa、ペアリードのForwardとReverseをR1.fq、R2.fqとしている。

 

 

<追記>dictファイルとfasta.faiファイルがない人は最初に作っておきます。

0、リファレンスのindexとdictファイルの作成

samtools faidx ref.fa #ref.fa.faiができる。
picard CreateSequenceDictionary R=ref.fa O=ref.dict #ref.dictができる

 

1、リファレンスへのリードのmapping。

bwa mem -t 8 -M -R '@RG	id:sample_1	LB:sample_1	PL:ILLUMINA	PM:HISEQ	SM:sample_1' ref.fa R1.fq R2.fq > aligned_reads.sam

-t: スレッド数 (デファルト1)

-R: リードのヘッダー情報。この後必要。

-M: mark shorter split hits as secondary(このあと必要)

マッピング情報を記載したsamファイルが出力される。

 

2、bam変換とソート(リードの名前順に並べ替える)

java -jar picard.jar SortSam INPUT=aligned_reads.sam OUTPUT=sorted_reads.bam SORT_ORDER=coordinate

sorted_reads.bamが出力される。samtoolsでソートした時との違いは、アメリエフさんのブログで丁寧に解説されています。

 

3、アライメント情報の分析

java -jar picard.jar CollectAlignmentSummaryMetrics R=ref.fa I=sorted_reads.bam O=alignment_metrics.txt

CollectAlignmentSummaryMetricsコマンドはマッピングのクオリティとノイズフィルターをパスしたリードの割合などが記された基準ファイル(詳細は

Tool documentationを参照)。イルミナリードのみ対応しているらしい。alignment_metrics.txtが出力される。

 

4、3でソートしたbamを使いインサートサイズを分析

java -jar picard.jar CollectInsertSizeMetrics INPUT=sorted_reads.bam OUTPUT=insert_metrics.txt HISTOGRAM_FILE=insert_size_histogram.pdf

このコマンドでインサートサイズを分析したテキストファイル (insert_metrics.txt) とpdfファイル (insert_size_histogram.pdf) が出力される。

 

5、リファレンス全部位のリードのデプスを分析。

samtools depth -a sorted_reads.bam > depth_out.txt

-a: 全部位を出力

samtools depthの出力はタブ区切りのreference name, position, coverage depthからなるファイル。depth_out.txtが出力される。

 

6、duplicateリードの分析

java -jar picard.jar MarkDuplicates INPUT=sorted_reads.bam OUTPUT=dedup_reads.bam METRICS_FILE=metrics.txt

dedup_reads.bamが出力される。duplicateリードについてはPicardのページに説明があるが、ペアリードの5'側の配列を指標に探すらしい。duplicate候補リードは、base call qualityの合計値からプライマリー/duplicateかの判定をつけられ、出力されるbamではタグがつけられる。

 

7、bamのインデックスファイルの作成

java -jar picard.jar BuildBamIndex INPUT=dedup_reads.bam

duplicateチェック済みのbamファイルにインデックスをつける。インデックスファイル (bam.bai) はbam参照に使われるファイル。

 

8、リアライメントのターゲットを見つける(第1ステップ)

java -jar GenomeAnalysisTK.jar -T RealignerTargetCreator -R ref.fa -I dedup_reads.bam -o realignment_targets.list

リアライメントの重要性はCLCのサポートなどに丁寧に説明されている。 次世代データのアライメント手法は非常に高速だが精度の面でイマイチのため、indel周辺でミスアライメントが起きてしまう。それを丁寧にアライメントして修復する。この第一ステップではリアライメントのリストを作る。

* 古いイルミナデータの場合、クオリティ定義が今の定義と異なるためエラーになることがある。その場合は、--fix_misencoded_quality_scoresオプションをつけてステップ8、9をランする。他のデータにつけてはいけない。

 

9、リアライメントする(第2ステップ)

java -jar GenomeAnalysisTK.jar -T IndelRealigner -R ref -I dedup_reads.bam -targetIntervals realignment_targets.list -o realigned_reads.bam

第一ステップで見つけたリストを元にリアライメントを実行する。補正されたrealigned_reads.bamが出力される(実行前後のbamをIGVに読み込んでindel周辺を確認してみるとよい)。

 

10、バリアントのコール(1回目) 

java -jar GenomeAnalysisTK.jar -T HaplotypeCaller -R ref.fa -I realigned_reads.bam -o raw_variants.vcf

リアライメントされたbamを使って 一回目の変異検出を実行する。変異を記載したVCF (variant call format) ファイルが出力される。

ここで検出されたバリアントを元に以下の作業が行われる。

 

11、SNVとindelの仕分け。#SNV

java -jar GenomeAnalysisTK.jar -T SelectVariants -R ref.fa -V raw_variants.vcf -selectType SNP -o raw_snps.vcf

文字通り塩基置換とindelを分ける。

-selectType SNP: 塩基置換を抽出。

 

12、SNVとindelの仕分け。#indel

java -jar GenomeAnalysisTK.jar -T SelectVariants -R ref.fa -V raw_variants.vcf -selectType SNP -o raw_snps.vcf

-selectType indel: indelを抽出。

 11と12の作業は後半でリアライメントした後にもう一度行われる(以降の修正により結果が変わる可能性があるため)。

 

13、 SNVのフィルタリング(1回目) 

java -jar GenomeAnalysisTK.jar -T VariantFiltration -R ref.fa -V raw_snps.vcf --filterExpression 'QD < 2.0 || FS > 60.0 || MQ < 40.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0 || SOR > 4.0' --filterName "basic_snp_filter" -o filtered_snps.vcf

VariantFiltrationの詳細は公式ページを参照。--filterExpressionで条件を指定する。上記のフィルタリング条件は

QD < 2.0 #QualByDepth。QUAL(マッピングクオリティに近い)をカバレッジで正規化

FS > 60.0 #FisherStrand。フィッシャー検定してリードの方向のbiasを評価したもの。biasなければ0近い。

MQ < 40.0 #RMSMappingQuality。その部位の全リードのマッピングクオリティ二乗平均平方根

MQRankSum < -12.5 #MappingQualityRankSumTest マッピングクオリティを ウィルコンの Rank Sum Test で評価したもの。

ReadPosRankSum < -8.0

となっている。コマンド中の2重のパイプ" || "はorを意味する。4つのパラメータについてはこちらを参照。上記の数値にした理由が書かれているが、結論としてはかなり寛大な条件となっている。

 

14、 indelのフィルタリング(1回目) 

java -jar GenomeAnalysisTK.jar -T VariantFiltration -R ref.fa -V raw_indels.vcf --filterExpression 'QD < 2.0 || FS > 200.0 || ReadPosRankSum < -20.0 || SOR > 10.0' --filterName "basic_indel_filter" -o filtered_indels.vcf

フィルタリング条件は

QD < 2.0

FS > 200.0

ReadPosRankSum < -20.0

SOR > 10.0

4つのパラメータについてはこちらを参照。ここでは、SORはindelのフィルタリングにのみ使っている。

 

15、ベースクオリティスコアの再較正(第1ステップ)

java -jar GenomeAnalysisTK.jar -T BaseRecalibrator -R ref.fa -I realigned_reads.bam -knownSites filtered_snps.vcf -knownSites filtered_indels.vcf -o recal_data.table

塩基のクオリティスコアはシーケンサーブラックボックスの1つで、算出の正確な方法は公開されていないが、それなりに信頼できる値ではある。だが各部位をおそらく独立に評価しており、複数塩基のパターンから信頼度を評価すると、例えばポリAの部位などのように、2つ目以降の塩基のクオリティが高すぎる部位があるとされる。そこで開発チームが機械学習や統計手法で算出したデータベースをを元に、各々のデータで検出された変異(ポジコン)とそれ以外の部位のミスマッチ(ネガコンとみなす)を比較し、全部位のクオリティスコアを再較正する。これをGATKのチームはBase Quality Score Recalibration (BQSR)と呼んでいる。(詳細は論文とHP参照)。

BQSR補正はfalse positiveを減らしtrue callを増やすために必要なstepと考えられている

 (indel検出のワークフローを比較した論文 (図1) でも実施されている)。

 

16、 ベースクオリティスコアの再較正(第2ステップ)

java -jar GenomeAnalysisTK.jar -T BaseRecalibrator -R ref -I realigned_reads.bam -knownSites filtered_snps.vcf -knownSites filtered_indels.vcf -BQSR recal_data.table -o post_recal_data.table

 

17、 ベースクオリティスコアの再較正レポート出力

java -jar GenomeAnalysisTK.jar -T AnalyzeCovariates -R ref.fa -before recal_data.table -after post_recal_data.table -plots recalibration_plots.pdf

 

18、ベースクオリティスコア分析に基づく修正

java -jar GenomeAnalysisTK.jar -T PrintReads -R ref.fa -I realigned_reads.bam -BQSR recal_data.table -o recal_reads.bam

ここでBQSR補正されたbamが出力される。

 

 19、バリアントのコール(2回目) 

java -jar GenomeAnalysisTK.jar -T HaplotypeCaller -R ref.fa -I recal_reads.bam -o raw_variants_recal.vcf

補正したbamを使って改めて変異のコールを行う。VCFファイルが出力される。 

 

  20、SNVとindelの仕分け。#SNV

java -jar GenomeAnalysisTK.jar -T SelectVariants -R ref.fa -V raw_variants_recal.vcf -selectType SNP -o raw_snps_recal.vcf

改めて仕分けを行う。

 

21、SNVとindelの仕分け。#indel

java -jar GenomeAnalysisTK.jar -T SelectVariants -R ref.fa -V raw_variants_recal.vcf -selectType INDEL -o raw_indels_recal.vcf

 

22、 SNVのフィルタリング(2回目) 

java -jar GenomeAnalysisTK.jar -T VariantFiltration -R ref.fa -V raw_snps_recal.vcf --filterExpression 'QD < 2.0 || FS > 60.0 || MQ < 40.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0 || SOR > 4.0' --filterName "basic_snp_filter" -o filtered_snps_final.vcf

フィルタリング条件は

QD < 2.0

FS > 60.0

MQ < 40.0

MQRankSum < -12.5

ReadPosRankSum < -8.0

SOR > 4.0 

となっている。1回目になかったSOR > 4.0が追加されている(SQRについてはFisherStrand (FS)の説明ページを参照)。他はのパラメータは1回目と同じ条件になっている。

 

23、 indelのフィルタリング(2回目) 

java -jar GenomeAnalysisTK.jar -T VariantFiltration -R ref.fa -V raw_indels_recal.vcf --filterExpression 'QD < 2.0 || FS > 200.0 || ReadPosRankSum < -20.0 || SOR > 10.0' --filterName "basic_indel_filter" -o filtered_indels_recal.vcf

フィルタリング条件は

QD < 2.0

FS > 200.0

ReadPosRankSum < -20.0

SOR > 10.0

 

24、変異部位のアノテーション

java -jar snpEff.jar -v snpeff_db filtered_snps_final.vcf > filtered_snps_final.ann.vcf

snpeff_dbが自分の使ったリファレンスのデータベース。データベース名が不明だったりリファレンスとデータベースが合わない場合、下記のsnpEffのデータベースをダウンロードする例を参照してください。

vcfをbedに変換すれば、bedtoolsのannotateを使って似たことができると思います(bedtools説明)。

 

 

補足1

このワークフローはlarge indelを検出することはできないが、large indel検出ツールの多くはbamファイルをインプットに使う。そのためこのワークフローで修正して出力したbamファイルはlarge indel検出にも利用することができる。

lage indel検出ツールについては以下のエントリーにまとめています。

 

補足2

変異の検出後、複数サンプルデータの共通変異をまとめたり固有変異を抽出するなら、以下のエントリーを参考にしてください。

 

補足3

複数の同じサンプルファイルがある時、例えば同じサンプルをマルチプレックスや違うレーンで複数読んだりしている場合、どこでデータを結合するかということがしばしば話題になる。

データをmergeするステップは無数に存在するが、gatkのチームは、独立してBQSR補正を行い、gatkのコール時に複数サンプルを入れる方法を1つのやり方として提案している。理由の1つに、リアライメント補正は1サンプルあたりのリードが増えると重くなることをあげている。

http://gatkforums.broadinstitute.org/gatk/discussion/6057/at-what-point-should-i-merge-read-group-bam-files-belonging-to-the-same-sample-into-a-single-file

ただし、カバレッジが非常に薄いデータがある場合、リアライメントができない恐れがあり、そういったデータの場合は前もってマージして置いた方が良いらしい。

 

補足4

 コマンドラインからIGVを起動してダイレクトにVCFファイルを読み込むには、以下のようにします。gtfやgenebankのファイルもあるなら、一緒に読み込めます。

igv -g input.fa filtered_snps_final.vcf,recal_reads.bam,input.gtf

igvをmacに入れていない場合、

brew install igv

で入ります。

IGVについては以下の記事も参考にしてください。


補足5

indelは1つのファイルで出力されます。数を数えたりフィルタリングするならBEDOPSのコマンドを使えます。

vcf2bed --snvs < filtered_snps_final.vcf|wc -l        #SNVをカウント
vcf2bed --insertions < filtered_indels_recal.vcf|wc -l  #Insertionをカウント
vcf2bed --deletions < filtered_indels_recal.vcf|wc -l  #Deletionをカウント

 vcf2bedは以下で紹介しています。


補足6

small indel変異をシミュレートするならいくつかツールはありますが、 ここではwgsimを紹介します。シーケンスbiasは再現できませんが、シーケンスエラー、SNV、indelを一定の確率で発生できるツールで、変異があるfastqを人工的に作ることができます。

 

 

 

今回の手法は一般的なショートリードでのwhole genome DNA sequenceに適応できると思われるものです。ただし、データによりこのワークフローが常にベストと一概に言えるわけではないことに注意してください。データによってはもっと高感度な手法も考えられます。

 

 

 

--解析例--

最初のコールと比べてどのくらいの違いが生じるのか調べてみた。ある程度ゲノムサイズが大きくないとわかりにくいので、クラミドモナスのリシーケンスデータをSRA検索サイトから探してダウンロードした。    

 

データ

- クラミドモナス (75bp x 2)

・ERR022944_1.fastq 2600万リード(6.9GB)

ERR022944_2.fastq 2600万リード( 6.9GB)

v3.1.29.dna.genome 120 Mb (1.2億) 

 

結果

 - クラミドモナス

1回目のコール後(上記のstep10)

SNV 111768

indel 23933

 

1回目のフィルタリングでPASSできなかった部位(上記のstep13、14)

SNV 10804

indel    258

 ↓

BQSR補正後(上記のstep20、21) フィルタリングされた部位も含む

SNV 113745

indel 24673

 ↓ 

2回目のフィルタリングでPASSできなかった部位(上記のstep22、23)

SNV 11163

indel 276

 ↓ 

Haplotypecallerの最終コール (PASSしている部位のみ)

SNV 102582

indel 24397

 

 

 

他のツールで検出した部位と比較する。bamはBQSR補正後のものを使い、フィルターも同様の条件で設定した。

FreeBayesで検出

SNV 118652 

UnifiedGenotyperで検出

SNV 139751

 

3ツールで共通の部位(bcftools isecで計算)

58816部位

 

ベン図は書かないが、UnifiedGenotyperとはhaplotypecallerは共通部分が多く、Freebayesはhaplotypercallerとはかなり食い違った。数が増えた減っただけでは直接の参考にはならないが、変化に対して大雑把な感覚を掴んでください。

 

 

 

 

 

1-24の作業をシェルスクリプトにしてまとめておけば、アノテーションまでほぼ自動で進めることができる。補足4のようなigvの起動まで組み込んでおくとちょっとかっこいいかもしれない。

今回は単体サンプルを想定したが、大規模データ、例えば人のexomeやゲノムの大規模プロジェクトなどでは複数サンプルを同時に解析したりSNPデータベース情報を使いVQSR補正を行いエラーを排除してtrue callを絞っていくやり方もある。そのあたりはまたブログに書いていく。

 

 

gatkのワークフローについてもっと知りたければ、gatkのホワイトペーパーを参照してください。

リンクをせておきます。

https://www.intel.com/content/dam/www/public/us/en/documents/white-papers/deploying-gatk-best-practices-paper.pdf

https://www.intel.com/content/dam/www/public/us/en/documents/white-papers/deploying-gatk-best-practices-paper.pdf