Extract aligned reads from a BAM file above a certain threshold

[The Snow]

Hello everybody,

I have the following problem: I have a bam file with all the aligned sequences (it was filtered with bamtools filter -isMapped true).

However I would like to filter not only the aligned reads BUT the aligned reads above a certain threshold (e.g. 200) in this way I can filter those sequences that align to a ref with just 1 or two reads

There is a way to do it?

Thanks!

[The Snow]

No one?

[dpryan]

That would be pretty easy to do in a small perl/python/whatever script. There's also a pretty easy to use C api for samtools. I imagine you could also do that in a bash script, but I would personally find it simpler to just write a quick program.

[GW_OK]

Could you filter using awk?

[dlawrence]

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

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