SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
error during GATK indel realigner david.tamborero Bioinformatics 6 07-18-2012 05:30 AM
Dindel last stage error nothing in the var.VCF file pankajsvats Bioinformatics 3 02-02-2012 05:00 AM
Dindel: Accurate indel calls from short-read data lh3 Literature Watch 0 10-29-2010 06:50 AM
MIRA output for Illumina mapping giving 100% coverage! Kasycas Bioinformatics 1 09-10-2010 03:45 AM
PubMed: Preparing a re-sequencing DNA library of 2 cancer candidate genes using the l Newsbot! Literature Watch 0 05-28-2009 05:00 AM

Reply
 
Thread Tools
Old 12-19-2010, 09:04 AM   #1
gaffa
Member
 
Location: Gothenburg/Uppsala, Sweden

Join Date: Oct 2010
Posts: 82
Default Dindel giving error for every candidate indel

I'm running Dindel 1.01 on bam files produced by BWA and Stampy, looking for small indels. All the steps run to completion but every candidate indel encounters an error, so the final vcf file is completely empty. The errors are the following (output from cat *.glf.txt | grep "chr" | cut -d " " -f 1 | sort | uniq -c):

Code:
     14 error_above_read_count_threshold
      3 error_Cannot_find_reference_sequence.
    598 error_too_few_reads
These proportions are typical of all runs I do - "error_too_few_reads" is dominating, though I highly doubt that there are actually too few reads - my coverage is ~25x, and through samtools pileup and manual inspection I know that there are plenty of indels that are well covered.

When I run the "getCIGARindels" step, I get some alarming messages printed to stdout:
Code:
Parsing indels from CIGAR strings...
Wrote indels in CIGARS for target chr1 to file candidate_dindels
Wrote indels in CIGARS for target chr2 to file candidate_dindels
error: faidx error, len==0
start: -24 end: 
error: faidx error, len==0
start: -11 end:
I get about ~20 of these "faidx error" in a typical run (though not for every chromosome). I can't find anything wrong with my reference genome file or it's index file (indexed using samtools faidx - appended at end of post). I don't know if these errors are related to the above ones.

My commands (running on 64bit Ubuntu):
Code:
dindel-1.01-linux-64bit --analysis getCIGARindels --bamFile input.bam --outputFile candidate_dindels --ref ref_genome.fasta

python makeWindows.py --inputVarFile candidate_dindels.variants.txt --windowFilePrefix realign_windows --numWindowsPerFile 2000

dindel-1.01-linux-64bit --analysis indels --bamFile input.bam --ref ref_genome.fasta --libFile candidate_dindels.libraries.txt --varFile $infile --outputFile $outfile
My fasta.fai file:
Code:
chr1	230208	6	100	101
chr2	813178	232523	100	101
chr3	316617	1053839	100	101
chr4	1531919	1373629	100	101
chr5	576869	2920874	100	101
chr6	270148	3503518	100	101
chr7	1090947	3776374	100	101
chr8	562643	4878237	100	101
chr9	439885	5446513	100	101
chr10	745741	5890804	100	101
chr11	666454	6644010	100	101
chr12	1078175	7317136	100	101
chr13	924429	8406100	100	101
chr14	784334	9339781	100	101
chr15	1091289	10131966	100	101
chr16	948062	11234175	100	101
chr17	85779	12191725	100	101
chr18	6318	12278369	100	101
gaffa is offline   Reply With Quote
Old 01-11-2011, 07:31 AM   #2
druano
Junior Member
 
Location: Netherlands

Join Date: Mar 2010
Posts: 1
Default

Hi,
I have exactly the same problem.
Did you ever figure out what the problem was?
druano is offline   Reply With Quote
Old 01-11-2011, 01:08 PM   #3
gaffa
Member
 
Location: Gothenburg/Uppsala, Sweden

Join Date: Oct 2010
Posts: 82
Default

Yes in fact - Kees Albers has confirmed that this is caused by the omission of the --doDiploid (or --doPooled) flag. I left the --doDiploid flag out because I ran Dindel on haploid samples, but that will not work. You need one of the two "--do*" flags. Might it be your case too that you are not specifying this flag?

Also the faidx errors are apparently nothing to worry about.
gaffa is offline   Reply With Quote
Old 03-09-2011, 06:02 AM   #4
colindaven
Senior Member
 
Location: Germany

Join Date: Oct 2008
Posts: 415
Default

Interesting. Thanks for that. Now why is it Dindel doesn't seem to have a mode for haploid samples ?
colindaven is offline   Reply With Quote
Old 03-14-2011, 11:18 AM   #5
gaffa
Member
 
Location: Gothenburg/Uppsala, Sweden

Join Date: Oct 2010
Posts: 82
Default

Quote:
Originally Posted by colindaven View Post
Interesting. Thanks for that. Now why is it Dindel doesn't seem to have a mode for haploid samples ?
Well, it is common for software in the NGS field to not support haploid samples. In earlier versions of Dindel there was a "force homozygous" option, but not anymore.

However, what I do is to run Dindel in pooled mode, and then use the "makeGenotypeLikelihoodFilePooled.py" script that comes with the program. This script prints a file with the likelihoods for homozygous ref/ref, heterozygous ref/alt and homozygous ref/ref - grabbing only the two homozygous likelihoods (ignoring the heterozygous likelihood) from this file allows me to force homozygosity (though then I have to compute some measure of confidence myself, e.g. a likelihood ratio or a posterior probability).
gaffa is offline   Reply With Quote
Reply

Tags
dindel, indel, stampy

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 09:41 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO