Extract reads from multiple regions from bam file

[meher]

HI,

samtools view aln.sorted.bam chr2:20,100,000-20,200,000 is used to extract reads from specific regions.

Is there a way to extract if we have multiple regions specified in a bed file format?


Regards
Mehar

[krobison]

You can continue specifying regions after the first one

samtools view aln.sorted.bam chr2:20,100,000-20,200,000 chr2:30,100,000-30,200,000

(BTW, will it really tolerate the commas?)

It's pretty easy to convert BED format to this using your favorite text mangling language; mine is Perl:

samtools view aln.sorted.bam `perl -p -e 'chomp; @a=split(/\t/); print "$a[0]:$a[1]-$a[2]"' sleepy.bed`

[meher]

You can continue specifying regions after the first one

samtools view aln.sorted.bam chr2:20,100,000-20,200,000 chr2:30,100,000-30,200,000

(BTW, will it really tolerate the commas?)

It's pretty easy to convert BED format to this using your favorite text mangling language; mine is Perl:

samtools view aln.sorted.bam `perl -p -e 'chomp; @a=split(/\t/); print "$a[0]:$a[1]-$a[2]"' sleepy.bed`

Thanks for the answer.

I have tried it but, the samtools is throwing an error "Argument list too long".
Any other way of doing this?

[dariober]

It's pretty easy to convert BED format to this using your favorite text mangling language; mine is Perl:

samtools view aln.sorted.bam `perl -p -e 'chomp; @a=split(/\t/); print "$a[0]:$a[1]-$a[2]"' sleepy.bed`

Hi- You don't need to convert BED format to string. samtools view accepts a bed file directly (see -L option "output alignments overlapping the input BED FILE")

Dario

[meher]

Extract reads from multiple regions from bam file

Hi- You don't need to convert BED format to string. samtools view accepts a bed file directly (see -L option "output alignments overlapping the input BED FILE")

Dario

Hi,

I tried to do this using the following command

samtools view -L sample.bed test.bam > test1.bam

and surprisingly the output files is thrice in size than the original and i tried to count the number of reads on the output bam using samtools view -c test1.bam but, it throws the error

[bam_header_read] EOF marker is absent. The input is probably truncated.
[bam_header_read] invalid BAM binary header (this is not a BAM file).
[main_samview] fail to read the header from "test1.bam"


Any help??

[dariober]

Hi,

I tried to do this using the following command

samtools view -L sample.bed test.bam > test1.bam

and surprisingly the output files is thrice in size than the original and i tried to count the number of reads on the output bam using samtools view -c test1.bam but, it throws the error

[bam_header_read] EOF marker is absent. The input is probably truncated.
[bam_header_read] invalid BAM binary header (this is not a BAM file).
[main_samview] fail to read the header from "test1.bam"


Any help??

Hello,
By default, samtools view expect bam as input and produces sam as output. This should explain why you get a very large output (uncompressed sam) and a complain about BAM binary header.
To fix it use the -b option. This should work:

samtools view -b -L sample.bed test.bam > test1.bam
samtools view -c test1.bam

Hope this helps
Dario

[wanfahmi]

Hi,

How could we extract reads for specific locus / gene from accepted_hits.bam file? Is that possible? Could we get the fastq / fasta from it? Thank you

[ersgupta]

try out "intersectBed" in bedtools...

http://code.google.com/p/bedtools/wiki/Usage#intersectBed