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 on defined set of SNPs

meduamedua InnsbruckMember

I have a set of sequenced PCR amplicons in bam files (one per samples) where I need to call defined SNP positions to get a combined genotype, no matter if it is reference or alternative homozygous or heterozygous.

I am running the following commands:

java -jar GATK/GenomeAnalysisTK.jar \
-T HaplotypeCaller \
-R $1 \
-I $f \
-o $out.vcf \
-dontUseSoftClippedBases \
--dbsnp workspace/vcforiginal.vcf \
-ERC BP_RESOLUTION \
-L projects/ION-Analysis-files/TSAncestryPanel/PrecisionID_AncestryPanel_hotspots.bed \
--downsampling_type none \
--genotyping_mode DISCOVERY \

Here my three questions:

1) Any idea why it is unable to update the index for the .vcf file ? Generally, why would it want to do that ?

Done. There were 3 WARN messages, the first 3 are repeated below.
WARN 11:32:04,723 RMDTrackBuilder - Unable to update index with the sequence dictionary for file /home/GMI/q023me/workspace/vcforiginal.vcf.idx; this will not affect your run of the GATK
WARN 11:32:04,946 InbreedingCoeff - Annotation will not be calculated. InbreedingCoeff requires at least 10 unrelated samples.
WARN 11:32:05,440 HaplotypeScore - Annotation will not be calculated, must be called from UnifiedGenotyper, not org.broadinstitute.gatk.tools.walke

2)
The results file shows the dbSNP annotation only for heterozygous calls, any way I can change this to show it for all positions ?

3)
Would -ERC BP_RESOLUTION be the most appropriate option ? Why can I not use genotype_given_allels option with this ? (I know I'm missing a point here)

Thank you in advance and kind regards,

Mayra

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    1) It might be a permissions issue, for example if you do not have write permissions to the directory where the file lives.

    2) It should show also for homozygous variants, but it will not show for homozygous reference. It should be possible to add this as post-processing of your final VCF (after GenotypeGVCFs and filtering) using VariantAnnotator.

    3) Yes, -ERC BP_RESOLUTION is great for this use case. Genotyping given alleles is an old functionality that was designed for UnifiedGenotyper; with HaplotypeCaller it has some problems so I wouldn't recommend it.

  • meduamedua InnsbruckMember

    Thank you for the quick answer, Geraldine!
    I have attached a test output file. Some 1/1 variants are not annotated, some are, don't quite get why ?
    I also don't quite understand what GenotypeGVCFs does, I ran it but since I don't want to combine samples, not sure what is the purpose otherwise ?
    I am trying to get all dbSNP rsno. annotation and if possible BaseCountsBySample and StrandAlleleCountsBySample for all positions,(including the reference genotype posititons). VariantAnnotator on it's own does not do those annotations for neither, it's telling me, so I called them from Haplotypecaller, but then only get them for homo and heterozygous variants.
    I m just wondering if it is even possible to do that, otherwise I'll do it manually or write a script.

    thanks again for the quick help and kind regards,

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    Hi @medua,

    Can you help us dig into this by posting example records where you expect a dbSNP annotation but the site does not have one?

    Thanks,
    Soo Hee

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @medua
    Hi,

    Have a look at this article for why you should analyze your samples together.

    -Sheila

  • meduamedua InnsbruckMember

    Hi guys, thanks for the link and I will re-analyse and send examples once I get some time. I had to finish analysis on the samples for now.

Sign In or Register to comment.