Heads up:
We’re moving the GATK website, docs and forum to a new platform. Read the full story and breakdown of key changes on this blog.
Notice:
If you happen to see a question you know the answer to, please do chime in and help your fellow community members. We encourage our fourm members to be more involved, jump in and help out your fellow researchers with their questions. GATK forum is a community forum and helping each other with using GATK tools and research is the cornerstone of our success as a genomics research community.We appreciate your help!

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.
Attention:
We will be out of the office for a Broad Institute event from Dec 10th to Dec 11th 2019. We will be back to monitor the GATK forum on Dec 12th 2019. In the meantime we encourage you to help out other community members with their queries.
Thank you for your patience!

Long tips in phylogeny after GATK best practices for variant calling

I have a population of 60 relatively closely related samples (non-human). After running through the GATK best practices to identify variants (with hard-filtering instead of VQSR), the phylogenies I reconstruct with the SNP VCF file have huge tips relative to internal branches, including the same sample that was sequenced twice. 99% of the variation between A and B is the same variation that exists between A and any other sample. For example, say Sample01A and Sample01B are differentiated by 1000 SNPs, while samples 01A and 60 are differentiated by 1050 SNPs, where 01A and 01B are two preps of the same sample.

Is this a symptom of a typical mis-use of GATKs pipeline? Or is this more likely occurring because of my samples specifically, such as my samples are too closely related and sequencing error is drowning out real variation?

Thanks for any suggestions, if possible,

Alex

Best Answer

Answers

  • aberry814aberry814 Member

    Hi Redmar, thanks for the response. I actually have hundreds of SNPs between the same isolate, however the genome is 40Mb and diploid, so it may be a similar percentage of base pairs as your Salmonella isolates. That's a great suggestion, though. Just by quickly looking at a couple of the variants, most do seem to be what I consider "sketchy" calls. That is, making a confident call that it is homozygous when 20% of the reads call it ALT, and vice versa. Your 90% cutoff seems appropriate. I'll read through the JEXL documentation to figure out how to set up that AD filter and report back.

  • aberry814aberry814 Member

    I tried the following command to extract AD/DP, copied verbatim from the JEXL documentation, and it yielded a blank VCF file:

    java -jar GenomeAnalysisTK.jar -T SelectVariants -R reference.fasta -V 2016.08.08_TC02_only.vcf -select "vc.getGenotype("TC002").getAD().0 / vc.getGenotype("TC002").getDP() > 0.50" -o test02.vcf

    I tested many variations. Single quotes, for example, gave me a "divide by" error.

    I've used variant context commands before, but for some reason I can't get this one to cooperate. Any suggestions?

  • aberry814aberry814 Member

    In case anybody has the "Divide by" error, try this:

    First, remove all sites where depth is zero:

    java -jar GenomeAnalysisTK.jar -T SelectVariants -R reference.fasta -V input.vcf -select 'vc.getGenotype("Sample01").getDP() > 0 && vc.getGenotype("Sample02").getDP() > 0' -o DPnonzero.vcf

    Then, in order to keep only heterozygous calls in which at least 20% of the total calls were "minor allele" calls, I needed to multiply the .getAD() and .getDP() values by 1.00 (1.0 didn't work, and not multiplying didn't work):

    java -jar GenomeAnalysisTK.jar -T SelectVariants -R /reference.fasta -V input.vcf -select '(vc.getGenotype("Sample01").isHet() && ((1.00vc.getGenotype("Sample01").getAD().0) / (1.00vc.getGenotype("Sample01").getDP()) < 0.80) && ((1.00vc.getGenotype("Sample01").getAD().0) / (1.00vc.getGenotype("Sample01").getDP()) > 0.20))' -o output.vcf

Sign In or Register to comment.