Looking for easy way to find RIGHT-most end of alignment in bam

[jmartin]

Is there a way to report the right-most end of an alignment in a bam file? The position given is the left-most part of the alignment, but I need to know where the right-end falls in the reference. I know this can be done with the CIGAR string, but I'm concerned that I don't fully understand how the CIGAR string works with regards to hard & soft clipping, as well as what its reporting with regards to chimeric alignments. I think it might be easy to get this wrong if I'm not careful, and if there is an already existing, robust & convenient tool or method I'd prefer using that.

[Brian Bushnell]

The easiest way to do this is to map with BBMap, using the "stoptag" flag, which will make it write an extra tag to each line prefixed by YS:i: which gives the read's stop location. I don't currently have a way to do that for existing sam/bam files, but I may incorporate it into reformat.

This is the code I use:

public static int calcCigarLength(String cigar, boolean includeHardClip){
	if(cigar==null){return 0;}
	int len=0;
	int current=0;
	for(int i=0; i<cigar.length(); i++){
		char c=cigar.charAt(i);
		if(Character.isDigit(c)){
			current=(current*10)+(c-'0');
		}else{
			if(c=='M' || c=='=' || c=='X' || c=='D' || c=='N' || c=='S'){
				len+=current;
			}else if (c=='H'){ 
				if(includeHardClip){len+=current;}
			}else if(c=='I'){
				//do nothing
			}else if(c=='P'){
				throw new RuntimeException("Unhandled cigar symbol: "+c+"\n"+cigar+"\n");
				//'P' is currently poorly defined
			}else{
				throw new RuntimeException("Unhandled cigar symbol: "+c+"\n"+cigar+"\n");
			}
			current=0;
		}
	}
	return len;
}

But, as you said, it gets more complicated when you have clipped reads:

public static String makeStopTag(int pos, int seqLength, String cigar, boolean perfect){
	return "YS:i:"+(pos+((cigar==null || perfect) ? seqLength : -countLeadingClip(cigar, false)+calcCigarLength(cigar, false))-1);
}

public static int countLeadingClip(String cigar, boolean includeHardClip){
	if(cigar==null){return 0;}
	int len=0;
	int current=0;
	for(int i=0; i<cigar.length(); i++){
		char c=cigar.charAt(i);
		if(Character.isLetter(c) || c=='='){
			if((includeHardClip && c=='H') || (SUBTRACT_LEADING_SOFT_CLIP && c=='S')){
				len+=current;
			}else{
				break;
			}
			current=0;
		}else{
			current=(current*10)+(c-'0');
		}
	}
	return len;
}

Different aligners calculate pos differently for clipped reads which is why I have flags for counting or not counting hard- and soft-clipped bases.

[dariober]

Is there a way to report the right-most end of an alignment in a bam file?

This script below based on python pysam should do it. It adds a YS tag with the end position of the read (mocking Brian's answer).

Save the script as addReadEndToBam.py (or whatever). It will print to stdout in BAM format.

Example usage:

addReadEndToBam.py in.bam > out.bam

# Example output alignment:
M00886:29:000000000-A95CG:1:2103:5504:6001	89	LmjF.01	4	0	34M	=	4	0	CCCTAACCCTAACCTTGACCCTAACCCTATCCCT	FB0AA0AAB1BBBA11GFCGFAB1>B1>1>>>>>	XT:A:R	NM:i:2	SM:i:0	AM:i:0	X0:i:5	X1:i:0	XM:i:2	XO:i:0	XG:i:0	MD:Z:14C14A4	YS:i:37

This is based on pysam 0.8.1:

#!/usr/bin/env python

import pysam
import sys

insam= sys.argv[1]
samfile = pysam.AlignmentFile(insam, "rb")
outfile = pysam.AlignmentFile("-", "wb", template=samfile)

for aln in samfile:
    ys= aln.reference_end
    if not ys:
        ys= -1
    aln.setTag('YS', ys)
    outfile.write(aln)

samfile.close()
outfile.close()
sys.exit()