Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • gatk indel calling problem

    I'm trying to get variant calling to work using a small set of simulated reads of the phi X genome. The problem I'm seeing is that indels are assigned a quality of zero and their vcf entries are incorrect.

    I used wgsim from samtools to generate simulated paired-end fastq read files of varying depths of coverage, then used bwa to align the simulated reads to the phi X reference. I followed the GATK recommendations to realign around possible indels, then to recalibrate, then I ran the UnifiedGenotyper using the -glm INDEL option along with --output-mode EMIT_ALL_SITES.

    The analysis info shows something like this:

    INFO 12:20:08,764 UnifiedGenotyper - Visited bases 5416
    INFO 12:20:08,765 UnifiedGenotyper - Callable bases 7
    INFO 12:20:08,765 UnifiedGenotyper - Confidently called bases 0
    INFO 12:20:08,765 UnifiedGenotyper - % callable bases of all loci 0.129
    INFO 12:20:08,765 UnifiedGenotyper - % confidently called bases of all loci 0.000
    INFO 12:20:08,766 UnifiedGenotyper - % confidently called bases of callable loci 0.000
    INFO 12:20:08,766 UnifiedGenotyper - Actual calls made 7
    INFO 12:20:08,766 TraversalEngine - Total runtime 1.08 secs, 0.02 min, 0.00 hours

    The corresponding vcf file shows (I added a space after the "GT:" to keep from happening):

    #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT somesample
    phix.fa 869 . T . 0 LowQual DP=183;MQ=48.18;MQ0=0;QD=0.00 GT: DP:PL ./.:183:0,0,0
    phix.fa 875 . GCT . 0 LowQual DP=190;MQ=49.29;MQ0=0;QD=0.00 GT: DP:PL ./.:190:0,0,0
    phix.fa 1408 . T . 0 LowQual DP=209;MQ=52.38;MQ0=0;QD=0.00 GT: DP:PL ./.:209:0,0,0
    phix.fa 1409 . A . 0 LowQual DP=213;MQ=52.17;MQ0=0;QD=0.00 GT: DP:PL ./.:213:0,0,0
    phix.fa 1949 . G . 0 LowQual DP=205;MQ=57.70;MQ0=0;QD=0.00 GT: DP:PL ./.:205:0,0,0
    phix.fa 1950 . A . 0 LowQual DP=205;MQ=57.82;MQ0=0;QD=0.00 GT: DP:PL ./.:205:0,0,0
    phix.fa 4225 . AC . 0 LowQual DP=179;MQ=57.51;MQ0=0;QD=0.00 GT: DP:PL ./.:179:0,0,0
    The vcf entries correspond to actual indels in the simulated data. For example, at position 870 there is an insertion of CC, and 876 there is a deletion of CT. The problem is the the "ALT" column shows a "." where it should have the actual alternate allele (TCC and G in those two cases).

    What can I do to 1) make the quality not zero; and 2) get the vcf file to correctly report the alternative allele?

  • #2
    Originally posted by Alex Renwick View Post
    I'm trying to get variant calling to work using a small set of simulated reads of the phi X genome. The problem I'm seeing is that indels are assigned a quality of zero and their vcf entries are incorrect.

    I used wgsim from samtools to generate simulated paired-end fastq read files of varying depths of coverage, then used bwa to align the simulated reads to the phi X reference. I followed the GATK recommendations to realign around possible indels, then to recalibrate, then I ran the UnifiedGenotyper using the -glm INDEL option along with --output-mode EMIT_ALL_SITES.

    The analysis info shows something like this:




    The corresponding vcf file shows (I added a space after the "GT:" to keep from happening):



    The vcf entries correspond to actual indels in the simulated data. For example, at position 870 there is an insertion of CC, and 876 there is a deletion of CT. The problem is the the "ALT" column shows a "." where it should have the actual alternate allele (TCC and G in those two cases).

    What can I do to 1) make the quality not zero; and 2) get the vcf file to correctly report the alternative allele?
    Hmmm.. you want to call all variants regardless of quality? try removing EMIT ALL SITES and lower the stand_emit_conf to 1

    Comment


    • #3
      Originally posted by NGSfan View Post
      Hmmm.. you want to call all variants regardless of quality? try removing EMIT ALL SITES and lower the stand_emit_conf to 1
      I tried that--lowering the standard to zero, in fact--and without EMIT_ALL_SITES there is no output. It's like the indels are treated as nonvariant sites.

      Comment


      • #4
        Hi guys,
        I have this exact problem. And I have tried -mmq 0 -mbq 0 -stand_emit_conf 0 -stand_call_conf 0 without seeing anything. I welcome any verification of my (non) findings with these parameters.

        I also tried using real base qualities from real data without seeing any change... still actively working on this, and I would like to hear any solution or suggestion people have.
        Last edited by oiiio; 11-18-2011, 03:33 PM.

        Comment


        • #5
          I've done some experimenting and discovered that the problem comes from the amount of sequencing error in the input data.

          The data I am using is simulated using wgsim (distributed with samtools). The default error rate is .02. With that rate, UnifiedGenotyper makes no indel calls, irrespective of setting call / emit confidences to zero and minimum base quality to zero. Snp calls work just fine, though.

          UnifiedGenotyper is happy to make indel calls when the error rate is set to zero, and it works about half the time with an error rate of 0.01.

          I haven't found any parameters that allow the user to specify sensitivity to sequencing errors. Anybody else have ideas?

          Comment


          • #6
            Nice find Alex,

            When you say it works half the time, do you mean that it is calling about half the indels in the sample? Or that half the samples you run UnifiedGenotyper on are making calls at all?

            Comment


            • #7
              I've done around a dozen trials, and it seems to be all or nothing. Either all candidate sites are called, or none of them are.

              edit: An error rate of 0.005 seems ok, too.
              Last edited by Alex Renwick; 12-02-2011, 12:28 PM.

              Comment

              Latest Articles

              Collapse

              • seqadmin
                Best Practices for Single-Cell Sequencing Analysis
                by seqadmin



                While isolating and preparing single cells for sequencing was historically the bottleneck, recent technological advancements have shifted the challenge to data analysis. This highlights the rapidly evolving nature of single-cell sequencing. The inherent complexity of single-cell analysis has intensified with the surge in data volume and the incorporation of diverse and more complex datasets. This article explores the challenges in analysis, examines common pitfalls, offers...
                06-06-2024, 07:15 AM
              • seqadmin
                Latest Developments in Precision Medicine
                by seqadmin



                Technological advances have led to drastic improvements in the field of precision medicine, enabling more personalized approaches to treatment. This article explores four leading groups that are overcoming many of the challenges of genomic profiling and precision medicine through their innovative platforms and technologies.

                Somatic Genomics
                “We have such a tremendous amount of genetic diversity that exists within each of us, and not just between us as individuals,”...
                05-24-2024, 01:16 PM

              ad_right_rmr

              Collapse

              News

              Collapse

              Topics Statistics Last Post
              Started by seqadmin, Yesterday, 07:24 AM
              0 responses
              10 views
              0 likes
              Last Post seqadmin  
              Started by seqadmin, 06-13-2024, 08:58 AM
              0 responses
              11 views
              0 likes
              Last Post seqadmin  
              Started by seqadmin, 06-12-2024, 02:20 PM
              0 responses
              16 views
              0 likes
              Last Post seqadmin  
              Started by seqadmin, 06-07-2024, 06:58 AM
              0 responses
              184 views
              0 likes
              Last Post seqadmin  
              Working...
              X