Question: Deciphering Cigar
5
gravatar for Leszek
6.1 years ago by
Leszek3.9k
IIMCB, Poland
Leszek3.9k wrote:

Hi, I need to get 5'-end position of each read in SAM file. For reads aligned in forward direction it's provided in SAM file (POS). But for reads aligned onto reverse strand the one has to figure it out from the cigar string. You cannot simple add length of the read to left-most position, as reads may be soft-clipped or some indel may occure.

I have written simple python function in order to find right most position of aligned read. Could you have a look and comment on that as I'm not sure if it's correct?

import re

#should nicely separate CIGAR entries
cigar_pat = re.compile(r"\d+[MIDNSHP=X]{1}")

def cigar2end( left,cigar ):
  """Return right-most position of aligned read."""
  #store info about each CIGAR category
  counts={ "M":0, #M 0 alignment match (can be a sequence match or mismatch)
           "I":0, #I 1 insertion to the reference
           "D":0, #D 2 deletion from the reference
           "N":0, #N 3 skipped region from the reference
           "S":0, #S 4 soft clipping (clipped sequences present in SEQ)
           "H":0, #H 5 hard clipping (clipped sequences NOT present in SEQ)
           "P":0, #P 6 padding (silent deletion from padded reference)
           "=":0, #= 7 sequence match
           "X":0, #X 8 sequence mismatch
        }
  #split cigar entries
  for centry in cigar_pat.findall(cigar):
    ccount  = int(centry[:-1])
    csymbol = centry[-1]
    counts[csymbol] = ccount
  #get number of aligned 'reference' bases
  aligned = counts["M"] + counts["D"] + counts["N"] #+ counts["="] + counts["X"]
  right   = left + aligned
  return right
ADD COMMENTlink modified 6.1 years ago • written 6.1 years ago by Leszek3.9k
4

You might consider using pysam for this type of thing.

ADD REPLYlink written 6.1 years ago by Sean Davis24k

thanks, got it:)

ADD REPLYlink written 6.1 years ago by Leszek3.9k

I was writing my own codes to process the CIGAR line in order to extract different information that I needed. But now I has been using Pysam and I save a lot of time because it have very useful functions and transform the CIGAR line into a tuple very easy to process.

ADD REPLYlink written 6.1 years ago by Geparada1.3k
5
gravatar for Leszek
6.1 years ago by
Leszek3.9k
IIMCB, Poland
Leszek3.9k wrote:

As Sean suggested I use pysam:

import pysam
samfile = pysam.Samfile( "somefile.bam","rb" )
#process all reads in region
for sam in samfile.fetch( ref,start,end ):
    print sam.pos, sam.aend
ADD COMMENTlink written 6.1 years ago by Leszek3.9k

does your code return the same start/end ?

ADD REPLYlink modified 6.1 years ago • written 6.1 years ago by Pierre Lindenbaum111k

what do you mean by the same start/end? it returns first and last aligned base of the read on the ref

ADD REPLYlink written 6.1 years ago by Leszek3.9k

Thanks for sharing your code. I used them to find start and end position of reads but it asks me to define ref/start/end. Any suggestions?  

ADD REPLYlink written 2.7 years ago by Bionic0

You need to provide the values for "ref", "start", and "end".  The code, as written above, assumes that has been done elsewhere.  

ADD REPLYlink written 2.7 years ago by Sean Davis24k
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: 664 users visited in the last hour