Recovering the original 5' position of reverse strand reads from SAM files

[microphobe]

Hello everyone,

I need to get the precise 5' position of my mapped directional RNA-seq reads. Problem is, for reads that map to the reverse strand, the BAM/SAM format output by my mapping software (TMAP) seems to be the left-most position of the last *aligning* base relative to the reference. In other words, it's the 3'-most base (using the orientation of the read) after clipping non-aligning bases.

I've tried just adding the length of the read onto the position, but this fails for reads with bad 3' ends or several indels in the alignment. I suppose I could write a script to count up the number of "S"s etc. in the cigar, but I am a wet-lab person and the thought of that just makes me sad.

Is there a better way to recover the original 5' position of the seed alignment for reverse reads? Thanks in advance.

[maubp]

The left position is just the POS field, the right position is POS plus the sum of the M/=/X/D/N operations in the CIGAR string. The S operations are irrelevant here.

If you are using a SAM/BAM library to process your data, there will probably be a method/function call already available to do this for you.