Notice:
If you happen to see a question you know the answer to, please do chime in and help your fellow community members. We encourage our fourm members to be more involved, jump in and help out your fellow researchers with their questions. GATK forum is a community forum and helping each other with using GATK tools and research is the cornerstone of our success as a genomics research community.We appreciate your help!

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.
Attention:
We will be out of the office on November 11th and 13th 2019, due to the U.S. holiday(Veteran's day) and due to a team event(Nov 13th). We will return to monitoring the GATK forum on November 12th and 14th respectively. Thank you for your patience.

A question about SelectVariants

zhoujiayizhoujiayi CanadaMember
edited July 2015 in Ask the GATK team

I am using the SelectVariants command to compare to samples.
Such as
java -Xmx2g -jar /home/zhoujiayi/GenomeAnalysisTK/GenomeAnalysisTK.jar \
-R /home/REF_FILES/HG19/hg19.fasta \
-T SelectVariants \
-V /home/a.vcf \
--discordance /home/b.vcf \
-o /home/UniqueInA.vcf

java -Xmx2g -jar /home/zhoujiayi/GenomeAnalysisTK/GenomeAnalysisTK.jar \
-R /home/zhoujiayi/REF_FILES/HG19/hg19.fasta \
-T SelectVariants \
-V /home/b.vcf \
--discordance /b.vcf \
-o /home/UniqueInB.vcf

java -Xmx2g -jar /home/zhoujiayi/GenomeAnalysisTK/GenomeAnalysisTK.jar \
-R /home/zhoujiayi/REF_FILES/HG19/hg19.fasta \
-T SelectVariants \
-V /home/a.vcf \
--concordance /home/b.vcf \
-o /home/concordance.vcf

Then, I checked the number of variants in each vcf file and got the following result:
Number of variants in a.vcf is 53502.
Number of variants in b.vcf is 52665.
Number of variants in concordance.vcf is 48765.
Number of variants in UniqueInA.vcf is 4737.
Number of variants in UniqueInB.vcf is 3876.

53502=48765+4737, but 48765+3876=52641( should be 52665, 24 variants are missing.)

Have you got any ideal how did this happen?

Answers

  • SheilaSheila Broad InstituteMember, Broadie admin

    @zhoujiayi
    Hi,

    I suspect the missing sites have a mix of indel and SNP for alternate allele. Can you post a few records of the missing sites?

    Thanks,
    Sheila

  • zhoujiayizhoujiayi CanadaMember

    I cannot provide the missing sites. Because if I use GATK command selectVariants to find out the missing sites, the result is an empty vcf file. I think this means in a.vcf, there exists some duplicated variants? But if so, how come duplicated variants are there? @Sheila

  • SheilaSheila Broad InstituteMember, Broadie admin
    edited July 2015

    @zhoujiayi
    Hi,

    Can you try running Select Variants with -selectType MIXED on b.vcf? If there are 24 sites output. we will know the missing sites are indeed mixed record sites. https://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_gatk_tools_walkers_variantutils_SelectVariants.php#--selectTypeToInclude

    Thanks,
    Sheila

  • zhoujiayizhoujiayi CanadaMember
    edited July 2015

    Hi, Sheila:
    @Sheila
    I run the -selectType MIXED on b.vcf, but I got a vcf file with 406 variants... The command I ran is:
    java -Xmx2g -jar /home/zhoujiayi/GenomeAnalysisTK/GenomeAnalysisTK.jar -R /home/zhoujiayi/REF_FILES/HG19/hg19.fasta -T SelectVariants -V /home/b.vcf -selectType MIXED -o /home/BMIXED.vcf

  • SheilaSheila Broad InstituteMember, Broadie admin

    @zhoujiayi
    Hi,

    Did you do any filtering of the variants? Are there any filter names applied to your VCFs?
    I am wondering if the excluded sites are filtered out.

    Thanks,
    Sheila

  • zhoujiayizhoujiayi CanadaMember

    @Sheila
    I didn't do any filtering to these file. I got them form the Ion Torrent server. I am not sure if Ion Torrent did any filtering at the variant calling step. But I think it doesn't matter, right? Even these vcf files have been filtered, as long as they are filtered by the same options, it should be fine.
    I just want to compare the difference in detected variants between these 2 vcf files.

  • SheilaSheila Broad InstituteMember, Broadie admin

    @zhoujiayi
    Hi,

    It does matter because sometimes if a lowQual or filtered site is labeled as such, it will not be included in the considered sites. In your case, I was wondering if the B.vcf contained some filtered sites that were not filtered in A.vcf. Can you try one more thing? Check and see how many MIXED sites are in A?

    I have been trying this out myself on a small test set, and I get the correct numbers of variants in each output file. How did you calculate the number of variants in each file?

    Thanks,
    Sheila

  • zhoujiayizhoujiayi CanadaMember

    @Sheila
    I checked that there are 1033 MIXED variants in a.vcf.

    I have tested 2 sets of vcf files, the number of variants matches for the other set. So I am really curious what causes this problem.

    I just counted the lines in the vcf file. Each variant is listed in one line. Of course, I didn't count the headers.

    Thank you.

  • zhoujiayizhoujiayi CanadaMember
    edited July 2015

    @Sheila
    I wrote an R function to find the duplicated sites in b.vcf. I found 28 sites are both SNPs and INDELs. I think those missing 24 variants are among these 28. The reason MIXED found much more than this number is because I think it founds the variants like following:(not real data)

    chr1 45472818 . G T, TGC 452.58 PASS

    But I am looking for the variants like:
    chr1 45472818 . G T 452.58 PASS
    chr1 45472818 . GTGTGTGTGTGTT G 192.67 PASS

    I am not sure if my understanding is correct.

    But even for the above cases, I think the numbers should still match. I am now working on to get the 24 missing variants. I will let you know the result once I get them.

  • zhoujiayizhoujiayi CanadaMember
    edited July 2015

    @Sheila

    I found the variants that exists in b.vcf but not in either concordance.vcf nor UniqueInB.vcf. They are among the 28 mixed sites, such as:

    chr1 85331808 . G GAGA 32.05 PASS AF=0.444444;AO=27;DP=64;FAO=28;FDP=63;FR=.;FRO=35;FSAF=7;FSAR=21;FSRF=7;FSRR=28;FWDB=-0.00777824;FXX=0.0156226;HRUN=1;LEN=3;MLLD=246.5;PB=0.5;PBP=1;QD=2.03488;RBI=0.0205399;REFB=-0.0229589;REVB=-0.0190102;RO=36;SAF=7;SAR=20;SRF=8;SRR=28;SSEN=0;SSEP=0;SSSB=0.0500128;STB=0.538435;STBP=0.668;TYPE=ins;VARB=0.0408464;OID=.;OPOS=85331809;OREF=-;OALT=AGA;OMAPALT=GAGA GT:AF:AO:DP:FAO:FDP:FRO:FSAF:FSAR:FSRF:FSRR:GQ:RO:SAF:SAR:SRF:SRR 0/1:0.444444:27:64:28:63:35:7:21:7:28:32:36:7:20:8:28
    chr1 85331808 . G GAGA 307.44 PASS AF=0.966667;AO=29;DP=64;FAO=58;FDP=60;FR=.;FRO=2;FSAF=13;FSAR=45;FSRF=1;FSRR=1;FWDB=0.0151638;FXX=0.0624902;HRUN=1;LEN=3;MLLD=230.796;PB=0.5;PBP=1;QD=20.4957;RBI=0.0287885;REFB=-0.0487047;REVB=0.0244712;RO=31;SAF=6;SAR=23;SRF=9;SRR=22;SSEN=0;SSEP=0;SSSB=-0.108005;STB=0.513025;STBP=0.385;TYPE=ins;VARB=0.00192847;OID=.;OPOS=85331809;OREF=-;OALT=AGA;OMAPALT=GAGA GT:AF:AO:DP:FAO:FDP:FRO:FSAF:FSAR:FSRF:FSRR:GQ:RO:SAF:SAR:SRF:SRR 1/1:0.966667:29:64:58:60:2:13:45:1:1:48:31:6:23:9:22
    chr1 224302154 . T TGAG 43.05 PASS AF=0.390533;AO=66;DP=178;FAO=66;FDP=169;FR=.;FRO=103;FSAF=23;FSAR=43;FSRF=42;FSRR=61;FWDB=0.00387499;FXX=0;HRUN=1;LEN=3;MLLD=260.332;PB=0.5;PBP=1;QD=1.01893;RBI=0.0150409;REFB=-0.000127504;REVB=-0.0145332;RO=112;SAF=23;SAR=43;SRF=51;SRR=61;SSEN=0;SSEP=0;SSSB=-0.12452;STB=0.538835;STBP=0.433;TYPE=ins;VARB=0.000269841;OID=.;OPOS=224302155;OREF=-;OALT=GAG;OMAPALT=TGAG GT:AF:AO:DP:FAO:FDP:FRO:FSAF:FSAR:FSRF:FSRR:GQ:RO:SAF:SAR:SRF:SRR 0/1:0.390533:66:178:66:169:103:23:43:42:61:43:112:23:43:51:61
    chr1 224302154 . T TGAG 791.97 PASS AF=0.950617;AO=85;DP=178;FAO=154;FDP=162;FR=.;FRO=8;FSAF=57;FSAR=97;FSRF=5;FSRR=3;FWDB=0.00100543;FXX=0.0414177;HRUN=1;LEN=3;MLLD=200.978;PB=0.5;PBP=1;QD=19.5549;RBI=0.00101569;REFB=-0.00266943;REVB=0.000144001;RO=83;SAF=31;SAR=54;SRF=37;SRR=46;SSEN=0;SSEP=0;SSSB=-0.0767129;STB=0.513401;STBP=0.148;TYPE=ins;VARB=0.000132306;OID=.;OPOS=224302155;OREF=-;OALT=GAG;OMAPALT=TGAG GT:AF:AO:DP:FAO:FDP:FRO:FSAF:FSAR:FSRF:FSRR:GQ:RO:SAF:SAR:SRF:SRR 1/1:0.950617:85:178:154:162:8:57:97:5:3:99:83:31:54:37:46

    Above 4 variants only exists in b.vcf. However, among these 28 sites, only 24 are missing, other 4 are still listed in either concordance.vcf or UniqueInB.vcf. Such as:

    chr1 45472818 . G T 452.58 PASS AF=0.597403;AO=66;DP=128;FAO=92;FDP=154;FR=.;FRO=62;FSAF=83;FSAR=9;FSRF=53;FSRR=9;FWDB=0.0177401;FXX=0.0128197;HRUN=1;LEN=1;MLLD=126.56;PB=0.5;PBP=1;QD=11.7554;RBI=0.0265047;REFB=-0.0907561;REVB=0.0196923;RO=62;SAF=61;SAR=5;SRF=53;SRR=9;SSEN=0;SSEP=0;SSSB=0.178682;STB=0.549657;STBP=0.337;TYPE=snp;VARB=0.0782834;OID=.;OPOS=45472818;OREF=G;OALT=T;OMAPALT=T GT:GQ:DP:FDP:RO:FRO:AO:FAO:AF:SAR:SAF:SRF:SRR:FSAR:FSAF:FSRF:FSRR 0/1:99:128:154:62:62:66:92:0.597403:5:61:53:9:9:83:53:9
    chr1 45472818 . GTGTGTGTGTGTT G 192.67 PASS AF=0.692308;AO=60;DP=106;FAO=72;FDP=104;FR=.;FRO=32;FSAF=64;FSAR=8;FSRF=19;FSRR=13;FWDB=0.0141152;FXX=0.301993;HRUN=1;LEN=12;MLLD=289.773;PB=0.5;PBP=1;QD=7.41044;RBI=0.0389814;REFB=-0.0441051;REVB=-0.036336;RO=24;SAF=59;SAR=1;SRF=16;SRR=8;SSEN=0;SSEP=0;SSSB=0.670228;STB=0.669266;STBP=0.001;TYPE=del;VARB=0.0140532;OID=.;OPOS=45472819;OREF=TGTGTGTGTGTT;OALT=-;OMAPALT=G GT:AF:AO:DP:FAO:FDP:FRO:FSAF:FSAR:FSRF:FSRR:GQ:RO:SAF:SAR:SRF:SRR 0/1:0.692308:60:106:72:104:32:64:8:19:13:9:24:59:1:16:8

    The above 2 variants exists in b.vcf and UniqueInB.vcf.

    Do you know the reason for this?

    Thank you.

  • SheilaSheila Broad InstituteMember, Broadie admin

    @zhoujiayi
    Hi,

    Ah. I see. So, some of the sites may be duplicated because the Ion Torrent file has a separate line for SNP and indel at the same site? GATK seems to be doing the correct thing here. It only considers one site once in the analysis. Unfortunately, because the VCFs were not produced by GATK, we cannot really help you out more. Our tools are not designed to analyze VCFs that deviate from our VCF standard.

    -Sheila

  • zhoujiayizhoujiayi CanadaMember

    @Sheila

    I understand that GATK only considers one site once in the analysis. But why it ignores sites such as position 85331808 in both concordance.vcf and UniqueInB.vcf? If GATK consider it once, should it be at least in one vcf file?

    Thank you.

  • SheilaSheila Broad InstituteMember, Broadie admin

    @zhoujiayi
    Hi,

    What happens when you run Select Variants with -concordance and -discordance on that site alone?

    Thanks,
    Sheila

  • zhoujiayizhoujiayi CanadaMember

    @Sheila
    How to run Select Variants on that site alone? I don't find any parameter refer to specify a position.

    Thank you.

  • zhoujiayizhoujiayi CanadaMember

    @Sheila
    Thank you. I made a mistake. Yes, the missing variants are among those 28 mixed sites. For -concordance, GATK only count the same site once. But for -discordance, GATK will count all the same sites. That's why I only missed 24 variants. Those are the variants in concordance.

    Thank you.

  • SheilaSheila Broad InstituteMember, Broadie admin

    @zhoujiayi
    Hi,

    That is wonderful you figured it out!

    -Sheila

Sign In or Register to comment.