Samtools Depth - Filling in the gaps with 0 - Perl?

[2seq]

Hi everyone,
I thought that perhaps others have encountered and previously solved the following problem:

Samtools depth allows you to collect read depth information across a range in the genome using bed files to supply the coordinates.

My problem, is that in my output, the low quality regions or the elimination events appear as gaps rather than as zeros. Is there a way to change the parameters of samtools depth to report for all positions even if it just reports a 0? Or do I just need to write a script to fill in the blanks?

Current output:

19      80820312        1       0       0
19      80820313        1       0       0
19      80820314        1       0       0
19      80826337        2       0       2
19      80826338        2       0       3
19      80826339        3       0       3
19      80826340        3       0       4

Desired output:

19      80826312        1       0       0
19      80820313        1       0       0
19      80820314        1       0       0
19      80820315        0       0       0
19      80820316        0       0       0
..
19      80820335        0       0       0
19      80820336        0       0       0
19      80826337        2       0       2
19      80826338        2       0       3
19      80826339        3       0       3
19      80826340        3       0       4

**with the .. region filled in as well.

I have managed to write a code that identifies the breakpoints, but I'm not sure how to add additional lines that will run consecutively through the elimination site until it reaches the next set of depth values. Below is the code I've been working on (it's a perl one-liner and pretty raw):

perl -e '$last == 0; while (<>){@a = split(/\t/, $_, 5); $n = $a[1]; if ($last == 0) { $last = ($a[1] + 1) ;  } elsif ($n == $last){ $last = ($n + 1);  next; } elsif ($n > ($last + 1)){ print "$a[0]\t$last\t0\t0\t0\n"; $last = ($n + 1); print;}  }' depth1-2-3.bed

All help is greatly appreciated!

[WhatsOEver]

there is actually no need to write something new. You can get your desired output by using bedtools genomecov: http://bedtools.readthedocs.org/en/l...genomecov.html

In short:
using the "-d" parameter gives you a per base coverage
using the "-bga" parameters will group regions of identical coverage

Both will output bases/regions of zero coverage.

[biograd]

It doesn't look like genomecov will allow you to filter the output according to mapping quality thresholds. Anyone know what the cutoff is? It's not in the docs. thanks!

[WhatsOEver]

As far as I know, it doesn't apply a filter at all. You would need to pre-filter your samfile using your own quality threshold which is pretty straight forward using awk:

awk -F "\t" '{ if($5 > 20) print $0}' ./yourSamFile.sam

This returns only those lines where the mapping quality is greater than 20.