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.

HaplotypeCaller missing variants in highly divergent regions?

whp09whp09 EdinburghMember
edited January 2015 in Ask the GATK team

Hi all,
I am using Haplotype caller to call variants between an outgroup (my species of interest) aligned to a divergent reference of another species, then updating the reference and realigning and calling homozygous non-reference variants, until the reference converges and can be used for my species of interest. However, I am noticing that particularly divergent areas (which in my case are the most important to update!) are not being called as efficiently as other areas (see attached photo), the bottom is the bam file while the top is the called variants of the VCF file - on the left hand side of the screen, HC has done just fine, but all of a sudden, it stops making calls! My command line is the following:

HaplotypeCaller

java -jar ~/Programs/GenomeAnalysisTK.jar -T HaplotypeCaller --heterozygosity 0.2 -R ./Round2/Hmel-Hhec-round2.fa -I ERR260306_final.bam -o ERR260306_raw_calls.g.vcf --emitRefConfidence BP_RESOLUTION --variant_index_type LINEAR --variant_index_parameter 128000 -nct 8 --min_mapping_quality_score 10

Select variants from right chromosome and depth

java -jar ~/Programs/GenomeAnalysisTK.jar -T SelectVariants -R Round2/Hmel-Hhec-round2.fa -V ERR260306_raw_calls.g.vcf -L HE670299 -L HE668226 -L HE668231 -L HE668841 -L HE669356 -L HE669537 -L HE670312 -L HE670354 -L HE670553 -L HE670642 -L HE670722 -L HE670891 -L HE670902 -L HE670903 -L HE670915 -L HE670929 -L HE671059 -L HE671060 -L HE671164 -L HE671167 -L HE671168 -L HE671174 -L HE671186 -L HE671226 -L HE671241 -L HE671245 -L HE671341 -L HE671428 -L HE671449 -L HE671456 -L HE671516 -L HE671523 -L HE671529 -L HE671591 -L HE671606 -L HE671720 -L HE671721 -L HE671778 -L HE671914 -L HE671918 -L HE671949 -L HE671959 -L HE672008 -o ERR260306_raw_calls_group.g.vcf -select "DP > 5"

Genotype GVCFs

java -jar ~/Programs/GenomeAnalysisTK.jar -T GenotypeGVCFs -R ./Round2/Hmel-Hhec-round2.fa --heterozygosity 0.2 --variant ERR260306_raw_calls_group.g.vcf -o ERR260306_raw_calls_group.vcf

Select Variants SNPs but NOT Invariant sites, also AF = 1.0 for only homozygote nonref -> final

java -jar ~/Programs/GenomeAnalysisTK.jar -T SelectVariants -R ./Round2/Hmel-Hhec-round2.fa -V ERR260306_raw_calls_group.vcf -o final.vcf -selectType SNP -select 'vc.getGenotype("ERR260306").isHomVar()'

I'm sure I'm missing something, any pointers on how to make HC more sensitive?

Answers

Sign In or Register to comment.