We've moved!
This site is now read-only. You can find our new documentation site and support forum for posting questions here.
Be sure to read our welcome blog!

Missing PS field in the VCF file produced by GenotypeGVCFs

svitlanasvitlana Brussels, BelgiumMember
Hello,

I followed GATK best practices to produce a VCF file for 20 individuals. GATK version is 4.1.0.0.

The BAM files were all verified by ValidateSamFile, no errors or warnings were detected.

HaplotypeCaller was run with --emit-ref-confidence, on each file individually. All the files were then merged with CombineGVCFs and finally the resulting gVCF file was used to produce my VCF file with GenotypeGVCFs.

Ah the end, I tried to run Beagle on the VCF file but I got the following error:

IllegalArgumentException: ERROR: inconsistent number of alleles for sample SAMPLE at marker [ctg9 107774 . C T]

The VCF entry for this call is the following:

ctg9 107774 . C T 4311.43 . AC=14;AF=0.350;AN=40;BaseQRankSum=0.198;DP=295;ExcessHet=0.3091;FS=0.611;InbreedingCoeff=0.3322;MLEAC=14;MLEAF=0.350;MQ=59.85;MQRankSum=0.00;QD=27.29;ReadPosRankSum=0.361;SOR=0.618
GT:AD:DP:GQ:PGT:PID:PL:PS 1|1:0,19:19:66:1|1:107737_C_T:990,66,0:107737 0/0:18,0:18:54:.:.:0,54,691 1|1:0,18:18:54:1|1:107774_C_T:806,54,0:107774 1|1:0,13:13:39:1|1:107774_C_T:581,39,0:107774 0/0:10,0:10:30:.:.:0,30,321 0/0:10,0:10:27:.:.:0,27,405 0/0:9,0:9:27:.:.:0,27,294 0/0:18,0:18:54:.:.:0,54,704 0|1:16,2:18:36:0|1:107774_C_T:36,0,666:107774 0/0:13,0:13:39:.:.:0,39,509 0|1:5,9:14:99:0|1:107774_C_T:363,0,183:107774 0|1:19,8:27:99:0|1:107774_C_T:279,0,758:107774 0/0:12,0:12:36:.:.:0,36,501 0/0:16,0:16:45:.:.:0,45,675 0|1:7,4:11:99:0|1:107774_C_T:144,0,324:107774 0|1:6,7:13:99:0|1:107774_C_T:276,0,231:107774 0/0:13,0:13:39:.:.:0,39,497 0|1:5,3:8:99:0|1:107774_C_T:111,0,201:107774 1|1:1,16:17:12:1|1:107774_C_T:765,12,0:107774 0/0:9,0:9:27:.:.:0,27,342

The SAMPLE indicated by Beagle corresponds to the last entry (0/0:9,0:9:27:.:.:0,27,342), which doesn't correspond to the format (GT:AD:DP:GQ:PGT:PID:PL:PS), the last field is missing.

If it can help, below are shown the entries from the gVCF file for the last sample:
ctg9 107769 . T <NON_REF> . . END=107782 GT:DP:GQ:MIN_DP:PL 0/0:11:27:9:0,27,342
ctg9 107783 . A <NON_REF> . . END=107792 GT:DP:GQ:MIN_DP:PL 0/0:9:24:9:0,24,360

And here is the entry for the same variant in another sample, which has the correct number of fields in the VCF.
ctg9 107774 . C T,<NON_REF> 751.10 . BaseQRankSum=-0.182;DP=19;ExcessHet=3.0103;MLEAC=2,0;MLEAF=1.00,0.00;MQRankSum=0.000;RAW_MQandDP=68400,19;ReadPosRankSum=0.511 GT:AD:DP:GQ:PGT:PID:PL:PS:SB 1|1:1,16,0:17:12:0|1:107774_C_T:765,12,0,768,54,810:107774:1,0,5,11

I suppose that this problem is related to the physical phasing done by HaplotypeCaller.. so it could maybe be corrected by rerunning HaplotypeCaller with --do-not-run-physical-phasing option. But it will take several weeks to run on my data, so maybe there is another way to fix this problem?

Thank you in advance.

P.S. I am sorry, code formatting didn't work for me

Answers

  • DerekCADerekCA Member, Administrator, Broadie, Moderator admin

    I'm sorry @svitlana, but it looks like you asked your question at an inopportune time! We have just moved to our new site (accessible at gatk.broadinstitute.org) and have shut down posting here. That means means that you will not find the answer to your question, here!

    In the meantime, I've duplicated your question to the new forum at this url, where you can add your own comments and follow the thread for updates. Please be sure to register for a new forum account (this new account will extend to all the software we support, including Terra and Cromwell/WDL). We apologize for any inconvenience!

    You can find more information about the move on our blog post.

Sign In or Register to comment.