SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Subsampling using 'head -n #"? kga1978 Bioinformatics 32 05-04-2013 12:22 PM
Subsampling Fastq chrisbala Bioinformatics 9 05-02-2013 08:16 AM
0 depth in pileup files from BAM files ? jjlaisnoopy Bioinformatics 3 04-22-2012 07:38 AM
Reverse engineering BAM files: BAM -> FASTQ gene coder Bioinformatics 3 01-03-2012 02:42 PM
NEw to Chip-seq and have .bam/.sam/.bam.bai files... then what? NGS newbie Bioinformatics 11 05-25-2011 07:48 AM

Reply
 
Thread Tools
Old 08-24-2012, 03:04 AM   #1
cdias
Junior Member
 
Location: Canada

Join Date: Jun 2012
Posts: 8
Default Subsampling bam files

Hello,

I need to extract reads from a bam file - for each genomic position, I need to extract a different number of reads.

Eg.
CHR START END #READS
2 12345 12350 10
2 14567 14568 7

I have gone through samtools, and cannot see a way to sample a specified number of reads for each position. I have also ran into some elegant perl scripts that will do this by position and %, but not variable for each position.

Can anyone help?
Thanks in advance!
cdias is offline   Reply With Quote
Old 08-24-2012, 03:27 AM   #2
TiborNagy
Senior Member
 
Location: Budapest

Join Date: Mar 2010
Posts: 329
Default

First of all sort and index your bam file.
Then you can use a shell script like this:

samtools view input.bam 2:12345-12350 | head 10

The head command will get the desired number of reads.
TiborNagy is offline   Reply With Quote
Old 08-24-2012, 09:04 AM   #3
cdias
Junior Member
 
Location: Canada

Join Date: Jun 2012
Posts: 8
Smile header works

Thanks! For what I needed worked like a charm!

The one problem with this command, is that it includes all heads -> you have to use head n+y
This y is the number of lines before the actual reads (^@), and is specific to each bam file.
It works great for a couple of files (my case), but will not work if you want to automate it for a larger number of files (unless you code that in as well).

cheers
cdias is offline   Reply With Quote
Old 08-27-2012, 01:14 AM   #4
TiborNagy
Senior Member
 
Location: Budapest

Join Date: Mar 2010
Posts: 329
Default

It is strange. Samtools view show headers only, when you use -h parameter.
Or you can use grep -v "^@" to remove header.

You can use it in a bash scripts.
TiborNagy is offline   Reply With Quote
Reply

Tags
bam file, bed, extract, read count, samtools

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 03:21 AM.


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