SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



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 10:45 AM
Extract aligned sequence coordinates from SAM or BAM file pirates.genome Bioinformatics 5 08-20-2012 08:06 AM
Extract base from bam file empyrean Bioinformatics 4 07-03-2012 03:47 AM
Fastest way to extract differing positions from each alignment in a BAM file CHRYSES Bioinformatics 5 12-14-2011 11:28 AM
Extract perfectly mapped reads from SAM/BAM file Graham Etherington Bioinformatics 2 07-21-2011 07:27 AM

Reply
 
Thread Tools
Old 10-11-2012, 04:52 AM   #1
The Snow
Junior Member
 
Location: Italy

Join Date: Jun 2012
Posts: 7
Default 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 06:41 AM.
The Snow is offline   Reply With Quote
Old 10-17-2012, 04:45 AM   #2
The Snow
Junior Member
 
Location: Italy

Join Date: Jun 2012
Posts: 7
Default

No one?
The Snow is offline   Reply With Quote
Old 10-17-2012, 08:27 AM   #3
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

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.
dpryan is offline   Reply With Quote
Old 10-17-2012, 09:37 AM   #4
GW_OK
Senior Member
 
Location: Oklahoma

Join Date: Sep 2009
Posts: 411
Default

Could you filter using awk?
GW_OK is offline   Reply With Quote
Old 07-29-2013, 02:02 AM   #5
dlawrence
Junior Member
 
Location: South Australia

Join Date: Mar 2012
Posts: 8
Default

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
dlawrence is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 04:10 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO