How to pile up reads and get per base coverage?

[heso]

Hi,

I have list of sequences (17-35bases) with sequence read count in column 1.
I would like to see the end processing patterns of these sequences (to pile up similar ones and get the per-base coverage).

E.g. two sequences in my input file:

      4 TTGATCCTTCGATGTCGGCTCTTCCT
    30 ATCCTTCGATGTCGGCTCTTCCTATCATT

alignment of these sequences would look like that:

 TTGATCCTTCGATGTCGGCTCTTCCT
    ATCCTTCGATGTCGGCTCTTCCTATCAT

So the coverage of the first 3 nucleotides (TTG) is read count=4, for ATCCTTCGATGTCGGCTCTTCCT its 30+4=34 and for ATCAT it's 30


In principle, the end result/graph should look like this one here:

http://www.frontiersin.org/files/Articles/54421/fgene-04-00145-HTML/image_m/fgene-04-00145-g003.jpg

I know bedtools can output genomic coverage, but I'm afraid it won't be useful here...Does anybody have an idea how to solve this? I also have .sam files containing the same sequences, if that's of any help...


Thanks in advance

[aprice67]

Hi,

I have a script that does that from my PhD days. I actually used this for some RNA secondary structure projects, so it looks like you're working on something similar.

https://github.com/price0416/scripts/blob/master/fullbam2depth.py

Disclaimer: I wrote this like 4 years ago and don't remember exactly all the details of it, so be sure to examine the code for yourself, it's not very long.

Hope this helps!