SEQanswers

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 12:48 PM.

Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.