7
$\begingroup$

Is there a standard method in Python to extract a CIGAR string from the BAM?

There are great libraries which parse the CIGAR, e.g. https://pypi.python.org/pypi/cigar/0.1

>>> c = Cigar('10M20S10M')
>>> c.mask_left(10).cigar
'10S20S10M'
>>> c.mask_left(9).cigar
'9S1M20S10M'
>>> Cigar('10S').mask_left(10).cigar
'10S'
>>> Cigar('10H').mask_left(10).cigar
'10H'
>>> Cigar('10H').mask_left(11).cigar
'10H'
>>> Cigar('10H').mask_left(9).cigar
'10H'

It looks like pysam already parses the cigar string, if I'm not mistaken:

import pysam

bam = 'myfile.bam'
bamfile = pysam.AlignmentFile(bam, 'rb')

for read in bamfile:
    if not read.is_unmapped:
        cigar = read.cigar
        print(cigar)

This outputs lists of tuples of the parse output:

[(5, 61), (0, 30), (5, 198)]
[(4, 11), (0, 30), (4, 248)]
[(4, 11), (0, 30), (4, 248)]
....

How could I simply output the CIGAR from read? I would prefer to create a column in a pandas DataFrame

$\endgroup$
3
$\begingroup$

If you really do just want the cigar string then it's read.cigarstring. However, I'm not sure what you're trying to gain with the cigar package from Brent. Unless you want to get the string with the masking changed then the tuple you get from pysam is the same as what you get from cigar (with the exception of the numeric operations instead of character operations).

$\endgroup$
  • $\begingroup$ Oh, the API on the website had me confused with cigar string. Perhaps I could use the parsed cigar from pysam, but I'm confused: How do I parse the list of tuples in order count the total number of insertions, deletions, matches (and possible other information) from this? I apologise if this is documented, but the pysam manual is confusing imho. $\endgroup$ – ShanZhengYang Mar 9 '18 at 19:42
  • $\begingroup$ Should I read the above, e.g. the first line [(5, 61), (0, 30), (5, 198)] as 61 hard-clippings, 30 matches, and 198 hard clippings? $\endgroup$ – ShanZhengYang Mar 9 '18 at 19:52
  • $\begingroup$ Yes, that is the correct interpretation of that list of tuples, which is documented in the API $\endgroup$ – Devon Ryan Mar 9 '18 at 20:19
  • $\begingroup$ I see; this question is largely my misunderstanding of the pysam website. I'm happy to delete this question, unless it's helpful for the community. $\endgroup$ – ShanZhengYang Mar 10 '18 at 18:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.