Calculating average coverage wih Pysam

[Martino]

Hi,

I have to program a program, which calculates each average at each position.

The idea was that that program was good, but it is too slow. That is why the code must be adapted to get is faster.

Does someone knows code for that?

[arvid]

I'm not sure what you mean by average, if you talk about each position; do you mean each precise position or some kind of window?

Have a look at "bedtools coverage" and "bedtools genomecov"; I like to use

bedtools genomecov -bg -split -ibam file.bam > file_coverage.bedgraph

which gives me a BedGraph; if you use "-d" instead of "-bg" you get the read count at each and every position.

You can get BEDTools here: http://code.google.com/p/bedtools/
There is also a python interface for it (pybedtools), but I'm not sure whether it makes sense to use it in this case, haven't really tried it.