GATK4 最佳實踐-生殖細胞突變的檢測與識別

歡迎關注"生信修煉手冊"！

GATK4 對於體細胞突變和生殖細胞突變的檢測分別給出了對應的pipeline:

Germline SNPs+Indels

Somatic SNVs + Indels

本篇主要關注生殖細胞突變的分析流程Germline SNPs+Indels。示意圖如下：

圖中紅色方框部分的從Analysis-Ready Bam 到，主要包括以下4步

HaplotyperCaller in GVCF mode

ImportGenomicsDB Consolidate GVCFs

GenotypeGVCFs

Filter Variants by Variabt Recalibration

官網給出了6套用於參考的pipeline

主要分析步驟都差不多，這裡我選擇第4個通用的流程 ，網址如下

https://github.com/gatk-workflows/gatk4-germline-snps-indels

1.  HaplotyperCaller in GVCF mode

對於每個樣本，採用HaplotyperCaller計算突變位點，命令如下

gatk --java-options "-Xmx6G -XX:GCTimeLimit=50 -XX:GCHeapFreeLimit=10" \    HaplotypeCaller \    -R ${ref_fasta} \    -I ${input_bam} \    -L ${interval_list} \    -O ${output_filename} \    -contamination 0 -ERC GVCF

ref_fasta代表參考基因組的fasta文件;input_bam代表預處理階段產生的 bam文件；interval代表interval list文件，如果指定這個參數，只會輸出指定區域的突變信息。對於全基因組測序，不需要這個參數，對於外顯子和目的區域捕獲測序, 則需要這個參數；output_filename代表每個樣本輸出的gvcf文件的名字。

2. ImportGenomicsDB Consolidate GVCFs

將所有樣本的gvcf文件合併，產生一個總的gvcf文件，命令如下：

gatk --java-options -Xmx2G  \    MergeVcfs \    --INPUT ${sep=' --INPUT ' input_vcfs} \    --OUTPUT ${output_filename}

3. GenotypeGVCFs

包括兩個步驟，第一步，導入MergeVcfs合併好的gvcf文件， 命令如下

gatk --java-options "-Xmx4g -Xms4g" \    GenomicsDBImport \    --genomicsdb-workspace-path ${workspace_dir_name} \    --batch-size ${batch_size} \    -L ${interval} \    --sample-name-map ${sample_name_map} \    --reader-threads 5 \    -ip 500

workspace_dir_name代表輸出目錄;batch_size指定同時訪問的最大文件數，默認值為0，表示同時訪問所有樣本的文件；interval代表interval list文件，如果指定這個參數，只會輸出指定區域的突變信息。對於全基因組測序，不需要這個參數，對於外顯子和目的區域捕獲測序, 則需要這個參數；sampple_name_map是一個文件，記錄了樣本名稱和每個樣本對應的gvcf文件的路徑。格式如下

sample1      sample1.vcf.gz
sample2      sample2.vcf.gz
sample3      sample3.vcf.gz

第二步， 運行GenotypeGVCFs，命令如下

gatk --java-options "-Xmx5g -Xms5g" \    GenotypeGVCFs \    -R ${ref_fasta} \    -O ${output_vcf_filename} \    -D ${dbsnp_vcf} \    -G StandardAnnotation \    --only-output-calls-starting-in-intervals \    --use-new-qual-calculator \    -V gendb://$WORKSPACE \    -L ${interval}

需要注意-V 參數，指定的是GenomicsDBImport輸出目錄的路徑。

4. Filter Variants by Variabt Recalibration

第一步，過濾vcf文件，條件為ExcessHet大於給定的閾值，命令如下：

gatk --java-options "-Xmx3g -Xms3g" \    VariantFiltration \    --filter-expression "ExcessHet > ${excess_het_threshold}" \    --filter-name ExcessHet \    -O ${variant_filtered_vcf_filename} \    -V ${vcf}

excess_het_threshold指定ExcessHet的閾值；variant_filtered_vcf_filename代表輸出的vcf文件的名字；vcf代表GenotypeGVCFs 生成的vcf文件的名字。注意，不滿足條件的記錄也會出現在最終生成的vcf文件中, 只不過對應的Filter欄位的信息不是PASS。

第二步，刪除vcf文件中的基因型信息，命令如下

gatk --java-options "-Xmx3g -Xms3g" \    MakeSitesOnlyVcf \    --INPUT ${variant_filtered_vcf_filename} \    --OUTPUT ${sites_only_vcf_filename}

第三步，合併不同區間的vcf文件，並建立索引

gatk --java-options "-Xmx6g -Xms6g" \  GatherVcfsCloud \  --ignore-safety-checks \  --gather-type BLOCK \  --input ${inputs.list} \  --output ${output_vcf_name}gatk --java-options "-Xmx6g -Xms6g" \   IndexFeatureFile \   --feature-file ${output_vcf_name}

output_vcf_name代表輸出的vcf文件；inputs.list指定不同區間的vcf文件的路徑，格式如下

cohortA_chr1.vcf.gz
cohortA_chr2.vcf.gz

第四步，分別評估SNP和INDEL突變位點的質量, 命令如下

gatk --java-options "-Xmx24g -Xms24g" \    VariantRecalibrator \    -V ${sites_only_variant_filtered_vcf} \    -O ${recalibration_filename} \    --tranches-file ${tranches_filename} \    --trust-all-polymorphic \    -tranche ${sep=' -tranche ' recalibration_tranche_values} \    -an ${sep=' -an ' recalibration_annotation_values} \    -mode INDEL \    --max-gaussians 4 \    -resource mills,known=false,training=true,truth=true,prior=12:${mills_resource_vcf} \    -resource axiomPoly,known=false,training=true,truth=false,prior=10:${axiomPoly_resource_vcf} \    -resource dbsnp,known=true,training=false,truth=false,prior=2:${dbsnp_resource_vcf}gatk --java-options "-Xmx100g -Xms100g" \      VariantRecalibrator \      -V ${sites_only_variant_filtered_vcf} \      -O ${recalibration_filename} \      --tranches-file ${tranches_filename} \      --trust-all-polymorphic \      -tranche ${sep=' -tranche ' recalibration_tranche_values} \      -an ${sep=' -an ' recalibration_annotation_values} \      -mode SNP \      --sample-every-Nth-variant ${downsampleFactor} \      --output-model ${model_report_filename} \      --max-gaussians 6 \      -resource hapmap,known=false,training=true,truth=true,prior=15:${hapmap_resource_vcf} \      -resource omni,known=false,training=true,truth=true,prior=12:${omni_resource_vcf} \      -resource 1000G,known=false,training=true,truth=false,prior=10:${one_thousand_genomes_resource_vcf} \      -resource dbsnp,known=true,training=false,truth=false,prior=7:${dbsnp_resource_vcf}

resource指定建模時參考的vcf文件，可以看到對於indel和snp, 參考的資料庫不一樣。這一步只是建模，會輸出一個recalibration table文件，用於ApplyVQSR命令。

第五步，運行VQSR, 命令如下

gatk --java-options "-Xmx5g -Xms5g" \    ApplyVQSR \    -O tmp.indel.recalibrated.vcf \    -V ${input_vcf} \    --recal-file ${indels_recalibration} \    --tranches-file ${indels_tranches} \    --truth-sensitivity-filter-level ${indel_filter_level} \    --create-output-variant-index true \    -mode INDELgatk_path --java-options "-Xmx5g -Xms5g" \    ApplyVQSR \    -O ${recalibrated_vcf_filename} \    -V tmp.indel.recalibrated.vcf \    --recal-file ${snps_recalibration} \    --tranches-file ${snps_tranches} \    --truth-sensitivity-filter-level ${snp_filter_level} \    --create-output-variant-index true \    -mode SNP

input_vcf文件為GatherVcfsCloud生成的vcf文件，先校正INDEL位點，後校正SNP位點。

掃描關注微信號，更多精彩內容等著你！