Difference in GenotypeGVCFs generated VCF after consolidation with GenomicsDBimport and CombineGVCF

prasundutta87prasundutta87 EdinburghMember

Hi,

I had a set of total 81 GVCFs that I first consolidated using GenomicsDBimport and then using CombineGVCF and then GenotypeGVCF was run in both cases. For GenomicsDBimport, I ran the command per contig and then I ran GenotypeGVCF on each database to get the final VCF file. Then I used Picard GatherGVCF to make the final VCF. The commands I used are wriiten below:

Using GenomicsDBimport:
java -Xmx90g -jar gatk-package-4.0.4.0-local.jar GenomicsDBImport -R water_buffalo_re_arranged_chrom_ref_genome.fa --TMP_DIR ./tmp --sample-name-map sample_names_map_new.txt --reader-threads 2 --genomicsdb-workspace-path "$contig" -L "$contig"

java -Xmx8G -XX:ConcGCThreads=1 -jar gatk-package-4.0.4.0-local.jar GenotypeGVCFs -R /water_buffalo_re_arranged_chrom_ref_genome.fa -new-qual -V gendb://"$contig" -O "$contig"_variants.vcf.gz

java -jar picard.jar GatherVcfs INPUT=list.txt OUTPUT=Final_med_buffalo_variants_81_samples.vcf.gz

Using CombineGVCF:
java -Xmx200g -XX:ConcGCThreads=1 -jar gatk-package-4.0.4.0-local.jar CombineGVCFs -R water_buffalo_re_arranged_chrom_ref_genome.fa --variant All_gvcf_gz.list -O combined_81.g.vcf.gz

java -Xmx8G -XX:ConcGCThreads=1 -jar gatk-package-4.0.4.0-local.jar GenotypeGVCFs -R water_buffalo_re_arranged_chrom_ref_genome.fa -new-qual -V combined_81.g.vcf.gz -O Final_variants_81_samples_using_CombineGVCF.vcf.gz

The final VCF in both the cases should be the same. Unfortunately, it was not. On running bcftools isec, I found that some variants were common to one VCF and some were in other. What could be the reason behind this discrepancy?

Kindly let me know if you need more information.

Answers

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @prasundutta87
    Hi,

    Can you post some sites that were discrepant? How many sites were there that were different between the two?

    Thanks,
    Sheila

  • prasundutta87prasundutta87 EdinburghMember
    edited June 2018

    Sure..

    Records private to VCF generated using CombineGVCF: 2313
    Records private to VCF generated using GenomicsDBimport: 2312

    Some example descripant sites (just posting site information and INFO fields):

    From CombineGVCF:

    CHROM POS ID REF ALT QUAL FILTER INFO

    1 314788 . CTTTT CT,CTT,C,CTTTTTTTTT,CTTTTTTTTTTTTTTTT,CTTTTTT 22673.7 . AC=48,28,2,10,4,3;AF=0.393,0.23,0.016,0.082,0.033,0.025;AN=122;BaseQRankSum=0.454;ClippingRankSum=0;DP=1140;ExcessHet=5.2228;FS=28.037;InbreedingCoeff=0.18;MLEAC=68,39,3,13,6,4;MLEAF=0.557,0.32,0.025,0.107,0.049,0.033;MQ=58.08;MQRankSum=0;QD=33.89;ReadPosRankSum=0.449;SOR=2.457
    1 473088 . TAAAA TAAAAA,TAAAAAA,TA,TAA,TAAAAAAAAAAAAAA,T 19215.5 . AC=12,23,28,12,18,3;AF=0.094,0.18,0.219,0.094,0.141,0.023;AN=128;BaseQRankSum=0;ClippingRankSum=0;DP=1059;ExcessHet=4.0535;FS=0;InbreedingCoeff=0.396;MLEAC=12,33,40,17,25,5;MLEAF=0.094,0.258,0.313,0.133,0.195,0.039;MQ=58.15;MQRankSum=0;QD=31.84;ReadPosRankSum=0.172;SOR=0.63

    From GenomicsDBimport:

    CHROM POS ID REF ALT QUAL FILTER INFO

    1 314788 . CTTT C,CT,CTTTTTTTTTTTTTT,CTTTTTTTTTTTTTTT,CTTTTT,CTTTTTTTT 22806.1 . AC=50,28,3,2,3,10;AF=0.41,0.23,0.025,0.016,0.025,0.082;AN=122;BaseQRankSum=0.454;ClippingRankSum=0;DP=1140;ExcessHet=4.7598;FS=28.037;InbreedingCoeff=0.193;MLEAC=71,39,4,4,4,13;MLEAF=0.582,0.32,0.033,0.033,0.033,0.107;MQ=58.08;MQRankSum=0;QD=30.77;ReadPosRankSum=0.449;SOR=2.457
    1 473088 . TAAA TAAAAAAAAAAAAA,TAAAA,T,TA,TAAAAA,TAAAAAAAAAAAAAAAAAAAAAAAA 18343.1 . AC=21,8,28,12,22,5;AF=0.164,0.063,0.219,0.094,0.172,0.039;AN=128;BaseQRankSum=0;ClippingRankSum=0;DP=1059;ExcessHet=4.0535;FS=0;InbreedingCoeff=0.3945;MLEAC=25,12,39,17,32,5;MLEAF=0.195,0.094,0.305,0.133,0.25,0.039;MQ=58.15;MQRankSum=0;QD=33.15;ReadPosRankSum=0.172;SOR=0.63

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin
    edited June 2018

    @prasundutta87
    Hi,

    Those are some messy indel sites. Are the final genotypes the same or different? Are all the discrepant sites messy indels?

    Thanks,
    Sheila

  • prasundutta87prasundutta87 EdinburghMember
    edited June 2018

    By 'messy' you mean multiallelic indels?

    The genotypes are different..all seem to be multiallelic indels..

    some site have ref A and alt G in one VCF file and the other has ref A and alt G,* in the other VCF file for the same position..what does that mean?

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @prasundutta87
    Hi,

    Yes, messy is many possible indel alleles as alt alleles :smile:

    Can you confirm this occurs with the very latest version? If so, I may need you to submit a bug report.

    Thanks,
    Sheila

Sign In or Register to comment.