Calculating max coverage per gene from BAM + GFF

[nk]

I have a BAM file with alignments and a GFF file with non-overlapping gene annotation. Using these, I would like to find out what the highest coverage for each gene was.

In other words, I am looking for something like HTSeq-count (http://www-huber.embl.de/users/ander...doc/count.html), except I want the max coverage of reads to be returned instead of the total count of reads.

Is there any tool that can do something like this or am I stuck writing my own script?

[dpryan]

Do you mean that you want the maximum per-base coverage in a given gene? I'm not familiar with any tool specifically for that, but the simplest route to making one would probably involve parsing the output from "samtools mpileup", which can be told to only output regions of interest.

[nk]

Yes, exactly. Thanks for the suggestion to use mpileup - isn't that usually just used for SNP calling though?

[dpryan]

Yeah, but it also gives the per-base depth (the rest can just be ignored). The benefit of this is that it makes filtering by MAPQ and base phred score easy simple, since samtools will do that part for you. Just make sure to adjust the -d parameter to a large value!

[nk]

In case anybody ever stumbles upon this - my solution ended up being to use bedtools genomecov to calculate the per-base genome coverage, bedtools intersect to overlap the genome coverage with my gff file and then awk to find the max coverage per feature.