bam-readcount fails to report all positions

Question

[mark.rose]

Hello

I'm using bam-readcount (version 0.7.4-unstable-35-0f3b356-dirty). Though it's version name doesn't exactly instill confidence, it is apparently the latest version available on github. Anyway, my problem is that it is failing to report on a large swath of my 17738 bp reference sequence. As shown below, bam-readcount output skips from position 8062 to 13923.

my_ref      8062    g       2       =:0:0.00:0.00:0.00:0:0:0.00:0.00
:0.00:0:0.00:0.00:0.00  A:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00
C:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00  G:2:1.00:33.00:1.00:1:1:
0.04:0.01:35.00:1:0.01:75.00:0.50       T:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:
0.00:0.00:0.00  N:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00

my_ref      13923   c       1       =:0:0.00:0.00:0.00:0:0:0.00:0.00
:0.00:0:0.00:0.00:0.00  A:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00
C:1:40.00:33.00:40.00:1:0:0.00:0.00:0.00:1:0.98:126.00:0.98     G:0:0.00:0.00:0.
00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00  T:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:
0.00:0.00:0.00  N:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00
enter code here

I have checked the bam file used as input to bam-readcount and have found that, while there are some gaps in covearage in this unrepresented interval, most of it is covered at about as well as regions that are represented in the bam-readcount output.

I've used three commands (shown below) but all give the same result.

bam-readcount  -f  my_ref.fa sample1-gsnap.bam.sort.bam my_ref
bam-readcount -b 0 -q 0  -f  my_ref.fa sample1-gsnap.bam.sort.bam my_ref
bam-readcount  -i -f  my_ref.fa sample1-gsnap.bam.sort.bam my_ref

What am I doing wrong?

Thanks for your help

Answer 1

[mark.rose]

Now, interestingly, when I described the problem above, I had used gsnap to align the reads. I just tried bowtie2. Now the gaps in the bam-readcount output more closely align with releative gaps in coverage. I say relative because there are regions inside these gaps which still contain coverage, albeit much reduced. Here is a screen shot of IGV showing the bam-readcouunt gap from 14955 to 15523.

Answer 2

[mark.rose]

Well obviously the post of the image didn't work. What it showed was a region with approx. the above coordinates which was relatively free of alignments except for about 10 at position 15112.

Answer 3

[ernfrid]

1) If you'd prefer a stable release over the most recent code changes, you can checkout the v0.7.4 tag (git checkout v0.7.4) before following the build instructions.

2) Your command lines look correct to me. It's interesting that you see improvement when using bowtie. I have never used gsnap. Off the top of my head, I don't have any ideas as to what the issue might be. If you could share the or some of the problematic regions that might provide some insight.