随時更新
2019 11/24 condaインストール追記（v3.8）
2020 5/13 インストール追記
2021 4/24 GATKアーカイブのダウンロードリンク追加（archive）
2021 5/9 GATK for Microbesリンク追加
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参照)、妥当な方法と言えそうである。追記*2
解説付きで上記のワークフローを書いておく。
追記
BQSRについて: 使用している生物種のデータセットも実行する意味があるのか、Documentを読んでください。
https://gatkforums.broadinstitute.org/gatk/discussion/44/base-quality-score-recalibration-bqsr
レビューペーパー要約。ヒトゲノムがターゲットですが、NGS登場から比較的早くに確立されたバリアントコールのフローを知るには良いペーパーです。特にNGS解析の経験がなく、右も左もわからないという方は、王道（鉄板）のフローを知るためにも目を通してください。
2020 5/15 GATK4
勉強するにはこちらのペーパーも読んでください。
From FastQ data to high confidence variant calls: the Genome Analysis Toolkit best practices pipeline.
https://www.ncbi.nlm.nih.gov/pubmed/25431634
GATKのオフィシャルスライド
BEST PRACTICES FOR VARIANT CALLING WITH THE GATK
https://www.broadinstitute.org/partnerships/education/broade/best-practices-variant-calling-gatk-1
準備編はこちら（必要なツールなど）
condaで導入するならバージョン管理のために仮想環境にいれた方がよい。
#step1 bioconda
conda create -n GATK3.8 -c bioconda -y gatk samtools bwa
conda activate GATK3.8
#2020 5/13 picardも追加
conda install -c bioconda -y picard
#step2 ライセンスのactivate
#別途こちらからダウンロードしたGenomeAnalysisTK.jarのパスを指定する。例えば~/document/GATK3.8/GenomeAnalysisTK.jarにあるなら
gatk3-register ~/document/GATK3.8/GenomeAnalysisTK.jar
> ヘルプを表示してインストールチェック
GenomeAnalysisTK -T HaplotypeCaller -h
$ GenomeAnalysisTK -T HaplotypeCaller -h
------------------------------------------------------------------------------------
The Genome Analysis Toolkit (GATK) v3.8-1-0-gf15c1c3ef, Compiled 2018/02/19 05:43:50
Copyright (c) 2010-2016 The Broad Institute
For support and documentation go to https://software.broadinstitute.org/gatk
[Sun Nov 24 21:36:59 JST 2019] Executing on Mac OS X 10.14.6 x86_64
OpenJDK 64-Bit Server VM 1.8.0_192-b01
------------------------------------------------------------------------------------
------------------------------------------------------------------------------------
usage: java -jar GenomeAnalysisTK.jar -T <analysis_type> [-args <arg_file>] [-I <input_file>] [--showFullBamList] [-rbs
<read_buffer_size>] [-rf <read_filter>] [-drf <disable_read_filter>] [-L <intervals>] [-XL <excludeIntervals>] [-isr
<interval_set_rule>] [-im <interval_merging>] [-ip <interval_padding>] [-R <reference_sequence>] [-ndrs] [-maxRuntime
<maxRuntime>] [-maxRuntimeUnits <maxRuntimeUnits>] [-dt <downsampling_type>] [-dfrac <downsample_to_fraction>] [-dcov
<downsample_to_coverage>] [-baq <baq>] [-baqGOP <baqGapOpenPenalty>] [-fixNDN] [-fixMisencodedQuals]
[-allowPotentiallyMisencodedQuals] [-OQ] [-DBQ <defaultBaseQualities>] [-PF <performanceLog>] [-BQSR <BQSR>] [-qq
<quantize_quals>] [-SQQ <static_quantized_quals>] [-DIQ] [-EOQ] [-preserveQ <preserve_qscores_less_than>]
[-globalQScorePrior <globalQScorePrior>] [-secondsBetweenProgressUpdates <secondsBetweenProgressUpdates>] [-S
<validation_strictness>] [-rpr] [-kpr] [-sample_rename_mapping_file <sample_rename_mapping_file>] [-U <unsafe>]
[-jdk_deflater] [-jdk_inflater] [-disable_auto_index_creation_and_locking_when_reading_rods] [-no_cmdline_in_header]
[-sites_only] [-writeFullFormat] [-compress <bam_compression>] [-simplifyBAM] [--disable_bam_indexing] [--generate_md5]
[-nt <num_threads>] [-nct <num_cpu_threads_per_data_thread>] [-mte] [-rgbl <read_group_black_list>] [-ped <pedigree>]
[-pedString <pedigreeString>] [-pedValidationType <pedigreeValidationType>] [-variant_index_type <variant_index_type>]
[-variant_index_parameter <variant_index_parameter>] [-ref_win_stop <reference_window_stop>] [-l <logging_level>] [-log
<log_to_file>] [-h] [-version] [-mmq <min_mapping_quality_score>] [-o <out>] [-D <dbsnp>] [-dontTrimActiveRegions]
[-comp <comp>] [-A <annotation>] [-XA <excludeAnnotation>] [-G <group>] [-debug] [-useFilteredReadsForAnnotations] [-ERC
<emitRefConfidence>] [-bamout <bamOutput>] [-bamWriterType <bamWriterType>] [-edr] [-disableOptimizations] [-nda]
[-newQual] [-hets <heterozygosity>] [-indelHeterozygosity <indel_heterozygosity>] [-heterozygosityStandardDeviation
<heterozygosity_stdev>] [-stand_call_conf <standard_min_confidence_threshold_for_calling>] [-maxAltAlleles
<max_alternate_alleles>] [-maxGT <max_genotype_count>] [-maxNumPLValues <max_num_PL_values>] [-inputPrior <input_prior>]
[-ploidy <sample_ploidy>] [-gt_mode <genotyping_mode>] [-alleles <alleles>] [-contamination
<contamination_fraction_to_filter>] [-contaminationFile <contamination_fraction_per_sample_file>] [-out_mode
<output_mode>] [-allSitePLs] [-gcpHMM <gcpHMM>] [-globalMAPQ <phredScaledGlobalReadMismappingRate>] [-sn <sample_name>]
[-kmerSize <kmerSize>] [-dontIncreaseKmerSizesForCycles] [-allowNonUniqueKmersInRef] [-numPruningSamples
<numPruningSamples>] [-minDanglingBranchLength <minDanglingBranchLength>] [-consensus] [-maxNumHaplotypesInPopulation
<maxNumHaplotypesInPopulation>] [-minPruning <minPruning>] [-graph <graphOutput>] [-GQB <GVCFGQBands>] [-ERCIS
<indelSizeToEliminateInRefModel>] [-mbq <min_base_quality_score>] [-allelesTrigger] [-doNotRunPhysicalPhasing]
[-dontUseSoftClippedBases] [-pcrModel <pcr_indel_model>] [-maxReadsInRegionPerSample <maxReadsInRegionPerSample>]
[-minReadsPerAlignStart <minReadsPerAlignmentStart>] [-APO <activityProfileOut>] [-ARO <activeRegionOut>] [-AR
<activeRegionIn>] [-activeRegionExtension <activeRegionExtension>] [-forceActive] [-activeRegionMaxSize
<activeRegionMaxSize>] [-bandPassSigma <bandPassSigma>] [-maxReadsInMemoryPerSample <maxReadsInMemoryPerSample>]
[-maxTotalReadsInMemory <maxTotalReadsInMemory>] [-ActProbThresh <activeProbabilityThreshold>] [-filterRNC] [-filterMBQ]
[-filterNoBases]
O.K !
--検出の流れ--
ここではリファンレンスを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ファイルが出力される。
２、bam変換とソート（リードの名前順に並べ替える）
java -jar picard.jar SortSam INPUT=aligned_reads.sam OUTPUT=sorted_reads.bam SORT_ORDER=coordinate
sorted_reads.bamが出力される。samtoolsでソートした時との違いは、アメリエフさんのブログで丁寧に解説されています。
３、アライメント情報の分析
java -jar picard.jar CollectAlignmentSummaryMetrics R=ref.fa I=sorted_reads.bam O=alignment_metrics.txt
CollectAlignmentSummaryMetricsコマンドはマッピングのクオリティとノイズフィルターをパスしたリードの割合などが記された基準ファイル（詳細は
Tool documentationを参照）。イルミナリードのみ対応しているらしい。alignment_metrics.txtが出力される。
４、３でソートした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) が出力される。
５、リファレンス全部位のリードのデプスを分析。
samtools depth -a sorted_reads.bam > depth_out.txt
-a: 全部位を出力
samtools depthの出力はタブ区切りのreference name, position, coverage depthからなるファイル。depth_out.txtが出力される。
６、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ではタグがつけられる。
７、bamのインデックスファイルの作成
java -jar picard.jar BuildBamIndex INPUT=dedup_reads.bam
duplicateチェック済みのbamファイルにインデックスをつける。インデックスファイル (bam.bai) はbam参照に使われるファイル。
８、リアライメントのターゲットを見つける（第１ステップ）
java -jar GenomeAnalysisTK.jar -T RealignerTargetCreator -R ref.fa -I dedup_reads.bam -o realignment_targets.list
リアライメントの重要性はCLCのサポートなどに丁寧に説明されている。 次世代データのアライメント手法は非常に高速だが精度の面でイマイチのため、indel周辺でミスアライメントが起きてしまう。それを丁寧にアライメントして修復する。この第一ステップではリアライメントのリストを作る。
* 古いイルミナデータの場合、クオリティ定義が今の定義と異なるためエラーになることがある。その場合は、--fix_misencoded_quality_scoresオプションをつけてステップ８、９をランする。他のデータにつけてはいけない。
９、リアライメントする（第２ステップ）
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周辺を確認してみるとよい）。
１０、バリアントのコール（１回目）
java -jar GenomeAnalysisTK.jar -T HaplotypeCaller -R ref.fa -I realigned_reads.bam -o raw_variants.vcf
リアライメントされたbamを使って 一回目の変異検出を実行する。変異を記載したVCF (variant call format) ファイルが出力される。
ここで検出されたバリアントを元に以下の作業が行われる。
１１、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: 塩基置換を抽出。
１２、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の作業は後半でリアライメントした後にもう一度行われる（以降の修正により結果が変わる可能性があるため）。
１３、 SNVのフィルタリング（１回目）
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なければ０近い。
・MQ < 40.0 #RMSMappingQuality。その部位の全リードのマッピングクオリティ二乗平均平方根。
・MQRankSum < -12.5 #MappingQualityRankSumTest マッピングクオリティを ウィルコンの Rank Sum Test で評価したもの。
・ReadPosRankSum < -8.0
となっている。コマンド中の２重のパイプ" || "はorを意味する。４つのパラメータについてはこちらを参照。上記の数値にした理由が書かれているが、結論としてはかなり寛大な条件となっている。
１４、 indelのフィルタリング（１回目）
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
#2021 5/24 gatk4では微妙にオプション名が変わっている
gatk VariantFiltration \
-R ref.fa \
-V raw_indels.vcf \
-o filtered_indels.vcf \
--filter-name "basic_indel_filter"
--filter-expression xxx
フィルタリング条件は
QD < 2.0
FS > 200.0
ReadPosRankSum < -20.0
SOR > 10.0
４つのパラメータについてはこちらを参照。ここでは、SORはindelのフィルタリングにのみ使っている。
１５、ベースクオリティスコアの再較正（第１ステップ）
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
塩基のクオリティスコアはシーケンサーのブラックボックスの１つで、算出の正確な方法は公開されていないが、それなりに信頼できる値ではある。だが各部位をおそらく独立に評価しており、複数塩基のパターンから信頼度を評価すると、例えばポリAの部位などのように、２つ目以降の塩基のクオリティが高すぎる部位があるとされる。そこで開発チームが機械学習や統計手法で算出したデータベースをを元に、各々のデータで検出された変異（ポジコン）とそれ以外の部位のミスマッチ（ネガコンとみなす）を比較し、全部位のクオリティスコアを再較正する。これをGATKのチームはBase Quality Score Recalibration (BQSR)と呼んでいる。（詳細は論文とHP参照）。
BQSR補正はfalse positiveを減らしtrue callを増やすために必要なstepと考えられている
(indel検出のワークフローを比較した論文 (図1) でも実施されている)。
１６、 ベースクオリティスコアの再較正（第２ステップ）
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
１７、 ベースクオリティスコアの再較正レポート出力
java -jar GenomeAnalysisTK.jar -T AnalyzeCovariates -R ref.fa -before recal_data.table -after post_recal_data.table -plots recalibration_plots.pdf
１８、ベースクオリティスコア分析に基づく修正
java -jar GenomeAnalysisTK.jar -T PrintReads -R ref.fa -I realigned_reads.bam -BQSR recal_data.table -o recal_reads.bam
ここでBQSR補正されたbamが出力される。
１９、バリアントのコール（２回目）
java -jar GenomeAnalysisTK.jar -T HaplotypeCaller -R ref.fa -I recal_reads.bam -o raw_variants_recal.vcf
補正したbamを使って改めて変異のコールを行う。VCFファイルが出力される。
２０、SNVとindelの仕分け。#SNV
java -jar GenomeAnalysisTK.jar -T SelectVariants -R ref.fa -V raw_variants_recal.vcf -selectType SNP -o raw_snps_recal.vcf
改めて仕分けを行う。
２１、SNVとindelの仕分け。#indel
java -jar GenomeAnalysisTK.jar -T SelectVariants -R ref.fa -V raw_variants_recal.vcf -selectType INDEL -o raw_indels_recal.vcf
２２、 SNVのフィルタリング（２回目）
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
となっている。１回目になかったSOR > 4.0が追加されている（SQRについてはFisherStrand (FS)の説明ページを参照）。他はのパラメータは１回目と同じ条件になっている。
２３、 indelのフィルタリング（２回目）
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
２４、変異部位のアノテーション。
java -jar snpEff.jar -v snpeff_db filtered_snps_final.vcf > filtered_snps_final.ann.vcf
snpeff_dbが自分の使ったリファレンスのデータベース。データベース名が不明だったりリファレンスとデータベースが合わない場合、下記のsnpEffのデータベースをダウンロードする例を参照してください。
vcfをbedに変換すれば、bedtoolsのannotateを使って似たことができると思います（bedtools説明）。
補足１
このワークフローはlarge indelを検出することはできないが、large indel検出ツールの多くはbamファイルをインプットに使う。そのためこのワークフローで修正して出力したbamファイルはlarge indel検出にも利用することができる。
lage indel検出ツールについては以下のエントリーにまとめています。
補足２
変異の検出後、複数サンプルデータの共通変異をまとめたり固有変異を抽出するなら、以下のエントリーを参考にしてください。
補足３
複数の同じサンプルファイルがある時、例えば同じサンプルをマルチプレックスや違うレーンで複数読んだりしている場合、どこでデータを結合するかということがしばしば話題になる。
データをmergeするステップは無数に存在するが、gatkのチームは、独立してBQSR補正を行い、gatkのコール時に複数サンプルを入れる方法を１つのやり方として提案している。理由の１つに、リアライメント補正は1サンプルあたりのリードが増えると重くなることをあげている。
ただし、カバレッジが非常に薄いデータがある場合、リアライメントができない恐れがあり、そういったデータの場合は前もってマージして置いた方が良いらしい。
補足4
コマンドラインからIGVを起動してダイレクトにVCFファイルを読み込むには、以下のようにします。gtfやgenbankのファイルもあるなら、一緒に読み込めます。
igv -g input.fa filtered_snps_final.vcf,recal_reads.bam,input.gtfigvをmacに入れていない場合、
brew install igvで入ります。
IGVについては以下の記事も参考にしてください。
補足5
indelは１つのファイルで出力されます。数を数えたりフィルタリングするなら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は以下で紹介しています。
補足６
small indel変異をシミュレートするならいくつかツールはありますが、 ここではwgsimを紹介します。シーケンスbiasは再現できませんが、シーケンスエラー、SNV、indelを一定の確率で発生できるツールで、変異があるfastqを人工的に作ることができます。
補足7
再較正は、ゲノムがそこまで似ていない場合、時間がかかる割に改善はほとんどないかもしれません。
Evaluation of variant identification methods for whole genome sequencing data in dairy cattle.
https://www.ncbi.nlm.nih.gov/pubmed/25361890
Impact of post-alignment processing in variant discovery from whole exome data.
https://www.ncbi.nlm.nih.gov/pubmed/27716037
補足8
結果の可視化。htmlで結果をまとめられます。ただし万以上のコールがあるとかなり重くなります。
捕捉9
本家GATK best practice（ターゲットに合わせ、複数種類がある）は、経験者であれば、Google cloudでインスタントに使用が可能です（リンク）。試していませんが、5ドルでランできるようなパイプラインも用意されています（リンク）。
追記10
elprepを使えば上記の煩雑な処理のかなりの部分を代替できます。elrepに切り替えることも検討して下さい。
追記11
2021 5/9
Introducing #GATK for Microbes -- read all about it in the new blog post by @bhanu_gandham https://t.co/honuayNIey
— GATK Dev Team (@gatk_dev) 2021年4月27日
今回の手法は一般的なショートリードでのwhole genome DNA sequenceに適応できると思われるものです（ヒトなどの2倍体ゲノム）。ただし、データによりこのワークフローが常にベストと一概に言えるわけではないことに注意してください。データによってはもっと高感度な手法も考えられます。
1-24の作業をシェルスクリプトにしてまとめておけば、アノテーションまでほぼ自動で進めることができる。
今回は単体サンプルを想定したが、大規模データ、例えば人のexomeやゲノムの大規模プロジェクトなどでは複数サンプルを同時に解析したりSNPデータベース情報を使いVQSR補正を行いエラーを排除してtrue callを絞っていくやり方もある。そのあたりはまたブログに書いていく。
gatkのワークフローについてもっと知りたければ、gatkのホワイトペーパーを参照してください。
リンクをせておきます。
追記1
hard filteringのパラメータチューニング。
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5374681/
追記*2
GATK4では、Unifiedgenotyperは消えている。HaplotypeCallerやMutect2のバリアント検出コマンドは残っているが、この２つではhaplotype のローカルデノボアセンブリを行う関係で、リアライメントのstepは不要になっており、それに伴い、indel realignmentのコマンドも無くなっている。
https://software.broadinstitute.org/gatk/blog?id=7847
GATK4のベストプラクティスを使ったフローはまた紹介します。
2021 4/27
Introducing #GATK for Microbes -- read all about it in the new blog post by @bhanu_gandham https://t.co/honuayNIey
— GATK Dev Team (@gatk_dev) 2021年4月27日
