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.
2 answers
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}'
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.
Log in to answer this question.