I have been trying to use ReadBackedPhasing, but its is producing an output with no phasing tags (HP tag). It works with a vcf file that is been called on a single bam file, but does not produce any output from vcf files generated from vcf-subset. In short, I would like to know what INFO/FORMAT fields it looks for while assigning HP tag.
I tried to join chromosome datasets for my SNPs data separated in different chromosomes.
input command: java -Xmx160g -cp GenomeAnalysisTK.jar org.broadinstitute.gatk.tools.CatVariants -R /nobackup/data7/jose/Medicago_truncatula_genome/JCVI.Medtr.v4.20130313.fasta -V /nobackup/data7/jose/GBSv1/Msativa_diploids/hapmap/mergedSNPs/myGBSGenos_mergedSNPs_chr1_ReadDepthFiltered.vcf -V /nobackup/data7/jose/GBSv1/Msativa_diploids/hapmap/mergedSNPs/myGBSGenos_mergedSNPs_chr2_ReadDepthFiltered.vcf -V /nobackup/data7/jose/GBSv1/Msativa_diploids/hapmap/mergedSNPs/myGBSGenos_mergedSNPs_chr3_ReadDepthFiltered.vcf -V /nobackup/data7/jose/GBSv1/Msativa_diploids/hapmap/mergedSNPs/myGBSGenos_mergedSNPs_chr4_ReadDepthFiltered.vcf -V /nobackup/data7/jose/GBSv1/Msativa_diploids/hapmap/mergedSNPs/myGBSGenos_mergedSNPs_chr5_ReadDepthFiltered.vcf -V /nobackup/data7/jose/GBSv1/Msativa_diploids/hapmap/mergedSNPs/myGBSGenos_mergedSNPs_chr6_ReadDepthFiltered.vcf -V /nobackup/data7/jose/GBSv1/Msativa_diploids/hapmap/mergedSNPs/myGBSGenos_mergedSNPs_chr7_ReadDepthFiltered.vcf -V /nobackup/data7/jose/GBSv1/Msativa_diploids/hapmap/mergedSNPs/myGBSGenos_mergedSNPs_chr8_ReadDepthFiltered.vcf -out /nobackup/data7/jose/GBSv1/Msativa_diploids/hapmap/mergedSNPs/myGBSGenos_mergedSNPs_chr1-8_ReadDepthFiltered.vcf -assumeSorted
java.lang.NullPointerException at htsjdk.variant.variantcontext.writer.IndexingVariantContextWriter.setIndexSequenceDictionary(IndexingVariantContextWriter.java:178) at htsjdk.variant.variantcontext.writer.IndexingVariantContextWriter.close(IndexingVariantContextWriter.java:133) at htsjdk.variant.variantcontext.writer.VCFWriter.close(VCFWriter.java:210) at org.broadinstitute.gatk.tools.CatVariants.execute(CatVariants.java:303) at org.broadinstitute.gatk.utils.commandline.CommandLineProgram.start(CommandLineProgram.java:248) at org.broadinstitute.gatk.utils.commandline.CommandLineProgram.start(CommandLineProgram.java:155) at org.broadinstitute.gatk.tools.CatVariants.main(CatVariants.java:311)
Any help would be appreciated.
Hello! Is it possible add a tag in a vcf when I calculate the ratio of the second filed of AD and DP? for example in this case: chr13 32912299 . T C 2372.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=-2.631;ClippingRankSum=-0.154;DP=266;FS=1.546;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=-0.144;QD=8.92;ReadPosRankSum=0.162;SOR=0.596 GT:AD:DP:GQ:PL 0/1:155,111:266:99:2401,0,4319
I would like calculate 111/226 Thank you very much in advance
I've a problem with a VCF when combiningGVCFs, gatk/v3.4. "headerKey END found in VariantContext field INFO at chr2L:64003 but this key isn't defined in the VCFHeader"
-nct 30 --analysis_type HaplotypeCaller --reference_sequence ../../reference_sequences/dmel/v6.0/dm6.fa \ --input_file ../../read_mapping/reference_line/RGnb/RGnb.bam \ --variant_index_type LINEAR --variant_index_parameter 128000 -mbq 10 -L chr2L \ --out ../../read_mapping/reference_line/RGnb/RGnb.g.vcf
VCF CHROM POS ID REF ALT QUAL FILTER INFO FORMAT RGnb chr2L 43265 . A G 2147.77 . AC=2;AF=1.00;AN=2;DP=54;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=29.59;SOR=1.118 GT:AD:DP:GQ:PL 1/1:0,54:54:99:2176,162,0 chr2L 64003 . GT G 57.73 . AC=1;AF=0.500;AN=2;BaseQRankSum=0.987;ClippingRankSum=-0.189;DP=45;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=60.85;MQRankSum=0.525;QD=1.28;ReadPosRankSum=1.149;SOR=0.622 GT:AD:DP:GQ:PL 0/1:25,8:33:95:95,0,521
Combine three gvcfs
--analysis_type CombineGVCFs --reference_sequence ../../reference_sequences/dmel/v6.0/dm6.fa \ --variant rg_GVCF.list --out ./RG.g.vcf
Error message: java.lang.IllegalStateException: Key END found in VariantContext field INFO at chr2L:64003 but this key isn't defined in the VCFHeader. We require all VCFs to have complete VCF
Previous issues for "key isn't defined in the VCFHeader" http://gatkforums.broadinstitute.org/discussion/3962/genotypeandvalidate-error-key-callstatus-found-in-variantcontext-field-info#latest http://gatkforums.broadinstitute.org/discussion/2331/select-variants-says-vcf-header-is-missing https://github.com/charite/jannovar/issues/70 https://groups.google.com/forum/#!topic/gawed-awesome-analyzers/RjXtIuYQ_V8
Could anyone advise me on how to fix this?
I am trying to open small uncompressed VCFs via picard-tools-1.119, via SplitVcfs.jar.
File vcfFile = new File("test2.vcf"); VCFFileReader vcfReader = new VCFFileReader(vcfFile,false);
Even though I've specified that an index is not required in the constructor, this generates the following error:
htsjdk.tribble.TribbleException: Index not found for: MYLOCATION\test2.vcf at htsjdk.tribble.TribbleIndexedFeatureReader.query(TribbleIndexedFeatureReader.java:253) at htsjdk.variant.vcf.VCFFileReader.query(VCFFileReader.java:124) (if I don't specify the index parameter, then it generates a different error:"An index is required, but none found" )
I've attached a test file so you can see it for yourself. ( NOTE: I've had to compress it into a .zip, otherwise the forum didn't allow me to attach the file itself, but if you want to replicate my results you will need to uncompress it and try loading the plain file (test2.vcf).
My question: is there any way to open uncompressed VCF files in picard via VCFFileReader or otherwise?
PS: these are small VCFs that would be exported out from larger datasets, so compressing/indexing them would seem like an unnecessary extra step.
I'm doing a variant analysis of genomic DNA from 2 related samples. I followed the up-to-date Best practices using HaplotypeCaller in GVCF mode for both samples followed by GenotypeGVCF to compute a common vcf of variant loci. I'm looking at variants that would be sample2-specific (present in sample2 but not in sample1)
Here is a line of this file:
CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Sample1 Sample2 chrIII 91124 . A AATAAGAGGAATTAGGCT 1132.42 . AC=2;AF=0.500;AN=4;DP=47;FS=0.000;MLEAC=2;MLEAF=0.500;MQ=58.85;MQ0=0;QD=7.99 GT:AD:DP:GQ:PL 1/1:0,25:25:55:1167,55,0 0/0:.:22:33:0,33,495
In the Genotype Field, sample2.AD is a . (dot) meaning that no reads passed the Quality filters. However, sample2.DP=22 meaning that 22 reads covered this position. This line suggest that this variation is specific to sample1 (genotype HomVar 1/1) and is not present in sample2 (HomRef 0/0). But given the biological relationship between sample1 and 2 (the way they were generated), I doubt that this variation is true: it is very likely to be present in sample2 as well. It's a false
I have 416 loci like this. For the vast majority of them, sample1 and 2 likely share the same variation. But since it is not impossible that a very few of them are really sample1=HomVar and sample2=HomRef, could you suggest me a way to detect those guys? What about comparing sample1.PL(1/1) and sample2.PL(0/0) ? For example could you suggest a rule of thumb to determine their ratio ?
Hi, GATK Team
I've run into a strange case where a SNP called by HaplotypeCaller has been represented as if it were an indel:
6 16327909 . ATGCTGATGCTGC CTGCTGATGCTGC 1390.70 PASS AC=1;AC_Orig=2;AF=0.500;AF_Orig=0.040;AN=2;AN_Orig=50;BaseQRankSum=0.788;DP=10;FS=6.154;InbreedingCoeff=0.1807;MQ=59.86;MQ0=0;MQRankSum=0.406;QD=2.77;ReadPosRankSum=0.358;VQSLOD=2.78;culprit=FS GT:DP:GQ:PL 0/1:10:70:284,0,214
This VCF entry (for a single individual) comes from a multi-sample VCF that has multiple alternate "alleles" at that position:
6 16327909 . ATGCTGATGCTGC ATGC,CTGCTGATGCTGC,A 1390.70 . AC=12,3,1;AF=0.024,6.048e-03,2.016e-03;AN=496;BaseQRankSum=0.788;DP=2791;FS=6.154;InbreedingCoeff=0.1807;MLEAC=13,3,1;MLEAF=0.026,6.048e-03,2.016e-03;MQ=59.86;MQ0=0;MQRankSum=0.406;QD=2.77;ReadPosRankSum=0.358 GT:AD:DP:GQ:PL
However, this mode of representing a SNP is causing processing and analysis problems further downstream after I've split the multi-sample VCF into individual files. Is there a way to fix this problem such that variants are listed in the most parsimonious (and hopefully standard) way?
generated with gatk 2.8-1-g932cd3a
Although it is rare I see Genotype Fields that are inconsistent with the AD values (Read as table):
CHROM POS ID REF ALT FILTER QUAL ABHet ABHom AC AF AN BaseCounts BaseQRankSum DP Dels FS GC HRun HaplotypeScore LowMQ MLEAC MLEAF MQ MQ0 MQRankSum MeanDP MinDP OND PercentNBaseSolid QD ReadPosRankSum Samples Somatic VariantType cosmic.ID 1.AB 1.AD 1.DP 1.F 1.GQ 1.GT 1.MQ0 1.PL 1.Z 2.AB 2.AD 2.DP 2.F 2.GQ 2.GT 2.MQ0 2.PL 2.Z 3.AB 3.AD 3.DP 3.F 3.GQ 3.GT 3.MQ0 3.PL 3.Z 4.AB 4.AD 4.DP 4.F 4.GQ 4.GT 4.MQ0 4.PL 4.Z 5.AB 5.AD 5.DP 5.F 5.GQ 5.GT 5.MQ0 5.PL 5.Z 11 92616485 0 A C PASS 63.71 0.333 0.698 1 0.1 10 89,54,0,0 -5.631 143 0 49.552 71.29 2 4.4154 0.0000,0.0000,143 1 0.1 50.27 0 -1.645 28.6 16 0.242 0 2.36 2.125 R5_A3_1 NA SNP COSM467570 NA 24,9 33 0.2727272727 54 A/A 0 0,54,537 -1.3055824197 0.33 9,18 27 0.6666666667 96 A/C 0 96,0,178 0.8660254038 NA 21,11 32 0.34375 21 A/A 0 0,21,466 -0.8838834765 NA 12,4 16 0.25 27 A/A 0 0,27,272 -1 NA 23,12 35 0.3428571429 42 A/A 0 0,42,537 -0.9296696802
This shows that for example sample 5 has a AD value of '23,12' and a GT of 'A/A' aka homyzougous reference allele. I've included a screenshot wich shows low base quality and complete strand bias (Which I suspect to mis variants). So whats the prob? and how can i recalculate the GT's based on AD? because i cannot filter based on genotypes when they are buggy....