SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Extract aligned reads from a BAM file above a certain threshold The Snow Bioinformatics 4 07-29-2013 02:02 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-18-2012, 04:26 AM   #1
meher
Member
 
Location: helsinki

Join Date: Jun 2011
Posts: 54
Default Extract reads from multiple regions from bam file

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
meher is offline   Reply With Quote
Old 10-18-2012, 04:38 AM   #2
krobison
Senior Member
 
Location: Boston area

Join Date: Nov 2007
Posts: 747
Default

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:

Code:
samtools view aln.sorted.bam `perl -p -e 'chomp; @a=split(/\t/); print "$a[0]:$a[1]-$a[2]"' sleepy.bed`
krobison is offline   Reply With Quote
Old 10-18-2012, 04:43 AM   #3
meher
Member
 
Location: helsinki

Join Date: Jun 2011
Posts: 54
Default

Quote:
Originally Posted by krobison View Post
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:

Code:
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?
meher is offline   Reply With Quote
Old 10-18-2012, 04:50 AM   #4
dariober
Senior Member
 
Location: Cambridge, UK

Join Date: May 2010
Posts: 310
Default

Quote:
Originally Posted by krobison View Post
It's pretty easy to convert BED format to this using your favorite text mangling language; mine is Perl:

Code:
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
dariober is offline   Reply With Quote
Old 10-18-2012, 06:38 AM   #5
meher
Member
 
Location: helsinki

Join Date: Jun 2011
Posts: 54
Default Extract reads from multiple regions from bam file

Quote:
Originally Posted by dariober View Post
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??
meher is offline   Reply With Quote
Old 10-18-2012, 07:05 AM   #6
dariober
Senior Member
 
Location: Cambridge, UK

Join Date: May 2010
Posts: 310
Default

Quote:
Originally Posted by meher View Post
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:
Code:
samtools view -b -L sample.bed test.bam > test1.bam
samtools view -c test1.bam
Hope this helps
Dario
dariober is offline   Reply With Quote
Old 12-14-2012, 02:33 AM   #7
wanfahmi
Member
 
Location: North Sea

Join Date: Apr 2008
Posts: 34
Default

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
wanfahmi is offline   Reply With Quote
Old 12-14-2012, 03:18 AM   #8
ersgupta
Member
 
Location: India

Join Date: Jun 2011
Posts: 26
Default

try out "intersectBed" in bedtools...


http://code.google.com/p/bedtools/wi...e#intersectBed
ersgupta 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 10:51 AM.


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