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!

VCF output

I got VCF file for a sample (Bam file-RNAseq) using HaplotypeCaller. After extract known SNPs (with rs ID) I have the frequency of each genotype as follows:
0/0 = 6036
0/1 = 421480
1/1 = 777386
0/2 = 501
1/2 = 1571
2/2 = 113

./. = 369

Total: 1,207,456 SNPs
The questions are :
1) Why the frequency of reference homozygous (0/0) is so low (0.5%) and for the Alternative is very high (64%)?
2) I checked and compared some "Ref" and "ALT" alleles in VCF with NCBI website. They are exactly the opposite reported. Why?
3)The highest rate of callvariants achieved when DP=2. why? The range of DP in my VCF file is 1-to-446.

Thanks

Answers

  • SheilaSheila Broad InstituteMember, Broadie admin

    @mahyarhey
    Hi,

    Can you please tell us the exact command you used and which version of GATK you are using?

    Thanks,
    Sheila

  • mahyarheymahyarhey BostonMember

    Hi Sheila, I used the following commands for Calling variants using HaplotypeCaller tool:

    java -Xms4096m -Xmx4096m -jar /GenomeAnalysisTK-3.2-2/GenomeAnalysisTK.jar \
    -R ucsc.hg19.fasta \
    -T HaplotypeCaller \
    -I myfile.bam \
    -recoverDanglingHeads \
    -dontUseSoftClippedBases \
    -stand_call_conf 20.0 \
    -stand_emit_conf 20.0 \
    --min_base_quality_score 20 \
    --dbsnp dbsnp_137.hg19.vcf \
    --output_mode EMIT_VARIANTS_ONLY \
    --genotyping_mode DISCOVERY \
    --phredScaledGlobalReadMismappingRate 45 \
    --emitRefConfidence BP_RESOLUTION \
    --out output.vcf

  • SheilaSheila Broad InstituteMember, Broadie admin

    @mahyarhey

    Hi,

    1) Because you specifically selected only the variants that have been reported in dbSNP (most of which are very commonly found in people), getting a lot of 0/1 and 1/1 in that subset is expected. Also, you have have requested a GVCF, but also asked for variants only. This is incorrect and bound to mess with how hom-ref calls are handled.

    2) Where and how exactly did you look up your information, and can you post an actual example?

    3) The range of DP is not helpful to us. A histogram of coverage frequencies would be much more useful.

    Unfortunately, interpreting genotype frequencies is not something we can help you with, unless you have evidence that something is actually wrong.

    -Sheila

  • mahyarheymahyarhey BostonMember

    This is the frequency of DP across SNPs-Call:

    DP SNPs-Call

    1 400
    2 15,949
    3 12,219
    4 11,513
    5 8,994
    6 7,736
    7 6,099
    8 5,209
    9 4,379
    10 3,534
    11-20 14,338
    21-50 4,675
    51-100 449
    101-200 98
    201-300 18
    301-400 9

    > 400 4

    Total 95,623

    For GVCF, I have to use --emit-all-sites command?

  • SheilaSheila Broad InstituteMember, Broadie admin

    @mahyarhey
    Hi,

    Why do the total SNPs not match in your first post and this post?

    No, for GVCFs you do not need to use the --EMIT_ALL_SITES flag.
    http://gatkforums.broadinstitute.org/discussion/4017/what-is-a-gvcf-and-how-is-it-different-from-a-regular-vcf

    -Sheila

  • mahyarheymahyarhey BostonMember

    Hi, because the total SNPs here is just for Chromosome-1 which is the longest. For the other chromosomes I had the same pattern.The question is why with DP=1, the program is still able to call variants and I had the highest call when DP=2? it is not make sense for me!

  • SheilaSheila Broad InstituteMember, Broadie admin

    @mahyarhey
    Hi,

    The best thing to do is compare the histogram of DP against the histogram of calls. The distribution is probably very similar. The expectation is that if most sites are covered to a low depth, most variant calls will have low DP, unless there is some artifact screwing up the distributions.

    -Sheila

Sign In or Register to comment.