SEQanswers - View Single Post - Extract aligned reads from a BAM file above a certain threshold

dlawrence · 07-29-2013, 03:02 AM

Install bedtools and then a small script like the below should do it:

Code:

BAM=foo.bam
MIN_DEPTH=10

genomeCoverageBed -ibam ${BAM} -bg | awk -v "min_depth=${MIN_DEPTH}" '{split($0, cols,"\t"); if (int(cols[4]) > min_depth) { print $0; }}' > min_depth_regions.bedgraph

intersectBed -abam ${BAM} -b min_depth_regions.bedgraph > minimum_depth.bam