For a given set of coordinates is there a simple fast way of determining the 'summit' coordinate,

Question

[User 9088]

I have a BAM file from DNase-seq experiment and a result of peak calling analysis.

For a given set of coordinates is there a simple fast way of determining the 'summit' coordinate, i.e the coordinate with the highest read pileup.

The peak caller does report the summit from its own data, but i want to check out different regions.

Thanks.

Answer 1

[Devon Ryan]

The depth is the 4th column, though that may not be correct if your reads are paired-end. You could find the location of the maximum with a simple awk script:

samtools mpileup ...options... foo.bam | awk '{if($4 > maxdepth) {chr=$1; pos=$2; maxdepth=$4;}}END{print chr,pos,maxdepth}'

Answer 2

[Alex Reynolds]

One option is to use BEDOPS and format your regions-of-interest into a BED file called regions.bed:

$ bam2bed < foo.bam | bedmap --echo --max-element regions.bed - > answer.bed

The bam2bed script converts the BAM file to BED and pipes it to bedmap, which uses the score operator --max-element to return the maximum scoring BAM element over each region in regions.bed.