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.

Is it possible to create a VCF file with information from a specific set of sites?

RenanRenan São Paulo, BrazilMember
Hi,

Basically, I have a set of fastq files from ~50 samples and I need a VCF file containing information about ~1.1 million of markers of all samples.
I prepared all 50 filtered bam files using GATK best practices and now I'm trying to create the raw VCF file just for a subset of positions of interest. The great problem is that the fastq sequences were obtained from a capture kit with low coverage and the sites absent in the final VCF file can be non-variant or missing data.
I tried two different approaches using HaplotypeCaller:
(1) Using a VCF file containing the positions of interest in other samples:

gatk HaplotypeCaller \
--java-options '-Xmx32g -XX:ParallelGCThreads=1' \
-R ref.fa.gz \
-I sample1.sorted.dedup.bam \
-I sample2.sorted.dedup.bam \
...
-I sampleN.sorted.dedup.bam \
-O output.raw.vcf \
-L other_samples.vcf \
--output-mode EMIT_ALL_SITES

(2) Using a ".interval_list" file with the specific site positions:

gatk HaplotypeCaller \
--java-options '-Xmx32g -XX:ParallelGCThreads=1' \
-R ref.fa.gz \
-I sample1.sorted.dedup.bam \
-I sample2.sorted.dedup.bam \
...
-I sampleN.sorted.dedup.bam \
-O output.raw.vcf \
-L list_of_snp_positions.interval_list \
--output-mode EMIT_ALL_SITES

As a result I obtained two VCF files containing information from the same 576,468 sites (approximately half of the information I desired).

I'm using java version "1.8.0_201" and The Genome Analysis Toolkit (GATK) v4.1.0.0

Then, my question is: Is it possible to create a VCF file with information from a specific set of sites?

Thanks in advance!

Answers

  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    Hi @Renan

    I am not sure I understand your question.Do you want to know how to use HaplotypeCaller to call variants a on a specific set of sites? If that is your questions, then the answer is option1 you mentioned above. Take a look at this doc: https://software.broadinstitute.org/gatk/documentation/article?id=11009

  • RenanRenan São Paulo, BrazilMember
    Hi @bhanuGandham

    Thanks for your answer, but it is not exactly what I want to know. Probably my question was not clear, please let me explain it better.
    What I want as an outcome is a VCF file containing genotype information of a list of ~1.1 million of sites for all my 50 samples.
    I already executed both commands I wrote above and they produced identical VCF files with just 576,468 sites. I still need the information from the rest ~500.000 sites, even they being non-variant or missing.
    Both lists I gave to the software with "-L" argument in the command lines above have the complete list of ~1.1 million of markers; even so, my output VCF has only half of the information.
    How can I get the information from the sites that were not in the output VCF?

    Thanks in advance!
  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    Hi @Renan

    Ah I see. Sorry I misunderstood your question.

    So the argument --output-mode EMIT_ALL_SITES is a little misleading. It does not email all variant and non-variant sites. It is used to emit all variant sites irrespective of their quality.

    In the latest release , we tried to mitigate this ambiguity by using the argument `--output-mode EMIT_ALL_ACTIVE_SITES' instead. However, we now realize that this might still be a little ambiguous and so we are looking into fixing this. You can see the progress of this fix here:

    Solution: What you are looking for is a GVCF. GVCF has records for all sites, whether there is a variant call there or not. Take a look at this doc: https://software.broadinstitute.org/gatk/documentation/article?id=11004

  • RenanRenan São Paulo, BrazilMember
    Hi @bhanuGandham

    I already tried to generate per chromosome GVCFs and then extract positions of interest, but GVCF files are too big and take too long to be created. Since the information I need is just a tiny part of the GVCF file, I had to look for an alternative way to reach my goal.
    If you allow me, I would suggest implementing a way to create VCFs/GVCFs from a file containing a predefined list of positions (variant or not).

    I really appreciate your attention.
    Thanks a lot.
  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    Hi @Renan

    The way to create VCFs from a file containing a predefined list of positions is using -L argument with a vcf file. See section D in https://software.broadinstitute.org/gatk/documentation/article?id=11009. If you are looking for non-variant sites then as mentioned earlier, you are looking for a gvcf file. I hope this helps bring some clarity.

Sign In or Register to comment.