Question: bam-readcount fails to report all positions
0
gravatar for mark.rose
3.2 years ago by
mark.rose30
United States
mark.rose30 wrote:

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

ADD COMMENTlink modified 3.2 years ago by ernfrid200 • written 3.2 years ago by mark.rose30
0
gravatar for mark.rose
3.2 years ago by
mark.rose30
United States
mark.rose30 wrote:

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.

enter image description here

ADD COMMENTlink written 3.2 years ago by mark.rose30
0
gravatar for mark.rose
3.2 years ago by
mark.rose30
United States
mark.rose30 wrote:

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.

ADD COMMENTlink written 3.2 years ago by mark.rose30
0
gravatar for ernfrid
3.2 years ago by
ernfrid200
United States
ernfrid200 wrote:

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.

ADD COMMENTlink written 3.2 years ago by ernfrid200
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1697 users visited in the last hour