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

SEQanswers (http://seqanswers.com/forums/index.php)

- Bioinformatics (http://seqanswers.com/forums/forumdisplay.php?f=18)

- - Extract aligned reads from a BAM file above a certain threshold (http://seqanswers.com/forums/showthread.php?t=24023)

The Snow

10-11-2012 05:52 AM

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!

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

All times are GMT -8. The time now is 08:36 PM.