Extracting Reference Sequence from a bam File

[andy11]

I'm currently using the bamtools API for a program where I'm trying to extract and print the reference sequence from a .bam file, but can't figure out how I could get to the reference sequence. Is there any way to determine the reference sequence directly from the bam file, or is there a better way to do this?

I've searched but haven't come up with anything.

Thanks for any help,
Andy

[maubp]

You can't

If you read the SAM/BAM file format definition, you'll see they don't actually contain the reference sequences. All the BAM file header contains is a SAM header (optional chunk of embedded text), number of references, their names and lengths.

You would normally have a FASTA file to accompany the SAM/BAM file.

[andy11]

Should've seen that, thanks.

Do you have any advice on how to rapidly extract a subsequence from a fasta file using an index file?

[westerman]

I've never understood the decision not to (optionally) bundle the reference into the BAM/SAM file. It seems that for most secondgen datasets that the space that the reference would take up would be a trivial amount while the amount of headaches avoided by not having to chase down the proper reference would be large. But, alas, I was not consulted. :-)

Getting back to 'andy11's question "how to rapidly extract a subsequence from a fasta file using an index file?" we will need a bit more information. Which program created the index file? I am assuming that you are not talking about the BAM/SAM index since that had nothing to do with the reference file.

[maubp]

I've never understood the decision not to (optionally) bundle the reference into the BAM/SAM file. It seems that for most secondgen datasets that the space that the reference would take up would be a trivial amount while the amount of headaches avoided by not having to chase down the proper reference would be large. But, alas, I was not consulted. :-)

I agree with you that optionally including the reference sequence could be very useful, especially for non-model organisms. The SAM/BAM design is clearly designed more for mapping and re-sequencing than for de novo assembly.

[lh3]

There are many small alignment files far smaller than the compressed human reference genome (e.g. most of 1000 genomes alignments).

Yes, for users, it would always be good to have more options, but for developers, having more options is a burden which has to be evaluated carefully; the adoption of SAM would also be delayed that way.

[bioinfosm]

Should've seen that, thanks.

Do you have any advice on how to rapidly extract a subsequence from a fasta file using an index file?

I think USEQ has a utility to specify bed file with coordinates, and it extracts the fasta sequences from a larger fasta file!