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.

Fake GC content?

Hello!
I tried GATK DiagnoseTargets tool to analyze my intervals metrics and found that INFO field GC sometimes differ from actual GC content from reference genome. There are more than 1000 intervals with GC=0.00. I supposed that this field is somewhat constant for the interval and is measured on reference sequence, but it looks like DiagnoseTargets tool has another logic for it. Is it possible to explain, how does DieagnoseTargets compute GC content?

I am using GATK Version=3.1-1
Tool was stared with default options:
java -jar GenomeAnalysisTK.jar -T DiagnoseTargets -R hg19.fa -I sample1.bam -I sample2.bam -I sample3.bam -I sample4.bam -I sample5.bam -I sample6.bam -I sample7.bam -L panel.bed -o output.vcf -missing missing.intervals

Example row from output vcf where GC=0 and actual GC content is not null (there are 61 GC on this interval):
chr11 75112674 . C <DT> . NO_READS END=75112787;GC=0.00;IDP=0.00

Thank you for your consideration!

Comments

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @kirill_tsukanov
    Hi,

    Can you try with the latest version and let me know if you still get the same results? I am getting an error when I try running DiagnoseTargets from version 3.1-1.

    Thanks,
    Sheila

  • Hi,
    I've just repeated experiment with Version=3.4-0, output GC Content remains the same (many intervals with zero GC).

    Thank you for reply!

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Does this only happen in intervals with NO_READS?

  • @Geraldine_VdAuwera said:
    Does this only happen in intervals with NO_READS?

    No! Look at this one:
    chr1 245017744 . A <DT> . COVERAGE_GAPS END=245017815;IDP=0.014;IGC=0.00

    IGC field is GC Content field in version 3.4-0.

    On the other hand, there are bunch of intervals with GC>0 where GC still does not match real GC content, look:
    chr2 242157644 . T <DT> . COVERAGE_GAPS END=242157806;IDP=8.33;IGC=0.227 nGCs=119 width=163 GCContent=0.7301
    There i added fields with number of GC nucleotides from reference (nGCs), width of the interval and measured GCContent, wich is much differ from GATK IGC. Look like IGC depends on quality, IDP or smth..

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @kirill_tsukanov
    Hi,

    How are you calculating IGC? I tested this on a small interval with the latest version, and IGC is calculated by the number of G and C bases in the reference divided by the number of bases in your interval. So, if the reference bases in your interval are TTG, your IGC is 0.333. (1 G divided by 3 total bases). For me, it looks like the IGC is being output properly.

    -Sheila

  • @Sheila said:
    kirill_tsukanov
    Hi,

    How are you calculating IGC? I tested this on a small interval with the latest version, and IGC is calculated by the number of G and C bases in the reference divided by the number of bases in your interval. So, if the reference bases in your interval are TTG, your IGC is 0.333. (1 G divided by 3 total bases). For me, it looks like the IGC is being output properly.

    -Sheila

    Hi,

    IGC is calculated by GATK. As you can see for the interval chr2:242157644-242157806 from my example above IGC = 0.227. In fact there are 119 G/C bases in 163 length interval, which is much greater than 0.227.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hang on, have we checked whether IGC is supposed to measure the reference GC in the interval, or the GC found in the reads spanning the interval? What does the header line in the VCF say?

  • Header line is uncertain:
    ##INFO=<ID=IGC,Number=1,Type=Float,Description="GC Content of the interval">

    I try to realise, what does IGC realy measure.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hmm, we'll have to look at the code to clarify this.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    I checked the code, it's definitely measuring the reference GC content. But it's possible that there is some logic in there that causes some positions within the interval to be ignored depending on the dataset that you run on, but it's not documented anywhere so it might not be intentional. Can you test if you get the same results for those intervals if you run on different samples? If so we would want to use that as a test case for debugging.

  • @Geraldine_VdAuwera said:
    I checked the code, it's definitely measuring the reference GC content. But it's possible that there is some logic in there that causes some positions within the interval to be ignored depending on the dataset that you run on, but it's not documented anywhere so it might not be intentional. Can you test if you get the same results for those intervals if you run on different samples? If so we would want to use that as a test case for debugging.

    Sorry for long reply!
    I've tested on separate dataset, and the problem still here. Ready to provide my dataset, command line options and list of suspicious intervals for testing, if you wish!

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin
Sign In or Register to comment.