Samtools view bam parse region issues

[rkizen]

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!

[lh3]

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.

[rkizen]

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.

[dpryan]

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.

[rkizen]

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.

[dpryan]

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.

[lh3]

@dpryan your reheader is better in the long run. Haven't thought about that...