GATK SelectVariants behavior with gvcfs

Hi, I have a VCF resulting from running Mutect2 (joint calling) on multiple Tumor-only samples (using Sarek). The specific file has 2 samples and 1842870 calls. If I count all the PASS calls

cat test.vcf | grep PASS | wc -l

I get 39389

Now, if I run

gatk SelectVariants -V $dataFILE.vcf -O test.vcf

gatk SelectVariants -R $genFILE -V $dataFILE.vcf -O test.vcf

gatk SelectVariants -R $genFILE -V $dataFILE.vcf -O test.vcf --exclude-filtered

gatk SelectVariants -R $genFILE -V $dataFILE.vcf -O test.vcf --exclude-non-variants

gatk SelectVariants -R $genFILE -V $dataFILE.vcf -O test.vcf --exclude-filtered --exclude-non-variants

and I count all the PASS calls, I get, respectively 3326, 3533, 3407, 1116, 2975

If I run bcftools or vcftools (see below) I get 39390 in all the cases

bcftools view -o test1.vcf -f "PASS,." $dataFILE.vcf.gz

bcftools view -Ou -s SMP1 $dataFILE.vcf.gz | bcftools view -o test.vcf -f "PASS,."

bcftools view -Ou -s SMP2 $dataFILE.vcf.gz | bcftools view -o test.vcf -f "PASS,."

vcftools --gzvcf $dataFILE.vcf.gz --recode --recode-INFO-all --stdout --keep-filtered "PASS" --keep-filtered "." > test.vcf

vcftools --gzvcf $dataFILE.vcf.gz --recode --recode-INFO-all --stdout --keep-filtered "PASS" --keep-filtered "." --indv SMP1 > test.vcf

vcftools --gzvcf $dataFILE.vcf.gz --recode --recode-INFO-all --stdout --keep-filtered "PASS" --keep-filtered "." --indv SMP2 > test.vcf

I must be doing something stupid. Can someone comment ?

NOTE: I am running GATK v4.6.0.0

Thanks

gatk selectvariants gvcfs

• 918 views • link

updated 8 months ago by Jeremy Leipzig

• written 8 months ago by Gianfilippo •

cat test.vcf | grep PASS | wc -l

hum.. hard to say without seeing the VCF itself but the genotypes may contain PASS for FORMAT/FT. I would count the variants using

bcftools view --apply-filters 'PASS,.' -O u test.vcf | bcftools query -f '\n' | wc -l

• 0 views • link

written 8 months ago by Pierre Lindenbaum

1 answer

Sarek does not make it particularly easy to obtain gVCFs, so i wonder that's what you really meant.

Show us an output line from one of your failing commands you didn't expect (something that isn't either PASS or .)

• 0 views • link

written 8 months ago by Jeremy Leipzig

Hi, below are a few PASS calls from one individual with 2 samples. I should add that now the gatk selectvariants shows a consistent behavior. So I assume something was somehow wrong before. Nevertheless, now I get the same N of PASS calls in all runs of gatk selectvariants. There seems to be no difference in call set between the 2 samples, which seems unlikely.

chr1    17904   .   G   A   .   PASS    AS_FilterStatus=SITE;AS_SB_TABLE=0,43|0,5;DP=51;ECNT=1;ECNTH=1;GERMQ=37;MBQ=31,30;MFRL=471,475;MMQ=40,40;MPOS=4;POPAF=2.33;ROQ=60;TLOD=7.79GT:AD:AF:DP:F1R2:F2R1:FAD:SB 0/1:24,3:0.137:27:11,2:12,1:24,3:0,24,0,3   0/1:19,2:0.129:21:9,1:10,1:19,2:0,19,0,2

chr1    183189  .   G   C   .   PASS    AS_FilterStatus=SITE;AS_SB_TABLE=29,12|6,0;DP=51;ECNT=2;ECNTH=2;GERMQ=75;MBQ=30,24;MFRL=476,473;MMQ=39,40;MPOS=7;POPAF=7.30;ROQ=53;TLOD=9.7GT:AD:AF:DP:F1R2:F2R1:FAD:SB 0/1:23,3:0.141:26:16,3:7,0:23,3:17,6,3,0    0/1:18,3:0.172:21:11,0:7,3:18,3:12,6,3,0

chr1    183238  .   G   C   .   PASS    AS_FilterStatus=SITE;AS_SB_TABLE=121,15|16,2;DP=163;ECNT=2;ECNTH=2;GERMQ=93;MBQ=33,33;MFRL=470,470;MMQ=40,40;MPOS=24;POPAF=7.30;ROQ=93;TLOD=36.99   GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:72,10:0.130:82:48,4:24,6:72,10:61,11,8,2    0/1:64,8:0.121:72:33,6:31,2:64,8:60,4,8,0

chr1    185968  .   A   G   .   PASS    AS_FilterStatus=SITE;AS_SB_TABLE=7,0|10,0;DP=18;ECNT=1;ECNTH=1;GERMQ=13;MBQ=36,27;MFRL=473,469;MMQ=40,40;MPOS=16;POPAF=7.30;ROQ=85;TLOD=22.82   GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:2,3:0.572:5:2,3:0,0:2,3:2,0,3,0 0/1:5,7:0.571:12:2,3:3,4:5,7:5,0,7,0

chr1    188428  .   G   A   .   PASS    AS_FilterStatus=SITE;AS_SB_TABLE=0,0|0,0;DP=2;ECNT=1;ECNTH=1;GERMQ=20;MBQ=0,39;MFRL=0,465;MMQ=60,32;MPOS=17;POPAF=7.30;ROQ=64;TLOD=6.61 GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:0,1:0.667:1:0,1:0,0:0,1:0,0,0,1 0/1:0,1:0.667:1:0,0:0,1:0,1:0,0,0,1

chr1    268054  .   A   G   .   PASS    AS_FilterStatus=SITE;AS_SB_TABLE=63,94|21,11;DP=196;ECNT=1;ECNTH=1;GERMQ=93;MBQ=33,31;MFRL=463,466;MMQ=40,40;MPOS=17;POPAF=0.876;ROQ=93;TLOD=62.11  GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:93,25:0.216:118:52,9:40,16:93,25:38,55,15,10    0/1:64,7:0.109:71:30,3:33,4:64,7:25,39,6,1

chr1    268833  .   T   A   .   PASS    AS_FilterStatus=SITE;AS_SB_TABLE=44,108|10,26;DP=203;ECNT=1;ECNTH=1;GERMQ=93;MBQ=33,30;MFRL=469,468;MMQ=40,40;MPOS=20;POPAF=7.30;ROQ=93;TLOD=76.93  GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:108,21:0.167:129:50,8:58,13:108,21:37,71,7,14   0/1:44,15:0.261:59:21,6:23,9:44,15:7,37,3,12

chr1    517599  .   A   AAC .   PASS    AS_FilterStatus=SITE;AS_SB_TABLE=0,2|3,6;DP=14;ECNT=1;ECNTH=1;GERMQ=8;MBQ=33,35;MFRL=455,457;MMQ=47,47;MPOS=38;POPAF=7.30;ROQ=93;RPA=15,16;RU=AC;STR;STRQ=6;TLOD=25.69  GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:1,6:0.786:7:1,3:0,3:1,6:0,1,3,3 0/1:1,3:0.667:4:0,2:1,1:1,3:0,1,0,3

• 0 views • link

updated 8 months ago by GenoMax

• written 8 months ago by Gianfilippo •

I wouldn't expect differences in high-quality genotype calls from two dna samples from the same subject.

Please show us a line you expected to be filtered out that was not filtered out.

• 0 views • link

written 8 months ago by Jeremy Leipzig