Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • dariober
    replied
    Originally posted by jmartin View Post
    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:

    Code:
    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:

    Code:
    #!/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()
    Last edited by dariober; 03-20-2015, 06:06 AM. Reason: Use reference_end

    Leave a comment:


  • Brian Bushnell
    replied
    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:

    Code:
    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:

    Code:
    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.
    Last edited by Brian Bushnell; 03-18-2015, 01:21 PM.

    Leave a comment:


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

    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.

Latest Articles

Collapse

  • seqadmin
    Multiomics Techniques Advancing Disease Research
    by seqadmin


    New and advanced multiomics tools and technologies have opened new avenues of research and markedly enhanced various disciplines such as disease research and precision medicine1. The practice of merging diverse data from various ‘omes increasingly provides a more holistic understanding of biological systems. As Maddison Masaeli, Co-Founder and CEO at Deepcell, aptly noted, “You can't explain biology in its complex form with one modality.”

    A major leap in the field has
    ...
    02-08-2024, 06:33 AM

ad_right_rmr

Collapse

News

Collapse

Topics Statistics Last Post
Started by seqadmin, 02-21-2024, 08:52 AM
0 responses
30 views
0 likes
Last Post seqadmin  
Started by seqadmin, 02-20-2024, 08:57 AM
0 responses
18 views
0 likes
Last Post seqadmin  
Started by seqadmin, 02-14-2024, 09:19 AM
0 responses
52 views
0 likes
Last Post seqadmin  
Started by seqadmin, 02-12-2024, 03:37 PM
0 responses
451 views
0 likes
Last Post seqadmin  
Working...
X