This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Vcf produced after ApplyVSQR step

Hi,

I am using variant recalibration for the first time. I have performed both steps of VariantRecalibrator and ApplyVSQR. I want to know whether I am performing steps correctly.

Here are my commands

gatk VariantRecalibrator -R new_hg38.fa \
                         -V allsample.vcf \
                         --resource hapmap,known=false,training=true,truth=true,prior=15.0:hapmap_3.3.hg38.vcf.gz \
                         --resource omni,known=false,training=true,truth=false,prior=12.0:omni2.5.hg38.vcf.gz \
                         --resource 1000G,known=false,training=true,truth=false,prior=10.0:thouzndG_phase1.snps.high_confidence.hg38.vcf.gz  \
                         --resource dbsnp,known=true,training=false,truth=false,prior=2.0:All_20180418.chr.hg38.vcf.gz \
                         -an QD \
                         -an MQ \
                         -an MQRankSum \
                         -an ReadPosRankSum \
                         -an FS \
                         -an SOR \
                         -mode SNP \
                         -O output.snp.recal  \
                         --tranches-file output.snp.tranches \
                         --rscript-file output.snp..plots.R

And

gatk ApplyVQSR -R new_hg38.fa \
               -V allsample.vcf \
               -O output.snps.vcf \
               --truth-sensitivity-filter-level 99.0 \
               --tranches-file output.snp.tranches  \
               --recal-file output.snp.recal \
               -mode SNP

This is my VCF.

chr1    13273   rs531730856 G   C   1023.88 VQSRTrancheSNP99.90to100.00 AC=6;AF=0.250;AN=24;AS_InbreedingCoeff=0.6320;AS_QD=12.80;BaseQRankSum=0.551;DB;DP=177;ExcessHet=0.0742;FS=12.051;InbreedingCoeff=0.6320;MLEAC=6;MLEAF=0.250;MQ=31.37;MQRankSum=0.589;NEGATIVE_TRAIN_SITE;QD=12.80;ReadPosRankSum=1.08;SOR=2.124;VQSLOD=-6.853e+00;culprit=MQRankSum    GT:AD:DP:GQ:PL  0/1:36,27:63:99:647,0,878   0/0:29,0:29:81:0,81,1215    1/1:0,4:4:12:116,12,0   ./.:0,0:0:.:0,0,0   1/1:0,4:4:12:129,12,0   0/0:8,0:8:24:0,24,246   ./.:0,0:0:.:0,0,0   0/0:19,0:19:51:0,51,765 0/1:1,8:9:0:197,0,0 0/0:10,0:10:27:0,27,405 0/0:4,0:4:12:0,12,131   0/0:5,0:5:15:0,15,204   0/0:16,0:16:39:0,39,585 ./.:0,0:0:.:0,0,0   0/0:6,0:6:15:0,15,225
chr1    13417   rs777038595 C   CGAGA   1426.83 .   AC=4;AF=0.167;AN=24;AS_InbreedingCoeff=0.2256;AS_QD=19.55;BaseQRankSum=-1.620e-01;DB;DP=154;ExcessHet=0.7463;FS=0.000;InbreedingCoeff=0.2256;MLEAC=5;MLEAF=0.208;MQ=21.09;MQRankSum=0.666;QD=19.55;ReadPosRankSum=-6.660e-01;SOR=1.179  GT:AD:DP:GQ:PL  0/1:20,22:42:99:709,0,1128  0/0:17,0:17:45:0,45,675 1/1:0,8:8:24:280,24,0   0/0:11,0:11:27:0,27,405 ./.:4,0:4:.:0,0,0   0/0:5,0:5:15:0,15,182   ./.:0,0:0:.:0,0,0   0/0:16,0:16:42:0,42,630 0/1:9,14:23:99:501,0,464    0/0:4,0:4:12:0,12,141   0/0:9,0:9:21:0,21,315   0/0:3,0:3:0:0,0,41  0/0:7,0:7:21:0,21,285   ./.:0,0:0:.:0,0,0   0/0:3,0:3:0:0,0,41
chr1    13649   rs879707275 G   C   104.2   VQSRTrancheSNP99.00to99.90  AC=2;AF=0.100;AN=20;AS_InbreedingCoeff=-0.2109;AS_QD=5.48;BaseQRankSum=0.00;DB;DP=53;ExcessHet=3.2736;FS=0.000;InbreedingCoeff=-0.2109;MLEAC=2;MLEAF=0.100;MQ=17.80;MQRankSum=0.00;QD=5.48;ReadPosRankSum=-3.350e-01;SOR=0.446;VQSLOD=-2.378e-01;culprit=MQ GT:AD:DP:GQ:PL  0/0:5,0:5:15:0,15,189   0/1:3,4:7:53:78,0,53    ./.:2,0:2:.:0,0,0   ./.:3,0:3:.:0,0,0   0/0:1,0:1:3:0,3,19  0/1:8,4:12:63:63,0,162  ./.:0,0:0:.:0,0,0   0/0:1,0:1:3:0,3,19  ./.:0,0:0:.:0,0,0   ./.:0,0:0:.:0,0,0   0/0:4,0:4:0:0,0,82  0/0:1,0:1:3:0,3,19  0/0:8,0:8:24:0,24,304   0/0:3,0:3:9:0,9,96  0/0:6,0:6:18:0,18,218
chr1    14574   rs28503599  A   G   71.21   VQSRTrancheSNP99.90to100.00 AC=5;AF=0.208;AN=24;AS_InbreedingCoeff=0.2613;AS_QD=14.24;BaseQRankSum=0.00;DB;DP=39;ExcessHet=0.0824;FS=0.000;InbreedingCoeff=0.2613;MLEAC=4;MLEAF=0.167;MQ=18.97;MQRankSum=-9.670e-01;NEGATIVE_TRAIN_SITE;QD=14.24;ReadPosRankSum=0.967;SOR=2.303;VQSLOD=-4.449e+00;culprit=MQ    GT:AD:DP:GQ:PL  0/0:13,0:13:39:0,39,447 ./.:2,0:2:.:0,0,0   ./.:0,0:0:.:0,0,0   0/0:1,0:1:3:0,3,19  0/0:1,0:1:3:0,3,19  0/1:1,2:3:21:40,0,21    1/1:0,2:2:6:54,6,0  0/0:3,0:3:9:0,9,107 0/0:1,0:1:3:0,3,19  ./.:0,0:0:.:0,0,0   0/0:4,0:4:12:0,12,140   0/0:2,0:2:6:0,6,59  0/0:3,0:3:0:0,0,22  1/1:0,1:1:3:26,3,0  0/0:3,0:3:9:0,9,100
chr1    14590   rs707679    G   A   17.26   VQSRTrancheSNP99.90to100.00 AC=2;AF=0.077;AN=26;AS_InbreedingCoeff=0.2729;AS_QD=8.63;DB;DP=47;ExcessHet=0.1047;FS=0.000;InbreedingCoeff=0.2729;MLEAC=1;MLEAF=0.038;MQ=22.51;NEGATIVE_TRAIN_SITE;QD=8.63;SOR=2.303;VQSLOD=-3.776e+00;culprit=MQ  GT:AD:DP:GQ:PL  0/0:18,0:18:51:0,51,765 0/0:3,0:3:0:0,0,41  0/0:2,0:2:6:0,6,78  0/0:1,0:1:3:0,3,19  ./.:0,0:0:.:0,0,0   0/0:3,0:3:0:0,0,4   1/1:0,2:2:6:54,6,0  0/0:5,0:5:15:0,15,137   0/0:2,0:2:6:0,6,59  0/0:1,0:1:3:0,3,37  0/0:3,0:3:9:0,9,100 0/0:2,0:2:6:0,6,59  0/0:3,0:3:9:0,9,100 ./.:0,0:0:.:0,0,0   0/0:2,0:2:6:0,6,59

Can you please tell me if that should be the output of variantreclibration step as I do not have FILTER field annotated as mentioned on gatk website. Also how should I proceed next whether should I filter these variants?

gatk applyvsqr
• 4,623 views • link
updated by andrew.j.skelton73 andrew.j.skelton73 • written by BAGeno BAGeno

I've formatted your post so that your code blocks are more readable

• 0 views • link
written by andrew.j.skelton73 andrew.j.skelton73

Could you post your tranche plot ( VariantRecalibrator should output you a pdf containing it). Also is it WGS ? or targeted (exome, gene panel) ?

• 0 views • link
written by Nicolas Rosewick Nicolas Rosewick

This is tranches file.

# Variant quality score tranches file
# Version number 5
targetTruthSensitivity,numKnown,numNovel,knownTiTv,novelTiTv,minVQSLod,filterName,model,accessibleTruthSites,callsAtTruthSites,truthSensitivity
90.00,2302437,371395,2.1178,1.4051,2.9424,VQSRTrancheSNP0.00to90.00,SNP,1291939,1162745,0.9000
99.00,2806818,409886,2.0969,1.4087,0.9599,VQSRTrancheSNP90.00to99.00,SNP,1291939,1279019,0.9900
99.90,2907437,425455,2.0915,1.4023,-2.7966,VQSRTrancheSNP99.00to99.90,SNP,1291939,1290647,0.9990
100.00,2982702,454343,2.0834,1.3928,-36326.3108,VQSRTrancheSNP99.90to100.00,SNP,1291939,1291939,1.0000
• 0 views • link
written by BAGeno BAGeno

could you post the plot. It's easier to read ;)

• 0 views • link
written by Nicolas Rosewick Nicolas Rosewick

Could you please tell me how should I plot graphs? I have R script made by VariantRecalibrator and also downloaded ggplot2 in R.

• 0 views • link
written by BAGeno BAGeno

Can you clarify if this is WES / WGS / Targeted sequencing?

• 0 views • link
written by andrew.j.skelton73 andrew.j.skelton73

This data is of exome sequencing

• 0 views • link
written by BAGeno BAGeno

1 answer

I'd check out this post which tries to explain VQSR all in one. Particularly look at the How ApplyRecalibration works in a nutshell section:

During the first part of the recalibration process, variants in your callset were given a score called VQSLOD. At the same time, variants in your training sets were also ranked by VQSLOD. When you specify a tranche sensitivity threshold with ApplyRecalibration, expressed as a percentage (e.g. 99.9%), what happens is that the program looks at what is the VQSLOD value above which 99.9% of the variants in the training callset are included. It then takes that value of VQSLOD and uses it as a threshold to filter your variants. Variants that are above the threshold pass the filter, so the FILTER field will contain PASS. Variants that are below the threshold will be filtered out; they will be written to the output file, but in the FILTER field they will have the name of the tranche they belonged to. So VQSRTrancheSNP99.90to100.00 means that the variant was in the range of VQSLODs corresponding to the remaining 0.1% of the training set, which are basically considered false positives.

• 0 views • link
written by andrew.j.skelton73 andrew.j.skelton73

So should I filter any sites which are not in my defined range of tranches?

• 0 views • link
written by BAGeno BAGeno

Log in to answer this question.

More posts like this
  • Expanding gVCF box from GATK
    written by Axzd • Axzd
    Hello, I wanted to expand the gVCF block of my vcf but I got this as a result ##bcftools_convertCommand=convert --gvcf2vcf -f Adineta_vaga.fsa; Date=Tue Aug 1 …
  • why CombineGCVFs in gatk not showing all the samples name?
    written by rj.rezwan • rj.rezwan
    Hi, I combined the 64 `.vcf` files using the `CombineGVCFs` in `gatk`. The command was completed successfully but its showing the output with only one …
  • HaplotypeCaller GENOTYPE GIVEN ALLELES doesn't genotype given alleles
    written by kamanovae • kamanovae
    Hi! I am trying to run the Gatk HaplotypeCaller (human data): ``` ./gatk-4.1.2.0/gatk HaplotypeCaller\ --reference ref.fa \ --input file.bam \ --genotyping-mode GENOTYPE_GIVEN_ALLELES \ --alleles allele_chunk_file.vcf …
  • Impute genotype probabilities with BEAGLE
    written by dario.galanti92 • dario.galanti92
    I have a vcf file of a plant species obtained with GATK4 and I would like to impute few missing genotypes, including genotype probabilities or …
  • How to remove non chromosome (1-22) from a .VCF file
    written by jmukisa90 • jmukisa90
    Hi there, I am performing the GATK pipeline for converting fastq files to VCF format. How do I remove, non-chromosome 1-22 variants from the column …
  • GATK GenotypeGVCFs changes HET to REF_ALT
    written by FGV FGV
    Dear all, I've been using GATK HaplotypeCaller / GenotypGVFs (v4.2.4.0) for a while but, recently found something strange. There is a position (7063) with 8 …
  • Strange Allele Frequency in vcf file
    written by GP_PG • GP_PG
    Hi everyone. I performed a snp calling using GATK and following all the suggestions in the BestPractices sections. I performed single snp calling for two …
  • why is there "." in Filter?
    written by arshiya.g02 • arshiya.g02
    This was the command I used `sudo ./gatk HaplotypeCaller -R /home/arshiya/efs/FastQC/hg38.fa -I /home/arshiya/efs/FastQC/sample_sort_reads.bam -O /home/arshiya/efs/FastQC/sample_haplo.vcf.gz` and the following was the output I got chr1 7820343 …
  • Diploid calls on chrY?
    written by dec986 dec986
    Hello, I recently got my genome sequenced. chrY 59030922 . A G 1244.77 PASS AC=1;AF=0.500;AN=2;BaseQRankSum=-6.580e-01;ClippingRankSum=0.00;DP=41;ExcessHet=3.0103;FS=9.230;MLEAC=1;MLEAF=0.500;MQ=37.33;MQRankSum=-4.563e+00;QD=30.36;ReadPosRankSum=-7.570e-01;SOR=0.076;VQSLOD=0.452;culprit=FS GT:AD:DP:GQ:PL 0/1:9,32:41:99:1273,0,180 chrY 59032665 . A G 418.77 PASS AC=1;AF=0.500;AN=2;BaseQRankSum=-4.182e+00;ClippingRankSum=0.00;DP=52;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=45.59;MQRankSum=-6.082e+00;QD=8.21;ReadPosRankSum=-1.528e+00;SOR=0.756;VQSLOD=2.78;culprit=DP …
  • Large VCF QUAL Values
    written by pld pld
    I'm currently working on calling variants for a virus that was sequenced. For lack of a better solution, I'm using the GATK UnifiedGenotyper pipeline. I'm …