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.

genotypeConcordance for GVCF file

nitinCelmatixnitinCelmatix NYCMember
edited May 2016 in Ask the GATK team

Hi,
I need to compare two gVCF file using genotypeConcordance tool. In spite of them being the same, genotypeConcordance lists one difference between them:

EVAL.gVCF

1       11854456        rs1801131       T       G       1484.77 .       DB;DP=35;MLEAC=2,0;MLEAF=1.00,0.00;MQ=60.00;MQ0=0       GT:AD:DP:GQ:PGT:PID:PL:SB       1/1:0,35,0:35:99:0|1:11854457_G_A:1513,105,0,1513,105,1513:0,0,14,21
1       11856348        .       N       G       .       .       .       GT:DP:GQ:MIN_DP:PL      0/0:45:99:38:0,102,1141

Comp.gVCF

1       11854456        .       T       G       6189.77 .       DP=141;MLEAC=2,0;MLEAF=1.00,0.00;MQ=60.00       GT:AD:DP:GQ:PGT:PID:PL:SB       1/1:0,141,0:141:99:0|1:11854457_G_A:6218,424,0,6218,424,6218:0,0,66,75
1       11856348        .       N       G       .       .       .       GT:DP:GQ:MIN_DP:PL      0/0:57:99:48:0,102,1800
ALLELES_MATCH  EVAL_SUPERSET_TRUTH  EVAL_SUBSET_TRUTH  ALLELES_DO_NOT_MATCH  EVAL_ONLY  TRUTH_ONLY
            1                    0                  0                     0          0           1

I wonder why that is. Is there a way to run genotypeConcordance on the GVCF file?

Please note that I replaced the ref allele with 'N' and the alt allele with G (which was the ref allele) because otherwise GATK throws the error of duplicate alleles.

Any help would be appreciated.
Thanks!

Post edited by Geraldine_VdAuwera on

Issue · Github
by Sheila

Issue Number
924
State
closed
Last Updated
Assignee
Array
Milestone
Array
Closed By
vdauwera

Answers

  • nitinCelmatixnitinCelmatix NYCMember
    edited May 2016

    Another update:
    Even when I switch the EVAL and COMP files, the last line in the resulting output file remains the same:

    ALLELES_MATCH EVAL_SUPERSET_TRUTH EVAL_SUBSET_TRUTH ALLELES_DO_NOT_MATCH EVAL_ONLY TRUTH_ONLY
    1 0 0 0 0 1

    Should the last 2 fields (EVAL_ONLY and TRUTH_ONLY) switch upon switching the EVAL and COMP files?
    Thanks!

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    How are these GVCFs? There is no NON-REF allele? Also, what do you mean when you say you replaced the alleles?

  • nitinCelmatixnitinCelmatix NYCMember

    I need to prepare a VCF (or gVCF) for all variants of interest (even when they are not present) so that every file has equal number of lines (the number of variants of interest) and then I need to compare them for concordance. So, I extracted them from gVCF using gvcftools.
    I modified the GVCFs a little bit to make them closer to VCFs. So as you can see by the snapshot above, first, I removed the ', from the variant line. I also replaced the to '.' for the non-variant positions. That also did not give me correct number of matches from genotypeConcordance. So, then I replaced the '.' by the ref allele in the non-variant line. That gave the error of VCF being malformed due to duplicate allele for ref and alt. So then I changed the ref allele to 'N' and the alt allele to the ref allele for non-variant to trick genotypeConcordance into running. It did run but as you can see above, inspite of the 2 lines being the same, it only reports 1 matching allele.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @nitinCelmatix
    Hi,

    Perhaps the best thing to do is run GenotypeGVCFs with -allSites to produce a proper VCF with all the sites in it. Then, you can run GenotypeConcordance.

    -Sheila

  • nitinCelmatixnitinCelmatix NYCMember

    Thanks Sheila,
    That did not work either. I first created a smaller gvcf using gvcf tools, then created a VCF with GATK using -allSites and then ran genptypeConcordance but still the number of TRUTH_ONLY alleles is not correct.
    Any help would be greatly appreciated.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @nitinCelmatix
    Hi,

    I just tested this out, and I get the correct results. Can you please post the exact command you ran and what version of GATK you are using? Also, can you test this on a single site (from the final -allSites VCF) where you switch the eval and comp file and post the output here?

    Thanks,
    Sheila

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @nitinCelmatix We can't provide any support if you're running on files that have been modified in this manner. Also, you shouldn't be using gvcftools on our GVCFs. That package was developed to work with Illumina's version of GVCFs, which are quite different, and will not do the right thing with our GVCFs. If you want to generate all-sites VCFs, you need to run GenotypeGVCFs with the -allSites argument as Sheila recommended. What you're doing right now is simply not supported and we won't help you deal with what comes out of there.

  • nitinCelmatixnitinCelmatix NYCMember

    Hi Sheila,
    Here are the commands:

    first create a smaller gVCF that only contains the variant positions in the variants.bed file for sample1:

    gzip -dc sample1.g.vcf.gz | gvcftools-0.17.0/bin/break_blocks --region-file variants.bed --ref hg19.fa --exclude-off-target > sample1.variants.g.vcf

    Then create the VCF:

    java -jar GenomeAnalysisTK.jar -T GenotypeGVCFs -R hg19.fa -V sample1.variants.g.vcf -allSites -o sample1.vcf

    Similarly, create the smaller gVCF for sample2:

    cat sample2.g.vcf | gvcftools-0.17.0/bin/break_blocks --region-file variants.bed --ref hg19.fa --exclude-off-target > sample2.variants.g.vcf

    and ceate the VCF:
    java -jar GenomeAnalysisTK.jar -T GenotypeGVCFs -R hg19.fa -V sample2.variants.g.vcf -allSites -o sample2.variants.vcf

    Calculate concordance:
    java -jar GenomeAnalysisTK.jar -T GenotypeConcordance -R hg19.fa -eval sample2.variants.vcf -comp sample1.vcf -o test.txt

    $tail test.txt

    ALLELES_MATCH EVAL_SUPERSET_TRUTH EVAL_SUBSET_TRUTH ALLELES_DO_NOT_MATCH EVAL_ONLY TRUTH_ONLY
    23 0 0 0 0 26

    I am using GATK 3.4

    Then I only keep the following 2 entries in the VCF files (first one is a variant, the 2nd is a non-variant):

    1 11854446 . T G 6189.77 . AC=2;AF=1.00;AN=2;DP=141;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=34.24;SOR=0.827 GT:AD:DP:GQ:PGT:PID:PL 1/1:0,141:141:99:1|1:11854457_G_A:6218,424,0
    1 11856338 . G . . . AN=2;DP=48 GT:DP:RGQ 0/0:48:99

    they both occur in sample1 and sample2 vcfs.

    When I run either of the following commands by switching eval and comp, the output file is same:
    java -jar GenomeAnalysisTK.jar -T GenotypeConcordance -R hg19.fa -comp sample1.vcf -eval sample2.vcf -o test.txt

    java -jar GenomeAnalysisTK.jar -T GenotypeConcordance -R hg19.fa -comp sample2.vcf -eval sample1.vcf -o test.txt

    These are the last 2 lines of test.txt in either of the above 2 cases (it says 1 TRUTH_ONLY allele was found inspite of the 2 files being the same):

    ALLELES_MATCH EVAL_SUPERSET_TRUTH EVAL_SUBSET_TRUTH ALLELES_DO_NOT_MATCH EVAL_ONLY TRUTH_ONLY
    1 0 0 0 0 1

    Next, when I remove the 2nd of the two entries in the VCF file (the non-variant one), and only keep the variant one, I get correct number. Output looks like:
    ALLELES_MATCH EVAL_SUPERSET_TRUTH EVAL_SUBSET_TRUTH ALLELES_DO_NOT_MATCH EVAL_ONLY TRUTH_ONLY
    1 0 0 0 0 0

    So, I think the non-variant in the VCF file is throwing GATK genotypeConcordance off.

    Hope this helps you to diagnose the problem.
    Thanks for any help.

  • nitinCelmatixnitinCelmatix NYCMember

    Hi,
    I was wondering if you got a chance to look into the inaccurate number that I am getting from genotypeConcordance.
    One quick thing that you can do to reproduce this at your end is to just use the two lines I pasted above in your vcf files.
    Thanks for all your help!
    Best!

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    GenotypeConcordance only considers variant sites, ie sites where there is at least one ALT allele. It is not designed to operate on non-variant sites. To check whether non-variant sites are concordant between two callsets you can use the site-level concordance functionality of VariantEval.

  • robertbrobertb torontoMember

    Hi,

    I'm some concordance with GATK 3.6 and get the following error:

    ERROR MESSAGE: Concordance Metrics is currently only implemented for DIPLOID genotypes, found eval ploidy: 1, comp ploidy: 0

    Not sure how to fix this. Does someone know what's going on? Thanks - Robert

    [[email protected] concordance]$ ./concordance.sh
    INFO 15:50:21,801 HelpFormatter - ----------------------------------------------------------------------------------
    INFO 15:50:21,803 HelpFormatter - The Genome Analysis Toolkit (GATK) v3.6-0-g89b7209, Compiled 2016/06/01 22:27:29
    INFO 15:50:21,804 HelpFormatter - Copyright (c) 2010-2016 The Broad Institute
    INFO 15:50:21,804 HelpFormatter - For support and documentation go to https://www.broadinstitute.org/gatk
    INFO 15:50:21,804 HelpFormatter - [Mon Jul 09 15:50:21 EDT 2018] Executing on Linux 2.6.32-696.18.7.el6.x86_64 amd64
    INFO 15:50:21,805 HelpFormatter - Java HotSpot(TM) 64-Bit Server VM 1.8.0_91-b14 JdkDeflater
    INFO 15:50:21,808 HelpFormatter - Program Args: -R /hpf/tools/centos6/bcbio/20180613/genomes/Hsapiens/GRCh37/seq/GRCh37.fa --comp:VCF /hpf/largeprojects/pray/rbaldwin/WGS/concordance/202214_S1.genome.PASS.header.snps.vcf --eval:VCF /hpf/largeprojects/pray/rbaldwin/WGS/concordance/wgs_test-gatk-haplotype-annotated-snps.vcf.gz -o /hpf/largeprojects/pray/rbaldwin/WGS/concordance/CONORDANCE_.eval.snps.txt -T GenotypeConcordance -L /hpf/largeprojects/pray/llau/internal_databases/baits/SS_clinical_research_exomes/S06588914/S06588914_Covered.sort.merged.bed
    INFO 15:50:21,814 HelpFormatter - Executing as [email protected] on Linux 2.6.32-696.18.7.el6.x86_64 amd64; Java HotSpot(TM) 64-Bit Server VM 1.8.0_91-b14.
    INFO 15:50:21,815 HelpFormatter - Date/Time: 2018/07/09 15:50:21
    INFO 15:50:21,815 HelpFormatter - ----------------------------------------------------------------------------------
    INFO 15:50:21,816 HelpFormatter - ----------------------------------------------------------------------------------
    INFO 15:50:21,835 GenomeAnalysisEngine - Strictness is SILENT
    INFO 15:50:21,996 GenomeAnalysisEngine - Downsampling Settings: Method: BY_SAMPLE, Target Coverage: 1000
    INFO 15:50:51,579 RMDTrackBuilder - Writing Tribble index to disk for file /hpf/largeprojects/pray/rbaldwin/WGS/concordance/202214_S1.genome.PASS.header.snps.vcf.idx
    INFO 15:50:58,028 IntervalUtils - Processing 54098923 bp from intervals
    WARN 15:50:58,044 IndexDictionaryUtils - Track eval doesn't have a sequence dictionary built in, skipping dictionary validation
    INFO 15:50:58,152 GenomeAnalysisEngine - Preparing for traversal
    INFO 15:50:58,195 GenomeAnalysisEngine - Done preparing for traversal
    INFO 15:50:58,196 ProgressMeter - [INITIALIZATION COMPLETE; STARTING PROCESSING]
    INFO 15:50:58,196 ProgressMeter - | processed | time | per 1M | | total | remaining
    INFO 15:50:58,197 ProgressMeter - Location | sites | elapsed | sites | completed | runtime | runtime
    INFO 15:51:28,245 ProgressMeter - 7:106786952 18993.0 30.0 s 26.4 m 39.9% 75.0 s 45.0 s

    ERROR ------------------------------------------------------------------------------------------
    ERROR A USER ERROR has occurred (version 3.6-0-g89b7209):
    ERROR
    ERROR This means that one or more arguments or inputs in your command are incorrect.
    ERROR The error message below tells you what is the problem.
    ERROR
    ERROR If the problem is an invalid argument, please check the online documentation guide
    ERROR (or rerun your command with --help) to view allowable command-line arguments for this tool.
    ERROR
    ERROR Visit our website and forum for extensive documentation and answers to
    ERROR commonly asked questions https://www.broadinstitute.org/gatk
    ERROR
    ERROR Please do NOT post this error to the GATK forum unless you have really tried to fix it yourself.
    ERROR
    ERROR MESSAGE: Concordance Metrics is currently only implemented for DIPLOID genotypes, found eval ploidy: 1, comp ploidy: 0
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    I would recommend checking your input files to make sure they all have diploid genotypes in them.

  • robertbrobertb torontoMember

    brilliant! thanks.

Sign In or Register to comment.