Hi all
I am validating indels, that I call myself against 1000 genome vcf.
I am interested by using my bam file and reference sequence in extracting regions that have 0 coverage.
Is there an easy way to do it?
•
written
by
stolarek.ir
i just a wrote a picard/java based tool to answer that question. It produces a bed file:
See
https://github.com/lindenb/jvarkit#biostar78285
$ java -jar dist/biostar78285.jar I=input.bam VALIDATION_STRINGENCY=LENIENT | head
chr1 300 330
chr1 468 571
chr1 626 726
chr1 776 816
chr1 867 990
chr1 1158 1418
chr1 1521 1538
chr1 11241 11243
chr1 11456 11470
chr1 11621 11623
if read contains deletions, does your tool falls in the same trap as bed tools, which reports deletion as 0 coverage region?
I've added an option USECIGAR to my code: if set, the program will scan the cigar string and detect the deletions in reference. (slower & requires more memory )
cool thanks. It's nice to be an inspiration for writing new tool. A bot of background behind why I am doing this. I am assessing new data type - moleculo in picking up bigger indels (currently I am at 2000 bp of good quality indels). To assess more or less how good moleculo data is I wanted to see how many of indels from 1000 genome vcf file it can pick. First result was moderate, so i went into finding why. Turned out pretty quickly that the culprit is 20% lack of coverage in my genome (plus heterozygous events, for which I need even more coverage). So I filtered out 0 coverage indels positions from 1000 genome vcf file, and then did the correct assessment
Log in to answer this question.
