Question: Find reference base for a position in a given read from a BAM file
0
gravatar for Andrew
4.4 years ago by
Andrew50
Canada
Andrew50 wrote:

If I take a random read from a BAM file, and want to know what the reference base would be at position 1 of that read, is there a way to do this?  I was thinking I could do something with the POS  column of the BAM file and then add to it the position of the base in the read I am looking at to get the position on the reference sequence.

Maybe a different way to say this is does there exist a way to get the nth base from a specific chromosome in a reference sequence (Fasta).

 

bam • 3.1k views
ADD COMMENTlink modified 4.4 years ago by Matt Shirley9.2k • written 4.4 years ago by Andrew50
3
gravatar for Matt Shirley
4.4 years ago by
Matt Shirley9.2k
Cambridge, MA
Matt Shirley9.2k wrote:

There is another way to get the reference base from just the BAM file, by using the MD tag: https://github.com/vsbuffalo/devnotes/wiki/The-MD-Tag-in-BAM-Files. I've implemented this as the .parse_md method in the simplesam Sam class.  

ADD COMMENTlink written 4.4 years ago by Matt Shirley9.2k

Thanks for the explanation of the MD flag and the detailed examples Matt :)
It blows my mind that the MD tag doesn't include a way to encode insertions... but I cant say it surprises me :P
Your parser, if his data contains the MD flag, is the way to go I think!

ADD REPLYlink written 4.4 years ago by John12k
1
gravatar for Andrew
4.4 years ago by
Andrew50
Canada
Andrew50 wrote:

Found the solution

http://seqanswers.com/forums/showthread.php?t=17315

In case anyone finds this in the future just do the following

samtools faidx <file.fa> chr1:500-500

This will get you the 500th base in the given fasta file.

ADD COMMENTlink written 4.4 years ago by Andrew50
2

Take care that if the read aligns with indels then the n-th position in the read will not correspond to the the POS+n position in the reference.

ADD REPLYlink written 4.4 years ago by dariober10k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1322 users visited in the last hour