Last active
October 23, 2015 06:35
-
-
Save andreas-wilm/236be00ad1dc014ecf1b to your computer and use it in GitHub Desktop.
Example on how to reconstruct alignment from pysam aligned segment
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# assuming: | |
# - refsq is set and always the same | |
# - bamfh is an open pysam samfile | |
for r in bamfh: | |
if r.is_unmapped: | |
continue | |
for (qapos, rpos) in r.get_aligned_pairs(): | |
# qapos is the aligned index, i.e. this ignores clipping. add that | |
if qapos is None: | |
qpos = None | |
else: | |
qpos = qapos + r.query_alignment_start | |
# qpos and rpos now safe to use, but might be None (indel) | |
rbase = qbase = "-" | |
if rpos is not None: | |
rbase = refsq[rpos].upper() | |
if qpos is not None: | |
qbase = r.seq[qpos].upper() | |
# now you can use rbase and qbase | |
# ... |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment