This is a test version of Biostars. For the public version, visit https://www.biostars.org.
GATK: How to filter variant sites of sample VCF that are not present in a set of control VCFs

I have to select variants sites from my sample which are not present in the control samples. I have tried this:

java -jar GenomeAnalysisTK.jar \
  -R /Users/bhaskaran/Desktop/1data/GRCm38_68.fa \
  -T CombineVariants \
  -V:CON1 /Users/bhaskaran/Desktop/Tools/GenomeAnalysisTK-3.2-2/control/8404.vcf \
  -V:CON2 /Users/bhaskaran/Desktop/Tools/GenomeAnalysisTK-3.2-2/control/8405.vcf \
  -V:CON3 /Users/bhaskaran/Desktop/Tools/GenomeAnalysisTK-3.2-2/control/8406.vcf \
  -V:CON4 /Users/bhaskaran/Desktop/Tools/GenomeAnalysisTK-3.2-2/control/8407.vcf \
  -V:CON5 /Users/bhaskaran/Desktop/Tools/GenomeAnalysisTK-3.2-2/control/8408.vcf \
  -V:TED /Users/bhaskaran/Desktop/Tools/GenomeAnalysisTK-3.2-2/parental_22.vcf \
  -o result101.vcf && \
java -jar GenomeAnalysisTK.jar \
  -R /Users/bhaskaran/Desktop/1data/GRCm38_68.fa \
  -T SelectVariants -V result101.vcf \
  -select "set == 'TED'" -o result102.vcf

Unfortunately, this command outputs absolutely no variant. No error message occurred. Is some one could help me?

sequencing snp next-gen
• 6,406 views • link
updated by Ram Ram • written by ravast • ravast

2 answers

Does

awk '! /\#/' result101.vcf | grep -w "TED" | grep -vi "CON1\|CON2\|CON3\|CON4\|CON5\|Intersection" > TED_only.vcf

produce something?

• 0 views • link
updated by Ram Ram • written by Vivek Vivek

Vivek as per your suggestion I tried this:

java -jar GenomeAnalysisTK.jar \
  -R /Users/bhaskaran/Desktop/1data/GRCm38_68.fa \
  -T CombineVariants \
  -V:CON1 /Users/bhaskaran/Desktop/Tools/GenomeAnalysisTK-3.2-2/control/8404.vcf \
  -V:CON2 /Users/bhaskaran/Desktop/Tools/GenomeAnalysisTK-3.2-2/control/8405.vcf \
  -V:CON3 /Users/bhaskaran/Desktop/Tools/GenomeAnalysisTK-3.2-2/control/8406.vcf \
  -V:CON4 /Users/bhaskaran/Desktop/Tools/GenomeAnalysisTK-3.2-2/control/8407.vcf \
  -V:CON5 /Users/bhaskaran/Desktop/Tools/GenomeAnalysisTK-3.2-2/control/8408.vcf \
  -V:TED /Users/bhaskaran/Desktop/Tools/GenomeAnalysisTK-3.2-2/parental_22.vcf \
  -o result101.vcf -env

and then

awk '! /\#/' result101.vcf | grep -w "TED" | grep -vi "CON1\|CON2\|CON3\|CON4\|CON5\|Intersection" > TED_only.vcf

again I'm getting no result. But when I compared this file manually I could find many variants in my sample that are not found in controls

• 0 views • link
updated by Ram Ram • written by ravast • ravast

I'm not sure how you are doing a manual check. Do you want the variant to NOT be present in atleast one control or ANY control?

Can you post what you get from

awk '! /\#/' result101.vcf | cut -f8 | awk -F ';' '{print $NF}' | sort | uniq -c
• 0 views • link
written by Vivek Vivek

I want variant not be present in a tleast one control. I got this result..

3937677 set=Intersection
• 0 views • link
updated by Ram Ram • written by ravast • ravast

That pretty much tells you that you have 3,937,677 variants in the VCF file and all of them are common to all samples.

• 0 views • link
updated by Ram Ram • written by Vivek Vivek

OK. But I have a doubt. When I saw the result101.vcf for example at chrm 1 position -

3010834    .    A    G    63.13    .    AC=1;AF=0.500;AN=2;BaseQRankSum=-2.428;DP=182;Dels=0.00;FS=24.892;HaplotypeScore=1.6524;InbreedingCoeff=-0.2223;MQ=44.54;MQ0=1;MQRankSum=0.170;QD=5.26;ReadPosRankSum=0.852;set=Intersection    GT:AD:DP:GQ:PL    ./.    ./.    ./.    ./.    ./.    0/1:2,5:7:19:80,0,19.

Does that mean for my sample, the variant at position 3010834 is present whereas it is absent in the controls...sorry if I'm wrong...

• 0 views • link
updated by Ram Ram • written by ravast • ravast

That looks like that variant was not called in the first 5 samples but if that were the case, GATK shouldn't tell you that its in set=intersection. You could look at the individual control sample VCFs and check if that variant entry is present in them and if so what genotype it has.

• 0 views • link
updated by Ram Ram • written by Vivek Vivek

You are correct. I'm trying the way as you said..

• 0 views • link
updated by Ram Ram • written by ravast • ravast

Its not geting

• 0 views • link
updated by Ram Ram • written by ravast • ravast

I tried with individual file also ..But it is still showing no variants present.....

• 0 views • link
updated by Ram Ram • written by ravast • ravast

You can use bedools subtract.

If you have the control samples in multiple vcf files you can first combine those to one vcf file. And then use bedtools subtract.

• 0 views • link
updated by Ram Ram • written by William William

Williams i tried with bedttools..but its showing no new variant..

• 0 views • link
written by ravast • ravast

Log in to answer this question.

More posts like this

No similar posts found.