SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
samtools -view region (overlap?) question dietmar13 Bioinformatics 5 04-22-2013 11:44 AM
samtools: parse error in SAM to BAM conversion chrisW Bioinformatics 12 01-14-2013 07:16 PM
View BAM file syintel87 Bioinformatics 3 01-06-2013 02:40 PM
How to use the Bioperl to parse the parse flat file of UniProtKB database? bewlib Bioinformatics 1 11-29-2012 05:30 PM
Samtools view region ftp DavyK Bioinformatics 1 09-13-2012 04:05 AM

Reply
 
Thread Tools
Old 07-25-2013, 02:46 AM   #1
rkizen
Member
 
Location: Edinburgh

Join Date: Jun 2013
Posts: 20
Default Samtools view bam parse region issues

I have mapped reads from illumina sequencing onto a reference genome using bwa aln and I wanted to get some regions of the bam files using samtools view:

samtools view -o test_region.sam merged.sort.bam 1:100-1000

however I kept getting the error:

[bam_parse_region] fail to determine the sequence name.

I have realized that I used an NCBI version of the reference genome which resulted in the reference sequence name in the bam file to be that of the NCBI sequence names. For example, here is bit of the bam file:

EBRI093151_0050:4:117:5648:11807#0 99 gi|358485511|ref|NC_006088.3| 66454 60 101M = 66940 587 CTCCAAGATCATCCAGTCCACCCATCCACCCACCACCAATGTAACCCCATTAAAACACGTCCCTCAGTACCACATCTAAATGTTTCTTGAGCACTGGCAGG fffffdffffefffffefeffffeffffffffffffffcdefedffdfeffdfaeffffdeffffdedfffeffe`f`ccdfadfffdaYd^cd`bbd`dc XT:A:U NM:i:0 SM:i:37 AM:i:37 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:101


Notice how instead of "chr1" it has "gi|358485511|ref|NC_006088.3|" as the sequence name. I believe this is the cause of the parsing issue. However, now I cannot seem to figure out how to pull a region from the bam file. I tried using:

samtools view -o test_region.sam merged.sort.bam NC_006088.3:1-1000

But I get the same error as before:
[bam_parse_region] fail to determine the sequence name.

If I put this instead:

samtools view -o test_region.sam merged.sort.bam NC_006088.3

It seems to do as one would expect and pull all read mappings for that sequence region, however, I would like to pull specific coordinates from each chromosome.

Does anybody have some ideas of how I could do this? I was thinking about replacing all the reference sequence names with their generic equivalent (ie "chr1,chr2 etc") however I am not certain that would solve my issue here.

Thank you!
rkizen is offline   Reply With Quote
Old 07-25-2013, 06:20 AM   #2
lh3
Senior Member
 
Location: Boston

Join Date: Feb 2008
Posts: 693
Default

Try "gi|358485511|ref|NC_006088.3|:100-1000", though I don't know if that works.

The best solution is to replace the sequence names with something easier.
lh3 is offline   Reply With Quote
Old 07-25-2013, 06:29 AM   #3
rkizen
Member
 
Location: Edinburgh

Join Date: Jun 2013
Posts: 20
Default

Thanks lh3! That works. I should change the sequence names but I needed this for something a bit urgent and there are a lot of sequence names which would require a bit of time ensuring that a name changing script was not corrupting the file in any way. Thank you again and I am kicking myself for not using the quotation marks to get past the pipes in the sequence names.
rkizen is offline   Reply With Quote
Old 07-25-2013, 06:45 AM   #4
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,476
Default

The reheader script can simply be "samtools reheader"

Do a "samtools view -H file.bam > header.sam", edit the header.sam file by replacing the chromosome names (they need to be in the same order as the original, so don't reorder anything!), then use samtools reheader. Unless your alignments are to a large number of contigs, this should be quite quick.
dpryan is offline   Reply With Quote
Old 07-25-2013, 06:59 AM   #5
rkizen
Member
 
Location: Edinburgh

Join Date: Jun 2013
Posts: 20
Default

Thanks dpryan!

I have not used reheader before. Will it replace all the sequence names? or will it just replace the header info?

There are a lot of contigs and there are many mappings which may not have the same contigs represented. I will have to write a script to run "samtools view -H file.bam > header.sam" for each mapping and then replace the sequence names for each of the header sam files and apply those back onto the origin bam file using samtools reheader.

For now I just need some intervals quickly so I will use the method lh3 to get these intervals and then work on converting the headers later.
rkizen is offline   Reply With Quote
Old 07-25-2013, 07:03 AM   #6
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,476
Default

If you have a BAM file, the chromosome/contig names are only stored in the header (there's an ordered list and each read just has a numeric value saying where in that list it's chromosome name can be found). Yeah, if you have a bunch of different mappings, then following what Heng Li suggested might prove faster.
dpryan is offline   Reply With Quote
Old 07-25-2013, 03:51 PM   #7
lh3
Senior Member
 
Location: Boston

Join Date: Feb 2008
Posts: 693
Default

@dpryan your reheader is better in the long run. Haven't thought about that...
lh3 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 12:00 AM.


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