Extract aligned reads from a BAM file above a certain threshold

The Snow · 10-11-2012, 05:52 AM

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 · 10-17-2012, 05:45 AM

No one?

dpryan · 10-17-2012, 09:27 AM

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 · 10-17-2012, 10:37 AM

Could you filter using awk?

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

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
extract the positions from genome where reads from all the samples have aligned	Maulik23	Bioinformatics	1	09-11-2012 11:45 AM
Extract aligned sequence coordinates from SAM or BAM file	pirates.genome	Bioinformatics	5	08-20-2012 09:06 AM
Extract base from bam file	empyrean	Bioinformatics	4	07-03-2012 04:47 AM
Fastest way to extract differing positions from each alignment in a BAM file	CHRYSES	Bioinformatics	5	12-14-2011 12:28 PM
Extract perfectly mapped reads from SAM/BAM file	Graham Etherington	Bioinformatics	2	07-21-2011 08:27 AM

10-11-2012, 05:52 AM	#1
The Snow Junior Member Location: Italy Join Date: Jun 2012 Posts: 7	Extract aligned reads from a BAM file above a certain threshold 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! Last edited by The Snow; 10-11-2012 at 07:41 AM.

10-17-2012, 05:45 AM	#2
The Snow Junior Member Location: Italy Join Date: Jun 2012 Posts: 7	No one?

10-17-2012, 09:27 AM	#3
dpryan Devon Ryan Location: Freiburg, Germany Join Date: Jul 2011 Posts: 3,480	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.

10-17-2012, 10:37 AM	#4
GW_OK Senior Member Location: Oklahoma Join Date: Sep 2009 Posts: 411	Could you filter using awk?