Question: Vcf produced after ApplyVSQR step
0
gravatar for BAGeno
19 months ago by
BAGeno180
BAGeno180 wrote:

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?

applyvsqr gatk • 792 views
ADD COMMENTlink modified 19 months ago by andrew.j.skelton736.1k • written 19 months ago by BAGeno180

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

ADD REPLYlink written 19 months ago by andrew.j.skelton736.1k

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

ADD REPLYlink written 19 months ago by Nicolas Rosewick9.3k

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
ADD REPLYlink written 18 months ago by BAGeno180

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

ADD REPLYlink written 18 months ago by Nicolas Rosewick9.3k

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

ADD REPLYlink written 18 months ago by BAGeno180

This is my tranches pdf file. Can you please tell me if its fine?

ADD REPLYlink modified 18 months ago • written 18 months ago by BAGeno180

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

ADD REPLYlink written 18 months ago by andrew.j.skelton736.1k

This data is of exome sequencing

ADD REPLYlink written 18 months ago by BAGeno180
1
gravatar for andrew.j.skelton73
19 months ago by
London
andrew.j.skelton736.1k wrote:

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.

ADD COMMENTlink written 19 months ago by andrew.j.skelton736.1k
1

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

ADD REPLYlink written 18 months ago by BAGeno180
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1873 users visited in the last hour
_