Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • Alex Renwick
    replied
    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.

    Leave a comment:


  • oiiio
    replied
    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?

    Leave a comment:


  • Alex Renwick
    replied
    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?

    Leave a comment:


  • oiiio
    replied
    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.

    Leave a comment:


  • Alex Renwick
    replied
    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.

    Leave a comment:


  • NGSfan
    replied
    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

    Leave a comment:


  • Alex Renwick
    started a topic gatk indel calling problem

    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?

Latest Articles

Collapse

  • seqadmin
    Exploring the Dynamics of the Tumor Microenvironment
    by seqadmin




    The complexity of cancer is clearly demonstrated in the diverse ecosystem of the tumor microenvironment (TME). The TME is made up of numerous cell types and its development begins with the changes that happen during oncogenesis. “Genomic mutations, copy number changes, epigenetic alterations, and alternative gene expression occur to varying degrees within the affected tumor cells,” explained Andrea O’Hara, Ph.D., Strategic Technical Specialist at Azenta. “As...
    07-08-2024, 03:19 PM
  • seqadmin
    Exploring Human Diversity Through Large-Scale Omics
    by seqadmin


    In 2003, researchers from the Human Genome Project (HGP) announced the most comprehensive genome to date1. Although the genome wasn’t fully completed until nearly 20 years later2, numerous large-scale projects, such as the International HapMap Project and 1000 Genomes Project, continued the HGP's work, capturing extensive variation and genomic diversity within humans. Recently, newer initiatives have significantly increased in scale and expanded beyond genomics, offering a more detailed...
    06-25-2024, 06:43 AM

ad_right_rmr

Collapse

News

Collapse

Topics Statistics Last Post
Started by seqadmin, Today, 11:09 AM
0 responses
15 views
0 likes
Last Post seqadmin  
Started by seqadmin, 07-19-2024, 07:20 AM
0 responses
147 views
0 likes
Last Post seqadmin  
Started by seqadmin, 07-16-2024, 05:49 AM
0 responses
121 views
0 likes
Last Post seqadmin  
Started by seqadmin, 07-15-2024, 06:53 AM
0 responses
111 views
0 likes
Last Post seqadmin  
Working...
X