bamnostic/utils.py
from __future__ import print_function
from __future__ import absolute_import
from __future__ import division
"""Module utilities and constants used throughout bamnostic
Copyright (c) 2018, Marcus D. Sherman
This code is part of the bamnostic distribution and governed by its
license. Please see the LICENSE file that should have been included
as part of this package.
Some methods are modified versions of their counterparts
within the BioPython.bgzf module. Below is the Copyright and licensing
for those parts.
Copyright (c) 2010-2015 by Peter Cock.
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
@author: "Marcus D. Sherman"
@copyright: "Copyright 2018, University of Michigan, Mills Lab
@email: "mdsherman<at>betteridiot<dot>tech"
"""
import struct
from collections import OrderedDict, namedtuple
import collections
# Python 2 doesn't put abstract base classes in the same spot as Python 3
import sys
_PY_VERSION = sys.version_info
_is_pypy = '__pypy__' in sys.builtin_module_names
if _is_pypy:
import __pypy__
if _PY_VERSION[0] == 2:
from collections import Sequence
from collections import OrderedDict as dict
else:
from collections.abc import Sequence
import numbers
import warnings
import re
def format_warnings(message, category, filename, lineno, file=None, line=None):
r"""Sets STDOUT warnings
Args:
message: the unformatted warning message being reported
category (str): the level of warning (handled by `warnings` module)
filename (str): filename for logging purposes (defaults to STDOUT)
lineno (int): where the error occurred.
Returns:
formatted warning string
"""
return ' {}:{}: {}:{}'.format(filename, lineno, category.__name__, message)
warnings.formatwarning = format_warnings
# pre-compiled structures to reduce iterative unpacking
unpack_int32 = struct.Struct('<i').unpack
unpack_int32L = struct.Struct('<l').unpack
# Helper class for performant named indexing of region of interests
class Roi(object):
r"""Small __slots__ class for region of interest parsing"""
__slots__ = ['__contig', 'start', 'stop', '__tid']
def __init__(self, contig=None, start=None, stop=None, tid=None):
r""" Initialize the class
Args:
contig (str): string representation of chromosome/contig of interest
start (int): starting base position of region of interest
stop (int): ending base position of region of interest
tid (int): position of reference within the BAM header
"""
self.__contig, self.start, self.stop, self.__tid = contig, start, stop, tid
@property
def contig(self):
return self.__contig
@contig.setter
def contig(self, ref_name):
self.__contig = ref_name
@property
def tid(self):
return self.__tid
@tid.setter
def tid(self, value):
self.__tid = value
def __repr__(self):
if self.tid is not None and not self.contig:
return 'Roi(tid: {}, start: {}, stop: {})'.format(self.tid, self.start, self.stop)
elif self.contig and self.tid is None:
return 'Roi(contig: {}, start: {}, stop: {})'.format(self.contig, self.start, self.stop)
else:
return 'Roi(tid: {}, contig: {}, start: {}, stop: {})'.format(self.tid, self.contig, self.start, self.stop)
def __str__(self):
return self.__repr__()
def flag_decode(flag_code):
r"""Simple read alignment flag decoder
Every read within a BAM file ought to have an associated flag code. Theses
flags are used for read filtering and QC. The flags are described below.
Additionally, they can be found `here <https://samtools.github.io/hts-specs/SAMv1.pdf>`_
Any given read's flag is determined by the *or* (`|`) operand of all appropriate bit flags.
Args:
flag_code (int): either a standalone integer/bit flag or the read object itself
Returns:
(:obj:`list` of :obj:`tuple`): list of flag and flag description tuples.
Raises:
ValueError: if provided flag is not a valid entry
Example:
If a flag is 516 it is comprised of flag 4 and flag 512
>>> flag_decode(516)
[(4, 'read unmapped'), (512, 'QC fail')]
==== ===== ==================================================================
Int Bit Description
==== ===== ==================================================================
1 0x1 Template having multiple segments in sequencing
2 0x2 Each segment properly aligned according to the aligner
4 0x4 Segment unmapped
8 0x8 Next segment in the template unmapped
16 0x10 SEQ being reverse complemented
32 0x20 SEQ of the next segment in the template being reverse complemented
64 0x40 The first segment in the template
128 0x80 The last segment in the template
256 0x100 Secondary alignment
512 0x200 Not passing filters, such as platform/vendor quality controls
1024 0x400 PCR or optical duplicate
2048 0x800 Supplementary alignment
==== ===== ==================================================================
"""
flags = {0x1: 'read paired', 0x2: 'read mapped in proper pair',
0x4: 'read unmapped', 0x8: 'mate unmapped',
0x10: 'read reverse strand', 0x20: 'mate reverse strand',
0x40: 'first in pair', 0x80: 'second in pair',
0x100: 'secondary alignment', 0x200: 'QC fail',
0x400: 'PCR or optical duplicate', 0x800: 'supplementary alignment'}
if isinstance(flag_code, numbers.Integral):
code = flag_code
else:
code = flag_code.flag
if not isinstance(code, numbers.Integral):
raise ValueError('Provided flag is not a valid entry')
return sorted([(key, flags[key]) for key in flags if key & code])
def yes_no(message):
""" Simple prompt parser"""
yes = set(['yes', 'ye', 'y', ''])
no = set(['no', 'n'])
print(message)
while True:
answer = input('Would you like to continue? [y/n] ').lower()
if answer in yes:
return True
elif answer in no:
return False
else:
print('Please answer "Yes" or "No"')
def filter_read(read, read_callback='all'):
# Accept all reads
if read_callback == 'nofilter':
return True
# check the read flags against filter criteria
elif read_callback == 'all':
return not read.flag & 0x704 # hex for filter criteria flag bits
# custom filter
elif callable(read_callback):
return read_callback(read)
else:
raise RuntimeError('read_callback should be "all", "nofilter", or a custom function that returns a boolean')
def _parse_sam_region(region):
""" Splits and casts SAM-formatted regions"""
sam_region = ':'.join(region.split()).replace('-', ':').split(':')
# convert start and stop to integers
for i, arg in enumerate(sam_region[1:]):
sam_region[i + 1] = int(arg)
return sam_region
def _handle_split_region(split_roi, until_eof=False):
""" Checks format against `until_eof` and creates the Roi object
Args:
split_roi (:py:obj:`list` or :py::obj:`tuple`): the contig, start, and stop information.
until_eof (bool): whether or not to allow access to end of reference or file (whichever is first)
Returns:
(:py:obj:`bamnostic.utils.Roi`): region of interest formatted as an object with named attributes.
Raises:
ValueError: if `until_eof` is not set and region is open-ended or improper region format.
"""
split_roi = list(split_roi)
# make sure the user didn't put multiple positional arguments
if 1 <= len(split_roi) <= 3:
# if the user gives an integer description of chromosome, convert to string
if isinstance(split_roi[0], (str, int)):
split_roi[0] = str(split_roi[0])
if None in split_roi[1:]:
# make sure the user wants to continue if they have used an open-ended region
if not until_eof:
raise ValueError('Open-ended region while `until_eof` is set to False')
return Roi(*split_roi)
else:
raise ValueError('improper region format')
def parse_region(contig=None, start=None, stop=None, region=None,
tid=None, reference=None, end=None, until_eof=False):
r""" Parses region information from all user set parameters.
The main goal of this function is to handle the many different ways a user
can put in genomic region data. One way is through keyword arguments. This
is the most straight forward. However, due to Pysam's API, there are multiple
synonyms for reference/contig or stop/end. Additionally, if the user knows
the refID/TID of the reference they are interested in, they can input that way
as well.
The second form a submission can take is through positional arguments. Just like
keyword, but ordered such that it makes up a genomic region of interest.
Note:
Positional arguments make utilizing the `tid` parameter difficult since it is
the 5th argument of the function signature.
The third form a submission can take is through using a SAM-formatted string. An
example SAM-formatted string looks like 'chr1:10-100'. As many users also copy &
paste directly from tab-delimited files, such as BED files, a SAM-formatted string
can take the form of 'chr1\\t10\\t100' where '\\t' indicates a tab space.
Lastly, the `until_eof` switch allows users to take all items from their desired
start position (be it the whole reference or a specific spot on the reference). Setting
this to `True` (default: `False`) will pull all reads to the end of the reference or file,
whichever is first.
Args:
contig (str): name of reference/contig
start (int): start position of region of interest (0-based)
stop (int): stop position of region of interest (0-based)
region (str): SAM region formatted string. Accepts tab-delimited values as well
tid (int): the refID or target id of a reference/contig
until_eof (bool): iterate until end of file (default: False)
Returns:
query (:py:class:`bamnostic.utils.Roi`): region of interest formatted as an object with named attributes.
Raises:
ValueError: if two synonym keywords are set, but contradict each other
Examples:
.. code-block:: python
# Keyword-based
>>> parse_region(contig = 'chr1', start = 10, stop = 100)
Roi(contig: chr1, start: 10, stop: 100)
# Using `tid` instead of `contig`
>>> parse_region(tid = 0, start = 10, stop = 100)
Roi(tid: 0, start: 10, stop: 100)
# Positional arguments
>>> parse_region('chr1', 10, 100)
Roi(contig: chr1, start: 10, stop: 100)
# SAM-formatted string (keyword)
>>> parse_region(region = 'chr1:10-100')
Roi(contig: chr1, start: 10, stop: 100)
# SAM-formatted string (positional)
>>> parse_region('chr1:10-100')
Roi(contig: chr1, start: 10, stop: 100)
# Tab-delimited region string
>>> parse_region('chr1\t10\t100')
Roi(contig: chr1, start: 10, stop: 100)
# Contradictory synonyms
>>> parse_region(reference='chr1', contig='chr10', start=10, stop = 100)
Traceback (most recent call last):
...
ValueError: either contig or reference must be set, not both
"""
# Check synonyms for the reference sequence
if contig and reference:
if contig != reference:
raise ValueError('either contig or reference must be set, not both')
else:
contig = contig
elif contig or reference:
contig = contig if contig else reference
# make sure the same thing isn't computed twice
if type(contig) is Roi: # class defined in bamnostic.utils
query = contig
# check for SAM-formatted regions or bed file format
elif region or (contig is not None and (':' in contig or '\t' in contig)):
roi = region if region else contig
query = _handle_split_region(_parse_sam_region(roi), until_eof=until_eof)
else:
if tid and not contig:
contig = None
if (stop and end) and (stop != end):
raise ValueError('either stop or end must be set, not both')
else:
stop = stop if stop else end
query = _handle_split_region((contig, start, stop), until_eof=until_eof)
query.tid = tid
return query
def unpack(fmt, _io, is_array=False):
"""Utility function for unpacking binary data from file object or byte
stream.
The only difference between this method and `struct.unpack` is that
`unpack` dynamically determines the size needed to read in based on
the format string. Additionally, it can process a file object or byte
stream and implement a read or slice (respectively). Mainly, this is a
quality of life function.
Args:
fmt (str): the string format of the binary data to be unpacked
_io: built-in binary format reader (default: io.BufferedRandom)
Returns:
unpacked contents from _io based on fmt string
"""
size = struct.calcsize(fmt)
try:
# if it is byte object
out = struct.unpack(fmt, _io)
except:
# if it is a file object
out = struct.unpack(fmt, _io.read(size))
if len(out) > 1 or is_array:
return out
else:
return out[0]
def make_virtual_offset(block_start_offset, within_block_offset):
"""Compute a BGZF virtual offset from block start and within block offsets.
The BAM indexing scheme records read positions using a 64 bit
'virtual offset', comprising in C terms:
block_start_offset << 16 | within_block_offset
Here block_start_offset is the file offset of the BGZF block
start (unsigned integer using up to 64-16 = 48 bits), and
within_block_offset within the (decompressed) block (unsigned
16 bit integer).
>>> make_virtual_offset(0, 0)
0
>>> make_virtual_offset(0, 1)
1
>>> make_virtual_offset(0, 2**16 - 1)
65535
>>> make_virtual_offset(0, 2**16)
Traceback (most recent call last):
...
ValueError: Require 0 <= within_block_offset < 2**16, got 65536
"""
if within_block_offset < 0 or within_block_offset >= 65536:
raise ValueError("Require 0 <= within_block_offset < 2**16, got %i" %
within_block_offset)
if block_start_offset < 0 or block_start_offset >= 281474976710656:
raise ValueError("Require 0 <= block_start_offset < 2**48, got %i" %
block_start_offset)
return (block_start_offset << 16) | within_block_offset
def split_virtual_offset(virtual_offset):
"""Divides a 64-bit BGZF virtual offset into block start & within block offsets.
>>> (100000, 0) == split_virtual_offset(6553600000)
True
>>> (100000, 10) == split_virtual_offset(6553600010)
True
"""
coffset = virtual_offset >> 16
uoffset = virtual_offset ^ (coffset << 16)
return coffset, uoffset
class LruDict(OrderedDict):
"""Simple least recently used (LRU) based dictionary that caches a given
number of items.
"""
def __init__(self, *args, **kwargs):
""" Initialize the dictionary based on collections.OrderedDict. This
is built of the basic `OrderedDict`. The major difference in instantiation
is the usage of the `max_cache` argument. This sets the dictionary size
to be used.
Args:
items (iterable): an iterable object of key/value pairs
max_cache (int): integer divisible by 2 to set max size of dictionary
"""
try:
self.max_cache = kwargs.pop('max_cache', 128)
mode = kwargs.pop('mode', 'fifo')
if mode == 'fifo':
self.mode = False
else:
self.mode = True
except AttributeError:
self.max_cache = 128
self.mode = False # FIFO
super(LruDict, self).__init__(*args, **kwargs)
self.__map = {}
if _is_pypy:
self.move_to_end = self._pypy_move_to_end
elif _PY_VERSION[:2] <= (3,2):
self.move_to_end = self._py27_move_to_end
else:
self.move_to_end = super(LruDict, self).move_to_end
self.cull()
def get(self, key):
""" Basic getter that renews LRU status upon inspection
Args:
key (str): immutable dictionary key
"""
value = OrderedDict.__getitem__(self, key)
self.move_to_end(key)
return value
def _pypy_move_to_end(self, key, last=True):
__pypy__.move_to_end(self, key, last)
def _py27_move_to_end(self, key, last=True):
"""Move an existing element to the end (or beginning if last is false).
Raise KeyError if the element does not exist.
"""
link = self._OrderedDict__map[key]
link_prev, link_next, _ = link
soft_link = link_next[0]
link_prev[1] = link_next
link_next[0] = link_prev
root = self._OrderedDict__root
if last:
last = root[0]
link[0] = last
link[1] = root
root[0] = soft_link
last[1] = link
else:
first = root[1]
link[0] = root
link[1] = first
first[0] = soft_link
root[1] = link
def update(self, others):
""" Same as a regular `dict.update`, however, since pypy's `dict.update`
doesn't go through `dict.__setitem__`, this is used to ensure it does
Args:
others (iterable): a dictionary or iterable containing key/value pairs
"""
if type(others) == dict:
for k,v in others.items():
self.__setitem__(k,v)
elif type(others) in [list, tuple]:
for k,v in others:
self.__setitem__(k,v)
else:
raise ValueError('Iteratable/dict must be in key, value pairs')
def cull(self):
""" Main utility function for pruning the LruDict
If the length of the LruDict is more than `max_cache`, it removes the LRU item
"""
if self.max_cache:
overflow = max(0, len(self) - self.max_cache)
if overflow > 0:
for _ in range(overflow):
self.popitem(last=self.mode)
def __setitem__(self, key, value):
"""Basic setter that adds new item to dictionary, and then performs cull()
to ensure max_cache has not been violated.
Args:
key (str): immutable dictionary key
value (any): any dictionary value
"""
OrderedDict.__setitem__(self, key, value)
self.cull()
# The BAM format uses byte encoding to compress alignment data. One such
# compression is how operations are stored: they are stored and an
# array of integers. These integers are mapped to their respective
# operation identifier. Below is the mapping utility.
_CIGAR_OPS = {'M': ('BAM_CMATCH', 0),
'I': ('BAM_CINS', 1),
'D': ('BAM_CDEL', 2),
'N': ('BAM_CREF_SKIP', 3),
'S': ('BAM_CSOFT_CLIP', 4),
'H': ('BAM_CHARD_CLIP', 5),
'P': ('BAM_CPAD', 6),
'=': ('BAM_CEQUAL', 7),
'X': ('BAM_CDIFF', 8),
'B': ('BAM_CBACK', 9)}
def parse_cigar(cigar_str):
"""Parses a CIGAR string and turns it into a list of tuples
Args:
cigar_str (str): the CIGAR string as shown in SAM entry
Returns:
cigar_array (list): list of tuples of CIGAR operations (by id) and number of operations
Raises:
ValueError: if CIGAR operation is invalid
Examples:
>>> parse_cigar('3M1I3M1D5M') # doctest: +ELLIPSIS, +NORMALIZE_WHITESPACE
[(('BAM_CMATCH', 0), 3), ..., (('BAM_CMATCH', 0), 5)]
"""
cigar_array = []
for cigar_op in re.finditer(r'(?P<n_op>\d+)(?P<op>\w)', cigar_str):
op_dict = cigar_op.groupdict()
n_ops = int(op_dict['n_op'])
op = _CIGAR_OPS.get(op_dict['op'], -1)
if op == -1:
raise ValueError('Invalid CIGAR operation ({}).'.format(op_dict['op']))
cigar_array.append((op, n_ops))
return cigar_array
def check_cigar_arg(cigar):
""" Checks to make sure CIGAR arugment is valid.
Args:
argument (str or :py:obj:`list`): CIGAR string (pre-formatted or raw)
Returns:
(:py:obj:`list`): CIGAR re-formatted as a list
Raises:
ValueError: if CIGAR is not a string or pre-formatted list
"""
if type(cigar) == str:
cigar = parse_cigar(cigar)
elif type(cigar) == list:
pass
else:
raise ValueError('CIGAR must be string or list of tuples of cigar operations (by ID) and number of operations')
return cigar
def cigar_changes(seq, cigar):
"""Recreates the reference sequence to the extent that the CIGAR string can
represent.
Args:
seq (str): aligned segment sequence
cigar (list): list of tuples of cigar operations (by id) and number of operations
Returns:
cigar_formatted_ref (str): a version of the aligned segment's reference \
sequence given the changes reflected in the cigar string
Raises:
ValueError: if CIGAR operation is invalid
Examples:
>>> cigar_changes('ACTAGAATGGCT', '3M1I3M1D5M')
'ACTGAATGGCT'
>>> cigar_changes('ACTAGAATGGCT', '3V1I3M1D5M')
Traceback (most recent call last):
...
ValueError: Invalid CIGAR operation (V).
"""
cigar = check_cigar_arg(cigar)
cigar_formatted_ref = ''
last_cigar_pos = 0
for op, n_ops in cigar:
op_id = op[1]
if op_id in {0, 7, 8}: # matches (uses both sequence match & mismatch)
cigar_formatted_ref += seq[last_cigar_pos:last_cigar_pos + n_ops]
last_cigar_pos += n_ops
elif op_id in {1, 4}: # insertion or clips
last_cigar_pos += n_ops
elif op_id == 3: # intron or large gaps
cigar_formatted_ref += 'N' * n_ops
elif op_id in {2, 5}:
pass
else:
raise ValueError('Invalid CIGAR string: {}'.format(op))
return cigar_formatted_ref
def md_changes(seq, md_tag):
"""Recreates the reference sequence of a given alignment to the extent that the
MD tag can represent.
Note:
Used in conjunction with `cigar_changes` to recreate the
complete reference sequence
Args:
seq (str): aligned segment sequence
md_tag (str): MD tag for associated sequence
Returns:
ref_seq (str): a version of the aligned segment's reference sequence given \
the changes reflected in the MD tag
Raises:
ValueError: if MD tag is None
Example:
>>> md_changes('CTTATATTGGCCTT', '3C4AT4')
'CTTCTATTATCCTT'
"""
if md_tag is None:
raise ValueError('No MD tag found or given for sequence')
ref_seq = ''
last_md_pos = 0
for mo in re.finditer(r'(?P<matches>\d+)|(?P<del>\^\w+?(?=\d))|(?P<sub>\w)', md_tag):
mo_group_dict = mo.groupdict()
if mo_group_dict['matches'] is not None:
matches = int(mo_group_dict['matches'])
ref_seq += seq[last_md_pos:last_md_pos + matches]
last_md_pos += matches
elif mo_group_dict['del'] is not None:
deletion = mo_group_dict['del']
ref_seq += deletion[1:]
elif mo_group_dict['sub'] is not None:
substitution = mo_group_dict['sub']
ref_seq += substitution
last_md_pos += 1
else:
pass
return ref_seq
def ref_gen(seq, cigar_string, md_tag):
"""Recreates the reference sequence associated with the given segment.
Uses the CIGAR string and MD tag to recreate the reference sequence associated
with the aligned segment. This is done without the need for looking up
the reference genome. Example reads, MD tags, and CIGAR strings taken from
`David Tang's Blog`_.
.. _David Tang's Blog: https://davetang.org/muse/2011/01/28/perl-and-sam/
Returns:
(str): generated reference sequence
Raises:
KeyError: if read does not contain MD tag
Examples:
# Only mismatches
>>> seq = 'CGATACGGGGACATCCGGCCTGCTCCTTCTCACATG'
>>> cigar = '36M'
>>> md = '1A0C0C0C1T0C0T27'
>>> ref_gen(seq, cigar, md)
'CACCCCTCTGACATCCGGCCTGCTCCTTCTCACATG'
# Insertions and mismatches
>>> seq = 'GAGACGGGGTGACATCCGGCCTGCTCCTTCTCACAT'
>>> cigar = '6M1I29M'
>>> md = '0C1C0C1C0T0C27'
>>> ref_gen(seq, cigar, md)
'CACCCCTCTGACATCCGGCCTGCTCCTTCTCACAT'
# Deletion and mismatches
>>> seq = 'AGTGATGGGGGGGTTCCAGGTGGAGACGAGGACTCC'
>>> cigar = '9M9D27M'
>>> md = '2G0A5^ATGATGTCA27'
>>> ref_gen(seq, cigar, md)
'AGGAATGGGATGATGTCAGGGGTTCCAGGTGGAGACGAGGACTCC'
# Insertion, deletion, and mismatches
>>> seq = 'AGTGATGGGAGGATGTCTCGTCTGTGAGTTACAGCA'
>>> cigar = '2M1I7M6D26M'
>>> md = '3C3T1^GCTCAG26'
>>> ref_gen(seq, cigar, md)
'AGGCTGGTAGCTCAGGGATGTCTCGTCTGTGAGTTACAGCA'
"""
return md_changes(cigar_changes(seq, cigar_string), md_tag)
def cigar_alignment(seq=None, cigar=None, start_pos = 0, qualities=None, base_qual_thresh=0, query=False):
"""Use the CIGAR to filter out all unaligned data bases
Any clipping results in the removal of those bases. If an insertion is seen in
the CIGAR, those bases are removed from the sequence. If a deletion is seen in
the CIGAR, those bases are padded with a period ('.') symbol.
Args:
seq (str): string sequence of the aligned segment.
cigar (str): the cigar string or `cigartuple` of the aligned segment.
start_pos (int): the first aligned position of the read
qualities (:py:obj:`array.array`): base quality array from read
Yields:
(:py:obj:`tuple` of :py:obj:`str` and :py:obj:`int`): nucleotide base and index position of that base relative to reference
Example:
>>> seq = 'AGTGATGGGAGGATGTCTCGTCTGTGAGTTACAGCA'
>>> cigar = '2M1I7M6D26M'
>>> start_position = 95
>>> c_a = cigar_alignment(seq, cigar, start_position)
>>> next(c_a)
('A', 95)
"""
cigar = check_cigar_arg(cigar)
cigar_aligned = ''
algn_seg = {}
last_cigar_pos = 0
for op, n_ops in cigar:
op_id = op if type(op) is int else op[1]
if op_id == 5: # BAM_CHARD_CLIP: skip hard clip CIGAR ops
pass
elif op_id in {1, 4}: # BAM_CINS or BAM_CSOFT_CLIP: remove from sequence
if query and op_id == 1:
seg_seq = seq[last_cigar_pos:last_cigar_pos + n_ops]
if qualities is not None:
seg_qual = qualities[last_cigar_pos:last_cigar_pos + n_ops]
for index, base in enumerate(seg_seq):
if qualities is not None:
if seg_qual[index] >= base_qual_thresh:
yield base, start_pos
else:
yield base, start_pos
last_cigar_pos += n_ops
elif op_id == 3: # BAM_CREF_SKIP: intron or large gaps
start_pos += n_ops
elif op_id == 2: # BAM_CDEL: pad for deletions
start_pos += n_ops
elif op_id in {0, 7, 8}: # matches (uses both sequence match & mismatch)
seg_seq = seq[last_cigar_pos:last_cigar_pos + n_ops]
if qualities is not None:
seg_qual = qualities[last_cigar_pos:last_cigar_pos + n_ops]
for index, base in enumerate(seg_seq):
if qualities is not None:
if seg_qual[index] >= base_qual_thresh:
yield base, start_pos
else:
yield base, start_pos
start_pos += 1
last_cigar_pos += n_ops
else:
raise ValueError('Invalid CIGAR string: {}'.format(op))