Subsampling bam files

[cdias]

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!

[TiborNagy]

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.

[cdias]

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

[TiborNagy]

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.