Test-drive the GATK tools and Best Practices pipelines on Terra


Check out this blog post to learn how you can get started with GATK and try out the pipelines in preconfigured workspaces (with a user-friendly interface!) without having to install anything.

Picard vcf filtering

jahneldarjahneldar NorwayMember

Hi, I use Picard to filter a vcf file. It seems not to function propely. I am filtering on read depth using MIN_DP=50 but I still see a lot of genotypes in the resulting file with lower coverage than 50.

Here is my call:
java -jar /tools/picard-2.9.0/picard.jar FilterVcf I=output.vcf O=output2.vcf MIN_AB=0.4 MIN_DP=50 MIN_GQ=30

Issue · Github
by shlee

Issue Number
1921
State
closed
Last Updated
Closed By
sooheelee

Best Answers

  • jahneldarjahneldar Norway
    Accepted Answer

    Ahh, so another step required?
    I'm locked in vcftools.
    Thank you for getting back to me.

    jahn

  • shleeshlee Cambridge ✭✭✭✭✭
    Accepted Answer

    Depending on your downstream analysis, physically removing the filtered sites may not be necessary. That is, downstream tools will understand that these sites are filtered. Just annotating the FILTER field is the typical approach so that you don't lose data along the way. E.g. although sites are filtered, you may wish to look into borderline cases and so would want to keep them in your file.

Answers

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    Hi @jahneldar,

    Can you give me an example of a VCF record that you expect to be filtered but is not? Thanks.

    FYI I would love to visit Norway.

  • EADGEADG KielMember ✭✭✭
    edited April 2017

    Hi @jahneldar,

    can you provide a snippet of your vcf-file ? It will be easier to reproduce the behavior if you can provide one. In addition do you also try to filter your variants via gatk- FilterVaraints oder SelectVariants ? Did these tools behave the same ?

    Greetings EADG

    @shlee Sry I was just writing my post and dont see that you answered right before me :neutral:

  • jahneldarjahneldar NorwayMember
    edited April 2017

    Yes, I did try to use the gatk, but I could not figure out how to do it. That's why I swapped to Picard. Given better 'manual/tutorial' on the FilterVariants I guess I'd stuck with that.

    BTW: I have used vcftools previously, but I got some error messages that I could not resolve when I tried to get the vcftool-filtered files into the BaseCalibrator. FYI: I am trying to follow the 'Calling Variants in RNAseq' and I am in the 'Base Quality Recalibrartion Stage'.

    Anyway, here is a snippet of my vcf-file:

    ##contig=<ID=scaffold01528,length=11040>
    ##contig=<ID=scaffold01529,length=10270>
    ##contig=<ID=scaffold01530,length=10174>
    ##contig=<ID=scaffold01531,length=6587>
    ##reference=file:///data/kvijed/IPN_Berry/AnitraRaw/scaffolds.longreads.1.fa
    #CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  100
    scaffold00001   8021    .   A   G   58.74   .   AC=2;AF=1.00;AN=2;DP=2;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=29.37;SOR=2.303 GT:AD:DP:GQ:PL  1/1:0,2:2:6:86,6,0
    scaffold00001   18375   .   T   A   29.75   .   AC=2;AF=1.00;AN=2;DP=3;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=14.87;SOR=2.303 GT:AD:DP:GQ:PL  1/1:0,2:2:6:57,6,0
    scaffold00001   20004   .   GT  G   62.73   .   AC=1;AF=0.500;AN=2;BaseQRankSum=1.718;ClippingRankSum=0.000;DP=79;ExcessHet=3.0103;FS=8.451;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=8.96;ReadPosRankSum=1.036;SOR=2.788 GT:AD:DP:GQ:PL  0/1:4,3:7:99:100,0,122
    scaffold00001   30166   .   ATAT    A   1014.73 .   AC=1;AF=0.500;AN=2;BaseQRankSum=-1.309;ClippingRankSum=0.000;DP=226;ExcessHet=3.0103;FS=18.013;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=12.08;ReadPosRankSum=3.314;SOR=2.581 GT:AD:DP:GQ:PL  0/1:55,29:84:99:1052,0,2222
    scaffold00001   31251   .   G   GA  125.73  .   
    

    Comments/help is greatly appreciated.

    jahn

    Post edited by shlee on
  • jahneldarjahneldar NorwayMember
    edited April 2017

    Hi shlee: here is the top of my input vcf file

    ##contig=<ID=scaffold01531,length=6587>
    ##reference=file:///data/kvijed/IPN_Berry/AnitraRaw/scaffolds.longreads.1.fa
    #CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  400
    scaffold00001   9252    .   A   G   99.28   .   AC=2;AF=1.00;AN=2;DP=3;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=33.09;SOR=2.833 GT:AD:DP:GQ:PL  1/1:0,3:3:9:127,9,0
    scaffold00001   11391   .   AT  A   82.73   .   AC=1;AF=0.500;AN=2;BaseQRankSum=-0.524;ClippingRankSum=0.000;DP=6;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=13.79;ReadPosRankSum=-0.431;SOR=0.693   GT:AD:DP:GQ:PL  0/1:2,4:6:49:120,0,49
    scaffold00001   12369   .   T   C   686.77  .   AC=2;AF=1.00;AN=2;DP=19;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=29.25;SOR=1.292    GT:AD:DP:GQ:PL  1/1:0,19:19:57:715,57,0
    scaffold00001   12426   .   A   C   228.80  .   AC=2;AF=1.00;AN=2;DP=7;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=32.69;SOR=1.609 GT:AD:DP:GQ:PL  1/1:0,7:7:21:257,21,0
    scaffold00001   13279   .   G   C   27200.77    .   AC=1;AF=0.500;AN=2;BaseQRankSum=0.447;ClippingRankSum=0.000;DP=1721;ExcessHet=3.0103;FS=1.106;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=15.81;ReadPosRankSum=0.975;SOR=0.781  GT:AD:DP:GQ:PL  0/1:855,866:1721:99:27229,0,26794
    scaffold00001   14031   .   C   T   52001.77    .   AC=2;AF=1.00;AN=2;BaseQRankSum=1.916;ClippingRankSum=0.000;DP=1368;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;MQRankSum=0.000;QD=27.92;ReadPosRankSum=1.733;SOR=0.825    GT:AD:DP:GQ:PL  1/1:1,1366:1367:99:52030,4082,0
    scaffold00001   14289   .   T   C   22131.77    .   AC=1;AF=0.500;AN=2;BaseQRankSum=-0.784;ClippingRankSum=0.000;DP=1465;ExcessHet=3.0103;FS=1.101;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=15.11;ReadPosRankSum=1.319;SOR=0.580 GT:AD:DP:GQ:PL  0/1:745,720:1465:99:22160,0,23245
    scaffold00001   16150   .   C   T   607.77  .   AC=2;AF=1.00;AN=2;DP=16;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=33.51;SOR=0.941    GT:AD:DP:GQ:PL  1/1:0,16:16:48:636,48,0
    

    Here is the top of my output file:

    ##contig=<ID=scaffold01530,length=10174>
    ##contig=<ID=scaffold01531,length=6587>
    ##reference=file:///data/kvijed/IPN_Berry/AnitraRaw/scaffolds.longreads.1.fa
    #CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  400
    scaffold00001   9252    .   A   G   99.28   AllGtsFiltered  AC=2;AF=1.00;AN=2;DP=3;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=33.09;SOR=2.833 GT:AD:DP:FT:GQ:PL   1/1:0,3:3:LowDP;LowGQ:9:127,9,0
    scaffold00001   11391   .   AT  A   82.73   AllGtsFiltered;AlleleBalance    AC=1;AF=0.500;AN=2;BaseQRankSum=-0.524;ClippingRankSum=0.000;DP=6;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=13.79;ReadPosRankSum=-0.431;SOR=0.693   GT:AD:DP:FT:GQ:PL   0/1:2,4:6:LowDP:49:120,0,49
    scaffold00001   12369   .   T   C   686.77  AllGtsFiltered  AC=2;AF=1.00;AN=2;DP=19;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=29.25;SOR=1.292    GT:AD:DP:FT:GQ:PL   1/1:0,19:19:LowDP:57:715,57,0
    scaffold00001   12426   .   A   C   228.80  AllGtsFiltered  AC=2;AF=1.00;AN=2;DP=7;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=32.69;SOR=1.609 GT:AD:DP:FT:GQ:PL   1/1:0,7:7:LowDP;LowGQ:21:257,21,0
    scaffold00001   13279   .   G   C   27200.77    PASS    AC=1;AF=0.500;AN=2;BaseQRankSum=0.447;ClippingRankSum=0.000;DP=1721;ExcessHet=3.0103;FS=1.106;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=15.81;ReadPosRankSum=0.975;SOR=0.781  GT:AD:DP:GQ:PL  0/1:855,866:1721:99:27229,0,26794
    scaffold00001   14031   .   C   T   52001.77    PASS    AC=2;AF=1.00;AN=2;BaseQRankSum=1.916;ClippingRankSum=0.000;DP=1368;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;MQRankSum=0.000;QD=27.92;ReadPosRankSum=1.733;SOR=0.825    GT:AD:DP:GQ:PL  1/1:1,1366:1367:99:52030,4082,0
    scaffold00001   14289   .   T   C   22131.77    PASS    AC=1;AF=0.500;AN=2;BaseQRankSum=-0.784;ClippingRankSum=0.000;DP=1465;ExcessHet=3.0103;FS=1.101;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=15.11;ReadPosRankSum=1.319;SOR=0.580 GT:AD:DP:GQ:PL  0/1:745,720:1465:99:22160,0,23245
    scaffold00001   16150   .   C   T   607.77  AllGtsFiltered  AC=2;AF=1.00;AN=2;DP=16;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=33.51;SOR=0.941    GT:AD:DP:FT:GQ:PL   1/1:0,16:16:LowDP:48:636,48,0
    

    As far as I can see most of these genotypes have DP<50 and it was (my intention at least) to have these filtered out. That does not seem to be the case.

    I'd appreciate feedback here.

    Thank you.
    Best
    Jahn

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭
    edited April 2017

    Hi Jahn (@jahneldar), your sites are being filtered! They have the AllGtsFiltered annotation in the FILTER field. You can then remove these records entirely from the file using GATK SelectVariants. Add the --excludeFiltered option to the SelectVariants command.

  • jahneldarjahneldar NorwayMember
    Accepted Answer

    Ahh, so another step required?
    I'm locked in vcftools.
    Thank you for getting back to me.

    jahn

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭
    Accepted Answer

    Depending on your downstream analysis, physically removing the filtered sites may not be necessary. That is, downstream tools will understand that these sites are filtered. Just annotating the FILTER field is the typical approach so that you don't lose data along the way. E.g. although sites are filtered, you may wish to look into borderline cases and so would want to keep them in your file.

Sign In or Register to comment.