Question: For a given set of coordinates is there a simple fast way of determining the 'summit' coordinate,
0
gravatar for User 9088
4.9 years ago by
User 908820
United States
User 908820 wrote:

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.

dnase-seq • 1.3k views
ADD COMMENTlink modified 4.5 years ago by Biostar ♦♦ 20 • written 4.9 years ago by User 908820
0
gravatar for Devon Ryan
4.9 years ago by
Devon Ryan89k
Freiburg, Germany
Devon Ryan89k wrote:

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}'

ADD COMMENTlink modified 4.9 years ago • written 4.9 years ago by Devon Ryan89k
0
gravatar for Alex Reynolds
4.9 years ago by
Alex Reynolds28k
Seattle, WA USA
Alex Reynolds28k wrote:

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.

ADD COMMENTlink modified 4.9 years ago • written 4.9 years ago by Alex Reynolds28k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1828 users visited in the last hour